AliAODEvent::GetHeader() returns AliVHeader
[u/mrichter/AliRoot.git] / PWGGA / PHOSTasks / PHOS_embedding / AliPHOSEmbedding.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 // Analysis task to perform embedding of the simulated AOD
17 // into real data (ESD)
18 // Output conntainer contain AOD with
19 //   0. Data before embedding
20 //   1. Data+with embedded signal
21 //   2. MC information from signal
22 // Data for embedding are obtained using standard Manager
23 // Signal is read from TChain with aodTree
24 //
25 // Authors: Dmitri Peressounko
26 // Date   : 28.05.2011
27
28 #include "TChain.h"
29 #include "TFile.h"
30 #include "TROOT.h"
31 #include "TClonesArray.h"
32 #include "TH2.h"
33 #include "TGeoManager.h"
34 #include "TGeoGlobalMagField.h"
35
36 #include "AliGeomManager.h"
37 #include "AliGRPObject.h"
38 #include "AliCDBPath.h"
39 #include "AliCDBEntry.h"
40 #include "AliESDEvent.h"
41 #include "AliAODEvent.h"
42 #include "AliAnalysisManager.h"
43 #include "AliVEventHandler.h"
44 #include "AliInputEventHandler.h"
45 #include "AliMultiplicity.h"
46 #include "AliAODMCParticle.h"
47 #include "AliAODMCHeader.h"
48 #include "AliCDBManager.h"
49 #include "AliAODHandler.h"
50 #include "AliMagF.h"
51
52 #include "AliPHOSReconstructor.h"
53 #include "AliPHOSClusterizerv1.h"
54 #include "AliPHOSDigit.h"
55 #include "AliPHOSGeometry.h"
56
57 #include "AliPHOSEmbedding.h"
58
59 ClassImp(AliPHOSEmbedding)
60
61
62 //________________________________________________________________________
63 AliPHOSEmbedding::AliPHOSEmbedding(const char *name) 
64   : AliAnalysisTaskSE(name),
65     fAODChain(0x0),
66     fDigitsTree(0x0),
67     fClustersTree(0x0),
68     fTreeOut(0x0),
69     fDigitsArr(0x0),
70     fEmbeddedClusters(0x0),
71     fEmbeddedCells(0x0),
72     fCellsPHOS(0x0),
73     fClusterizer(0x0),
74     fPHOSReconstructor(0x0),
75     fNSignal(0),
76     fNCaloClustersOld(0),
77     fInitialized(0)   
78 {
79   // Constructor
80   for(int i=0;i<5;i++)fOldPHOSCalibration[i]=0x0 ;
81
82 }
83
84 //________________________________________________________________________
85 void AliPHOSEmbedding::UserCreateOutputObjects()
86 {
87  //prepare output containers
88   
89   //Create branch for output MC particles  
90   TClonesArray *mcparticles = new TClonesArray("AliAODMCParticle",0);
91   mcparticles->SetName(AliAODMCParticle::StdBranchName());
92   AddAODBranch("TClonesArray", &mcparticles);
93
94   //Branch with clusters after embedding
95   fEmbeddedClusters = new TClonesArray("AliAODCaloCluster",0);
96   fEmbeddedClusters->SetName("EmbeddedCaloClusters");
97   AddAODBranch("TClonesArray", &fEmbeddedClusters);
98   
99
100   //Branch with cells after embedding
101   fEmbeddedCells = new AliAODCaloCells();
102   fEmbeddedCells->SetName("EmbeddedPHOScells");
103   AddAODBranch("AliAODCaloCells", &fEmbeddedCells);
104
105   
106   fDigitsArr = new TClonesArray("AliPHOSDigit",3*56*64) ;
107
108   AliAODHandler* handler = (AliAODHandler*) ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
109   if (handler) {
110       fTreeOut = handler->GetTree();
111   }
112   else {
113        AliWarning("No output AOD Event Handler connected.") ;
114   }
115
116 }
117
118 //________________________________________________________________________
119 void AliPHOSEmbedding::UserExec(Option_t *) {
120   // Main loop, called for each event
121   // Perform embedding
122   
123   AliESDEvent *event = dynamic_cast<AliESDEvent*>(InputEvent());
124   if (!event) {
125     Printf("ERROR: Could not retrieve event");
126     PostData(0, fTreeOut);
127     return;
128   }
129
130   AliCentrality *centrality = event->GetCentrality(); 
131   Float_t acentrality=centrality->GetCentralityPercentileUnchecked("V0M");
132   if( acentrality <= 0. || acentrality>80.){
133     return;
134   }
135   
136   //Read next AOD event
137   //If necesary method checks if AOD event is good to embed
138   //e.g. there are PHOS clusters etc.
139   AliAODEvent * signal = GetNextSignalEvent() ;
140   if (!signal) {
141     Printf("ERROR: Could not retrieve signal");
142     PostData(0, fTreeOut);
143     return;
144   }
145   
146   Init() ;
147
148   //Remember PHOS digits before any modification
149   CopyRecalibrateDigits() ;
150
151   //First re-reconstruct existing digits to ensure identical reconstruction before and 
152   //after embeding
153   AliAODEvent * nullSignal=0x0 ;
154   MakeEmbedding(event,nullSignal) ;
155
156   //Convert ESD + MC information from signal to AOD
157   ConvertESDtoAOD(event);
158   
159   //Now perform real embedding  
160   //Embed signal clusters into event
161   MakeEmbedding(event,signal) ;
162
163   //Convert ESD to AOD
164   ConvertEmbeddedClusters(event) ;
165  
166   //Fill MC information
167   ConvertMCParticles(signal) ;
168
169   delete signal ;
170
171   PostData(0, fTreeOut);
172 }
173 //______________________________________________________________________________
174 void AliPHOSEmbedding::CopyRecalibrateDigits(){
175   //Recalibrate digits if there is better calibration ("new") 
176   //exists than one used in reconstruction ("ESD")
177   
178   AliESDEvent *event = dynamic_cast<AliESDEvent*>(InputEvent());
179   if(!event){
180     AliError("No ESD event!");
181     return ;
182   }
183   if(fCellsPHOS)
184     delete fCellsPHOS ;
185   fCellsPHOS = new AliESDCaloCells() ;
186   fCellsPHOS->CreateContainer(event->GetPHOSCells()->GetNumberOfCells());
187
188   //Apply recalibration if necessary to account for
189   //difference between calibration used in reconstruction and 
190   //final calibration
191 //  if(fOldPHOSCalibration[0]){ //there is a difference in calibrations
192     for (Short_t icell = 0; icell < event->GetPHOSCells()->GetNumberOfCells(); icell++) {
193       Short_t id=0;
194       Double_t time=0., amp=0. ;
195       Int_t mclabel; 
196       Double_t efrac =0. ;
197 //      if (fCellsPHOS->GetCell(icell, id, amp,time,mclabel,efrac) != kTRUE)
198       if (event->GetPHOSCells()->GetCell(icell, id, amp,time,mclabel,efrac) != kTRUE)
199         break;      
200       Int_t relId[4] ;
201       AliPHOSGeometry::GetInstance()->AbsToRelNumbering(id,relId);
202       Float_t calibESD=fOldPHOSCalibration[relId[0]-1]->GetBinContent(relId[2],relId[3]) ;
203       Float_t calibNew=fPHOSReconstructor->Calibrate(1.,id) ;
204       if(calibESD>0.)
205         amp=amp*calibNew/calibESD ;
206       fCellsPHOS->SetCell(icell, id, amp, time,mclabel,efrac);     
207     }
208 //  }  
209 }
210 //______________________________________________________________________________
211 void AliPHOSEmbedding::Init(){
212   if(fInitialized)
213     return ;
214   
215   AliESDEvent *event = dynamic_cast<AliESDEvent*>(InputEvent());
216   if(!event){
217     AliError("Can not obtain InputEvent!") ;
218     return ;
219   }
220   Int_t runNum = event->GetRunNumber();
221   AliCDBManager::Instance()->SetRun(runNum);
222
223   fPHOSReconstructor = new AliPHOSReconstructor() ;
224   AliCDBPath path("PHOS","Calib","RecoParam");
225   AliCDBEntry *entry=AliCDBManager::Instance()->Get(path.GetPath());
226   if(!entry){
227     AliError(Form("Can not get OCDB entry %s",path.GetPath().Data())) ;
228     return ;
229   }
230   
231   TObjArray* recoParamArray = (TObjArray*)entry->GetObject();
232   AliPHOSRecoParam* recoParam = (AliPHOSRecoParam*)recoParamArray->At(2); 
233   fPHOSReconstructor->SetRecoParam(recoParam) ;
234   
235   InitMF() ;  
236   InitGeometry() ;
237
238   fInitialized=kTRUE ;
239 }
240
241 //________________________________________________________________________
242 AliAODEvent* AliPHOSEmbedding::GetNextSignalEvent(){
243   //Read signal AOD event from the chain
244   
245   if(fAODChain==0){
246     AliError("No chain to read signal events") ;
247     return 0 ;
248   }
249   
250   AliAODEvent* aod = new AliAODEvent;
251   aod->ReadFromTree(fAODChain);
252   
253   if(fNSignal>=fAODChain->GetEntries()){
254     delete aod ;
255     return 0 ;
256   }
257   fAODChain->GetEvent(fNSignal) ;
258   fNSignal++ ;
259
260   //check if event is read and there is at least 2 PHOS clusters
261   //Otherwise fill histogram and read next event
262   //TODO
263
264   return aod ;
265   
266 }
267
268 //________________________________________________________________________
269 void AliPHOSEmbedding::MakeEmbedding(AliESDEvent *event,AliAODEvent * signal){
270   //Perform embedding of the signal to the event
271
272   gROOT->cd(); //make sure that the digits and RecPoints Trees are memory resident
273   if(fDigitsTree){
274     delete fDigitsTree ;
275   }
276   fDigitsTree = new TTree("digitstree","digitstree");
277   fDigitsTree->Branch("PHOS","TClonesArray", &fDigitsArr, 32000);
278   
279   if(fClustersTree){
280     delete fClustersTree ;
281   }
282   fClustersTree = new TTree("clustertree","clustertree");
283
284   //Remember number of Clusters before we added new ones...
285   fNCaloClustersOld=event->GetNumberOfCaloClusters() ;
286
287
288 /*
289   printf("---------------Before---------\n") ;
290   Int_t n=event->GetNumberOfCaloClusters() ;
291   for (Int_t iClust=0; iClust<event->GetNumberOfCaloClusters(); ++iClust) {
292     AliESDCaloCluster * cluster = event->GetCaloCluster(iClust);
293      if(!cluster->IsPHOS())
294       continue;
295    Float_t pos[3] = { 0.};
296     cluster->GetPosition(pos);
297     printf("Cluster(%d): E=%f, x=%f, z=%f, CPV=%f \n",iClust,cluster->E(),pos[0],pos[2],cluster->GetEmcCpvDistance()) ;
298     UShort_t * index    = cluster->GetCellsAbsId() ;
299     for(Int_t ic=0; ic < cluster->GetNCells(); ic++ )
300     printf("Dig(%d)=%d, ",ic,index[ic]) ;
301   printf("\n") ;
302         
303         
304   }
305   printf("---------------END Before---------\n") ;
306 */
307
308   
309   //create digits
310   MakeDigits(signal);     
311
312 /*  
313 for(Int_t i=0;i<fDigitsArr->GetEntriesFast();i++){  
314    AliPHOSDigit * d = (AliPHOSDigit *) fDigitsArr->At(i) ;
315    printf("  Digit(%d) = %d, E=%f \n",i,d->GetId(),d->GetEnergy()) ;
316 }
317 */
318   
319   //clusterize and make tracking
320   fPHOSReconstructor->Reconstruct(fDigitsTree,fClustersTree) ;
321
322   //Note that the cirrent ESDEvent will be modified!
323   fPHOSReconstructor->FillESD(fDigitsTree, fClustersTree, event) ;
324
325
326 /*
327    printf("---------------After---------\n") ;
328   for (Int_t iClust=n; iClust<event->GetNumberOfCaloClusters(); ++iClust) {
329     AliESDCaloCluster * cluster = event->GetCaloCluster(iClust);
330     if(!cluster->IsPHOS())
331       continue;
332     Float_t pos[3] = { 0.};
333     cluster->GetPosition(pos);
334     printf("Cluster(%d): E=%f, x=%f, z=%f, CPV=%f, Label=%d \n",iClust-n,cluster->E(),pos[0],pos[2],cluster->GetEmcCpvDistance(),cluster->GetLabel()) ;
335     UShort_t * index    = cluster->GetCellsAbsId() ;
336     for(Int_t ic=0; ic < cluster->GetNCells(); ic++ )
337       printf("Dig(%d)=%d, ",ic,index[ic]) ;
338   printf("\n") ;
339   }
340   printf("---------------END After---------\n") ;
341 */
342
343
344 }
345 //______________________________________________________________________________
346 void AliPHOSEmbedding::ConvertESDtoAOD(AliESDEvent* esd) 
347 {
348   // ESD Filter analysis task executed for each event
349   
350   if(!esd)return;
351     
352   AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kTRUE);
353   AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillExtension(kTRUE);
354       
355   ConvertHeader(*esd);
356   
357   // Fetch Stack for debuggging if available 
358   
359   // Update the header
360   Int_t nTracks    = esd->GetNumberOfTracks();
361   Int_t nV0s      = esd->GetNumberOfV0s();
362   Int_t nCascades = esd->GetNumberOfCascades();
363   Int_t nKinks    = esd->GetNumberOfKinks();
364   Int_t nVertices = nV0s + nCascades /*V0 wihtin cascade already counted*/+ nKinks + 1 /* = prim. vtx*/;
365   Int_t nPileSPDVertices=1+esd->GetNumberOfPileupVerticesSPD(); // also SPD main vertex
366   Int_t nPileTrkVertices=esd->GetNumberOfPileupVerticesTracks();
367   nVertices+=nPileSPDVertices;
368   nVertices+=nPileTrkVertices;
369   Int_t nJets     = 0;
370   Int_t nCaloClus = esd->GetNumberOfCaloClusters();
371   Int_t nFmdClus  = 0;
372   Int_t nPmdClus  = esd->GetNumberOfPmdTracks();
373     
374        
375   AODEvent()->ResetStd(nTracks, nVertices, nV0s, nCascades, nJets, nCaloClus, nFmdClus, nPmdClus);
376       
377   ConvertPrimaryVertices(*esd);
378   
379   ConvertCaloClusters(*esd);
380   
381   ConvertEMCALCells(*esd);
382   
383   ConvertPHOSCells(*esd);
384   
385 }
386 //______________________________________________________________________________
387 void AliPHOSEmbedding::ConvertPHOSCells(const AliESDEvent& esd)
388 {
389   // fill PHOS cell info
390   if (esd.GetPHOSCells()) { // protection against missing ESD information
391     AliESDCaloCells &esdPHcells = *(esd.GetPHOSCells());
392     Int_t nPHcell = esdPHcells.GetNumberOfCells() ;
393     
394     AliAODCaloCells &aodPHcells = *(AODEvent()->GetPHOSCells());
395     aodPHcells.CreateContainer(nPHcell);
396     aodPHcells.SetType(AliAODCaloCells::kPHOSCell);
397     for (Int_t iCell = 0; iCell < nPHcell; iCell++) {      
398       aodPHcells.SetCell(iCell,esdPHcells.GetCellNumber(iCell),esdPHcells.GetAmplitude(iCell),esdPHcells.GetTime(iCell),-1,0.);
399     }
400     aodPHcells.Sort();
401   }
402 }
403
404 //______________________________________________________________________________
405 void AliPHOSEmbedding::ConvertEMCALCells(const AliESDEvent& esd)
406 {
407   // fill EMCAL cell info
408   if (esd.GetEMCALCells()) { // protection against missing ESD information
409     AliESDCaloCells &esdEMcells = *(esd.GetEMCALCells());
410     Int_t nEMcell = esdEMcells.GetNumberOfCells() ;
411     
412     AliAODCaloCells &aodEMcells = *(AODEvent()->GetEMCALCells());
413     aodEMcells.CreateContainer(nEMcell);
414     aodEMcells.SetType(AliAODCaloCells::kEMCALCell);
415     for (Int_t iCell = 0; iCell < nEMcell; iCell++) {      
416       aodEMcells.SetCell(iCell,esdEMcells.GetCellNumber(iCell),esdEMcells.GetAmplitude(iCell),esdEMcells.GetTime(iCell));
417     }
418     aodEMcells.Sort();
419   }
420 }
421 //______________________________________________________________________________
422 void AliPHOSEmbedding::ConvertCaloClusters(const AliESDEvent& esd)
423 {
424   
425   // Access to the AOD container of clusters
426   TClonesArray &caloClusters = *(AODEvent()->GetCaloClusters());
427   Int_t jClusters(0);
428   
429   for (Int_t iClust=fNCaloClustersOld; iClust<esd.GetNumberOfCaloClusters(); ++iClust) {
430     
431     AliESDCaloCluster * cluster = esd.GetCaloCluster(iClust);
432     
433     Int_t  id        = cluster->GetID();
434     Int_t  nLabel    = cluster->GetNLabels();
435     Int_t *labels    = cluster->GetLabels();
436     
437     Float_t energy = cluster->E();
438     Float_t posF[3] = { 0.};
439     cluster->GetPosition(posF);
440     
441     AliAODCaloCluster *caloCluster = new(caloClusters[jClusters++]) AliAODCaloCluster(id,
442                                                                                       nLabel,
443                                                                                       labels,
444                                                                                       energy,
445                                                                                       posF,
446                                                                                       NULL,
447                                                                                       cluster->GetType(),0);
448     
449     Double_t dx=cluster->GetTrackDx() ;
450     Double_t dz=cluster->GetTrackDz() ;
451     Float_t cpv=99. ; //No track matched by default
452     TArrayI * itracks = cluster->GetTracksMatched() ;  
453     if(itracks->GetSize()>0){
454       Int_t iTr = itracks->At(0);
455       if(iTr>=0 && iTr<esd.GetNumberOfTracks()){
456         AliESDtrack *track = esd.GetTrack(iTr) ;
457         Double_t pt = track->Pt() ;
458         Short_t charge = track->Charge() ;
459         cpv=TestCPV(dx, dz, pt,charge) ;
460       }
461     }
462     
463     caloCluster->SetCaloCluster(cluster->GetDistanceToBadChannel(),
464                                 cluster->GetDispersion(),
465                                 cluster->GetM20(), cluster->GetM02(),
466                                 cpv,  
467                                 cluster->GetNExMax(),cluster->GetTOF()) ;
468     
469     caloCluster->SetPIDFromESD(cluster->GetPID());
470     caloCluster->SetNCells(cluster->GetNCells());
471     caloCluster->SetCellsAbsId(cluster->GetCellsAbsId());
472     caloCluster->SetCellsAmplitudeFraction(cluster->GetCellsAmplitudeFraction());
473
474     /*    
475     TArrayI* matchedT =         cluster->GetTracksMatched();
476     if (fNumberOfTracks>0 && matchedT && cluster->GetTrackMatchedIndex() >= 0) {        
477       for (Int_t im = 0; im < matchedT->GetSize(); im++) {
478         Int_t iESDtrack = matchedT->At(im);;
479         if (fAODTrackRefs->At(iESDtrack) != 0) {
480           caloCluster->AddTrackMatched((AliAODTrack*)fAODTrackRefs->At(iESDtrack));
481         }
482       }
483     }
484     */
485     
486   } 
487   caloClusters.Expand(jClusters); // resize TObjArray to 'remove' slots for pseudo clusters      
488 }
489 //______________________________________________________________________________
490 void AliPHOSEmbedding::ConvertEmbeddedClusters(const AliESDEvent* esd)
491 {
492    //Copy PHOS clusters and cells after embedding
493   
494   // Access to the AOD container of clusters
495   Int_t jClusters(0);
496   fEmbeddedClusters->Clear();
497   fEmbeddedClusters->Expand(esd->GetNumberOfCaloClusters()-fNCaloClustersOld) ;
498   for (Int_t iClust=fNCaloClustersOld; iClust<esd->GetNumberOfCaloClusters(); ++iClust) {
499     
500     AliESDCaloCluster * cluster = esd->GetCaloCluster(iClust);
501     
502     Int_t  id        = cluster->GetID();
503     Int_t  nLabel    = cluster->GetNLabels();
504     Int_t *labels    = cluster->GetLabels();
505     
506     Float_t energy = cluster->E();
507     Float_t posF[3] = { 0.};
508     cluster->GetPosition(posF);
509     
510     AliAODCaloCluster *caloCluster = new((*fEmbeddedClusters)[jClusters++]) AliAODCaloCluster(id,
511                                                                                       nLabel,
512                                                                                       labels,
513                                                                                       energy,
514                                                                                       posF,
515                                                                                       NULL,
516                                                                                       cluster->GetType(),0);
517
518     
519     Double_t dx=cluster->GetTrackDx() ;
520     Double_t dz=cluster->GetTrackDz() ;
521     Float_t cpv=99. ; //No track matched by default
522     TArrayI * itracks = cluster->GetTracksMatched() ;  
523     if(itracks->GetSize()>0){
524       Int_t iTr = itracks->At(0);
525       if(iTr>=0 && iTr<esd->GetNumberOfTracks()){
526         AliESDtrack *track = esd->GetTrack(iTr) ;
527         Double_t pt = track->Pt() ;
528         Short_t charge = track->Charge() ;
529         cpv=TestCPV(dx, dz, pt,charge) ;
530       }
531     }
532
533     caloCluster->SetCaloCluster(cluster->GetDistanceToBadChannel(),
534                                 cluster->GetDispersion(),
535                                 cluster->GetM20(), cluster->GetM02(),
536                                 cpv,  
537                                 cluster->GetNExMax(),cluster->GetTOF()) ;
538     
539     caloCluster->SetPIDFromESD(cluster->GetPID());
540     caloCluster->SetNCells(cluster->GetNCells());    
541     caloCluster->SetCellsAbsId(cluster->GetCellsAbsId());
542     caloCluster->SetCellsAmplitudeFraction(cluster->GetCellsAmplitudeFraction());            
543   } 
544   
545   for(Int_t i=0;i<fEmbeddedClusters->GetEntriesFast();i++){
546     AliAODCaloCluster *caloCluster =(AliAODCaloCluster *)fEmbeddedClusters->At(i) ;
547     caloCluster->GetID() ;
548     
549   }
550   
551   //Now cells
552   if (esd->GetPHOSCells()) { // protection against missing ESD information
553     AliESDCaloCells &esdPHcells = *(esd->GetPHOSCells());
554     Int_t nPHcell = esdPHcells.GetNumberOfCells() ;    
555     fEmbeddedCells->CreateContainer(nPHcell);
556     fEmbeddedCells->SetType(AliAODCaloCells::kPHOSCell);
557     for (Int_t iCell = 0; iCell < nPHcell; iCell++) {      
558       fEmbeddedCells->SetCell(iCell,esdPHcells.GetCellNumber(iCell),esdPHcells.GetAmplitude(iCell),esdPHcells.GetTime(iCell));
559     }
560     fEmbeddedCells->Sort();
561   }
562
563   
564 }
565 //______________________________________________________________________________
566 void AliPHOSEmbedding::ConvertMCParticles(const AliAODEvent* aod)
567 {
568   //Copy MC branches to new AOD
569   
570   TClonesArray *mcArray = (TClonesArray*)aod->FindListObject(AliAODMCParticle::StdBranchName());
571   TClonesArray * aodMcParticles = (TClonesArray*)AODEvent()->FindListObject(AliAODMCParticle::StdBranchName());
572   for(Int_t i=0;i<mcArray->GetEntriesFast();i++){
573     
574     AliAODMCParticle* aodpart =  (AliAODMCParticle*) mcArray->At(i);
575 //    AliAODMCParticle * mcp =new ((*aodMcParticles)[i]) AliAODMCParticle(*aodpart);
576     new ((*aodMcParticles)[i]) AliAODMCParticle(*aodpart);
577     
578   }
579
580   /*  
581   // MC header
582   AliAODMCHeader * mcHeader = (AliAODMCHeader *)aod.FindListObject(AliAODMCHeader::StdBranchName());
583   AliAODMCHeader * aodMcHeader = new AliAODMCHeader(*mcHeader);
584   aodMcHeader->SetName(AliAODMCHeader::StdBranchName());
585   AddAODBranch("AliAODMCHeader",&aodMcHeader);    
586   */
587
588 }
589 //______________________________________________________________________________
590 void AliPHOSEmbedding::ConvertHeader(AliESDEvent & esd){
591   
592   AliAODHeader* header = dynamic_cast<AliAODHeader*>(AODEvent()->GetHeader());
593   if(!header) AliFatal("Not a standard AOD");
594   
595   header->SetRunNumber(esd.GetRunNumber());
596   header->SetOfflineTrigger(fInputHandler->IsEventSelected()); // propagate the decision of the physics selection
597
598   TTree* tree = fInputHandler->GetTree();
599   if (tree) {
600     TFile* file = tree->GetCurrentFile();
601     if (file) header->SetESDFileName(file->GetName());
602   }
603   
604   header->SetBunchCrossNumber(esd.GetBunchCrossNumber());
605   header->SetOrbitNumber(esd.GetOrbitNumber());
606   header->SetPeriodNumber(esd.GetPeriodNumber());
607   header->SetEventType(esd.GetEventType());
608   
609   header->SetEventNumberESDFile(esd.GetHeader()->GetEventNumberInFile());
610   if(const_cast<AliESDEvent&>(esd).GetCentrality()){
611     header->SetCentrality(const_cast<AliESDEvent&>(esd).GetCentrality());
612   }
613   else{
614     header->SetCentrality(0);
615   }
616   if(const_cast<AliESDEvent&>(esd).GetEventplane()){
617     header->SetEventplane(const_cast<AliESDEvent&>(esd).GetEventplane());
618   }
619   else{
620     header->SetEventplane(0);
621   }
622
623   // Trigger
624   header->SetFiredTriggerClasses(esd.GetFiredTriggerClasses());
625   header->SetTriggerMask(esd.GetTriggerMask()); 
626   header->SetTriggerCluster(esd.GetTriggerCluster());
627   header->SetL0TriggerInputs(esd.GetHeader()->GetL0TriggerInputs());    
628   header->SetL1TriggerInputs(esd.GetHeader()->GetL1TriggerInputs());    
629   header->SetL2TriggerInputs(esd.GetHeader()->GetL2TriggerInputs());    
630   
631   header->SetMagneticField(esd.GetMagneticField());
632   header->SetMuonMagFieldScale(esd.GetCurrentDip()/6000.);
633   header->SetZDCN1Energy(esd.GetZDCN1Energy());
634   header->SetZDCP1Energy(esd.GetZDCP1Energy());
635   header->SetZDCN2Energy(esd.GetZDCN2Energy());
636   header->SetZDCP2Energy(esd.GetZDCP2Energy());
637   header->SetZDCEMEnergy(esd.GetZDCEMEnergy(0),esd.GetZDCEMEnergy(1));
638   
639   // ITS Cluster Multiplicty
640   const AliMultiplicity *mult = esd.GetMultiplicity();
641   for (Int_t ilay = 0; ilay < 6; ilay++) header->SetITSClusters(ilay, mult->GetNumberOfITSClusters(ilay));
642   
643   // TPC only Reference Multiplicty
644   //  Int_t refMult  = fTPCaloneTrackCuts ? (Short_t)fTPCaloneTrackCuts->GetReferenceMultiplicity(&esd, kTRUE) : -1;
645   header->SetTPConlyRefMultiplicity(-1);
646   
647   //
648   Float_t diamxy[2]={static_cast<Float_t>(esd.GetDiamondX()),static_cast<Float_t>(esd.GetDiamondY())};
649   Float_t diamcov[3]; 
650   esd.GetDiamondCovXY(diamcov);
651   header->SetDiamond(diamxy,diamcov);
652   header->SetDiamondZ(esd.GetDiamondZ(),esd.GetSigma2DiamondZ());
653
654   
655   //PHOS matrixes
656   char path[255] ;
657   TGeoHMatrix * m ;
658   for(Int_t mod=0; mod<5; mod++){
659     snprintf(path,255,"/ALIC_1/PHOS_%d",mod+1) ; //In Geometry modules numbered 1,2,.,5
660     if (gGeoManager->cd(path)){
661       m = gGeoManager->GetCurrentMatrix() ;
662       header->SetPHOSMatrix(new TGeoHMatrix(*m),mod) ;
663     }
664     else{
665       header->SetPHOSMatrix(NULL,mod) ;
666     }
667   } 
668   
669   //EventPlane
670   AliEventplane *eventplane = esd.GetEventplane();
671   Double_t epQ = eventplane->GetEventplane("Q"); 
672   
673   //Standard Eventplane setter is too complicated...
674   //use durty hack
675   header->SetZDCN1Energy(epQ) ;
676   
677   AliCentrality *centrality = esd.GetCentrality(); 
678   Double_t acentrality=centrality->GetCentralityPercentileUnchecked("V0M");
679   header->SetZDCN2Energy(acentrality) ;
680
681   
682 }//______________________________________________________________________________
683 void AliPHOSEmbedding::ConvertPrimaryVertices(AliESDEvent const&esd){
684   // Access to the AOD container of vertices
685
686   // Access to the AOD container of vertices
687   Int_t fNumberOfVertices = 0;
688     
689   Double_t pos[3] = { 0. };
690   Double_t covVtx[6] = { 0. };
691
692   // Add primary vertex. The primary tracks will be defined
693   // after the loops on the composite objects (V0, cascades, kinks)
694   const AliESDVertex *vtx = esd.GetPrimaryVertex();
695   
696   vtx->GetXYZ(pos); // position
697   vtx->GetCovMatrix(covVtx); //covariance matrix
698   
699   TClonesArray * vertexes=AODEvent()->GetVertices();
700
701   AliAODVertex* primaryVertex = new((*vertexes)[fNumberOfVertices++])
702   AliAODVertex(pos, covVtx, vtx->GetChi2toNDF(), NULL, -1, AliAODVertex::kPrimary);
703   primaryVertex->SetName(vtx->GetName());
704   primaryVertex->SetTitle(vtx->GetTitle());
705   
706   TString vtitle = vtx->GetTitle();
707   if (!vtitle.Contains("VertexerTracks")) 
708     primaryVertex->SetNContributors(vtx->GetNContributors());
709   
710   if (fDebug > 0) primaryVertex->Print();  
711   
712   // Add SPD "main" vertex 
713   const AliESDVertex *vtxS = esd.GetPrimaryVertexSPD();
714   vtxS->GetXYZ(pos); // position
715   vtxS->GetCovMatrix(covVtx); //covariance matrix
716   AliAODVertex * mVSPD = new((*vertexes)[fNumberOfVertices++])
717   AliAODVertex(pos, covVtx, vtxS->GetChi2toNDF(), NULL, -1, AliAODVertex::kMainSPD);
718   mVSPD->SetName(vtxS->GetName());
719   mVSPD->SetTitle(vtxS->GetTitle());
720   mVSPD->SetNContributors(vtxS->GetNContributors()); 
721   
722   // Add SPD pileup vertices
723   for(Int_t iV=0; iV<esd.GetNumberOfPileupVerticesSPD(); ++iV)
724   {
725     const AliESDVertex *vtxP = esd.GetPileupVertexSPD(iV);
726     vtxP->GetXYZ(pos); // position
727     vtxP->GetCovMatrix(covVtx); //covariance matrix
728     AliAODVertex * pVSPD =  new((*vertexes)[fNumberOfVertices++])
729     AliAODVertex(pos, covVtx, vtxP->GetChi2toNDF(), NULL, -1, AliAODVertex::kPileupSPD);
730     pVSPD->SetName(vtxP->GetName());
731     pVSPD->SetTitle(vtxP->GetTitle());
732     pVSPD->SetNContributors(vtxP->GetNContributors()); 
733   }
734   
735   // Add TRK pileup vertices
736   for(Int_t iV=0; iV<esd.GetNumberOfPileupVerticesTracks(); ++iV)
737   {
738     const AliESDVertex *vtxP = esd.GetPileupVertexTracks(iV);
739     vtxP->GetXYZ(pos); // position
740     vtxP->GetCovMatrix(covVtx); //covariance matrix
741     AliAODVertex * pVTRK = new((*vertexes)[fNumberOfVertices++])
742     AliAODVertex(pos, covVtx, vtxP->GetChi2toNDF(), NULL, -1, AliAODVertex::kPileupTracks);
743     pVTRK->SetName(vtxP->GetName());
744     pVTRK->SetTitle(vtxP->GetTitle());
745     pVTRK->SetNContributors(vtxP->GetNContributors());
746   }
747
748   
749 }
750 //__________________________________________________________________________________
751 void AliPHOSEmbedding::MakeDigits(AliAODEvent * signal){
752   
753   //-------------------------------------------------------------------------------------
754   //Transform CaloCells into Digits which can be used for standard reconstruction
755   //Add signal digits to the event
756   //-------------------------------------------------------------------------------------
757   
758   fDigitsArr->Clear() ;
759   fDigitsTree->Branch("PHOS","TClonesArray", &fDigitsArr, 32000);
760
761   //First copy data digits
762   Int_t ndigit=0 ;
763   for (Short_t icell = 0; icell < fCellsPHOS->GetNumberOfCells(); icell++) {
764     Short_t id=0;
765     Double_t time=0., amp=0. ;
766     Int_t mclabel;
767     Double_t efrac =0. ;
768     if (fCellsPHOS->GetCell(icell, id, amp, time,mclabel, efrac) != kTRUE)
769       break;
770         
771     new((*fDigitsArr)[ndigit]) AliPHOSDigit(-1,id,float(amp),float(time),ndigit);
772     ndigit++;
773   }
774
775
776   //Add Digits from Signal
777   TClonesArray sdigits("AliPHOSDigit",0) ;
778   Int_t isdigit=0 ;
779   if(signal){
780     AliAODCaloCells* cellsS = signal->GetPHOSCells();
781     Int_t cellLabels[1000]={0} ;       //1000 should be enough for simulated
782 //    Int_t cellSecondLabels[1000]={0} ; //low-statistics event.
783     for(Int_t i=0;i<cellsS->GetNumberOfCells();i++){
784       cellLabels[i]=-1 ;
785 //      cellSecondLabels[i]=-1;
786     }
787     //------------------------------------------------------------------------------------
788     //Ancestry information
789     //Celect digits contributing to signal clusters and add primary information
790     //(it is not stored in CaloCells)
791     //------------------------------------------------------------------------------------
792     sdigits.Expand(cellsS->GetNumberOfCells());
793     for(Int_t i=0; i<signal->GetNumberOfCaloClusters(); i++) {    
794       //cluster from embedded signal
795       AliVCluster *clus = signal->GetCaloCluster(i);  
796     
797       if(!clus->IsPHOS())
798         continue;
799     
800       Int_t label = clus->GetLabel();
801 //      Int_t label2 = -1 ;
802 //      if (clus->GetNLabels()>=2) label2 = clus->GetLabelAt(1) ;
803     
804       UShort_t * index    = clus->GetCellsAbsId() ;
805       for(Int_t ic=0; ic < clus->GetNCells(); ic++ ){
806         for (Int_t icell = 0; icell < cellsS->GetNumberOfCells(); icell++){
807            Short_t cellNumber;
808            Double_t cellAmplitude=0., cellTime=0. ;
809            Int_t mclabel;
810            Double_t efrac =0. ;
811            cellsS->GetCell(icell, cellNumber, cellAmplitude, cellTime,mclabel,efrac) ;
812            if(cellNumber==index[ic]){
813               cellLabels[icell]=label;
814 //              cellSecondLabels[icell]=label2;
815               break ;
816           }
817         }
818       }
819     }
820
821     for (Int_t icell = 0; icell < cellsS->GetNumberOfCells(); icell++) {
822       Short_t cellNumber;
823       Double_t cellAmplitude=0., cellTime=0. ;
824       Int_t mclabel;
825       Double_t efrac =0. ;
826       if (cellsS->GetCell(icell, cellNumber, cellAmplitude, cellTime,mclabel,efrac) != kTRUE)
827         break;
828       //Add only digits related to the cluster, no noisy digits...
829       if(cellLabels[icell]==-1){
830         continue ;
831       }    
832       new(sdigits[isdigit]) AliPHOSDigit(cellLabels[icell],cellNumber,float(cellAmplitude),float(cellTime),isdigit);    
833       isdigit++;
834     }
835   }
836
837   
838   //Merge digits  
839   Int_t icurrent = 0 ; //index of the last used digit in underlying event
840   fDigitsArr->Expand(ndigit+isdigit) ;
841   for(Int_t i=0; i<isdigit;i++){
842     AliPHOSDigit * sdigit=static_cast<AliPHOSDigit*>(sdigits.At(i)) ;
843     Bool_t added=kFALSE ;
844     for(Int_t id=icurrent;id<ndigit;id++){
845       AliPHOSDigit * digit=static_cast<AliPHOSDigit*>(fDigitsArr->At(id)) ;
846       if(sdigit->GetId() ==  digit->GetId() ){
847         *digit += *sdigit ;  //add energies
848         icurrent=id+1 ;
849         added=kTRUE ;
850         break ; //no more digits with same ID in the list
851       }
852       if(sdigit->GetId() < digit->GetId() ){
853         icurrent=id ;
854         break ; //no more digits with same ID in the list
855       }
856     }
857     if(!added){
858       new((*fDigitsArr)[ndigit]) AliPHOSDigit(*sdigit) ;
859       ndigit++ ;
860     }
861   }
862  
863   //Change Amp back from Energy to ADC counts
864   //Note that Reconstructor uses best ("new") calibration
865   for(Int_t i=0; i<ndigit;i++){
866      AliPHOSDigit * digit=static_cast<AliPHOSDigit*>(fDigitsArr->At(i)) ;
867      Float_t calib =fPHOSReconstructor->Calibrate(1.,digit->GetId()) ;  
868      if(calib>0.)
869        digit->SetEnergy(digit->GetEnergy()/calib) ;
870   }  
871      
872   fDigitsArr->Sort() ;
873   for (Int_t i = 0 ; i < ndigit ; i++) { 
874     AliPHOSDigit *digit = static_cast<AliPHOSDigit*>( fDigitsArr->At(i) ) ; 
875     digit->SetIndexInList(i) ;     
876   }
877   fDigitsTree->Fill() ;
878
879 }
880 //__________________________________________________________________________________
881 void AliPHOSEmbedding::SetOldCalibration(TH2F **calibH){ 
882   //Copy histograms with calibration coeff
883   gROOT->cd() ;
884   char name[55] ;
885   for(Int_t i=0;i<5;i++){
886     fOldPHOSCalibration[i] = new TH2F(*((TH2F*)calibH[i])) ;
887     snprintf(name,55,"OldPHOSCalibrationMod%d",i) ;
888     fOldPHOSCalibration[i]->SetName(name) ;
889   }
890 }  
891 //____________________________________________________________________________
892 void AliPHOSEmbedding::InitMF(){
893   
894   printf("............Init MF \n") ;
895   //------------------------------------
896   // Initialization of the Mag.Fiels from GRP entry
897   // Copied from AliReconstruction::InitGRP()
898   //------------------------------------
899   AliCDBEntry* entry = AliCDBManager::Instance()->Get("GRP/GRP/Data");
900   AliGRPObject * aGRPData=0 ;
901   if (entry) {
902
903     TMap* m = dynamic_cast<TMap*>(entry->GetObject());  // old GRP entry
904
905     if (m) {
906        AliInfo("Found a TMap in GRP/GRP/Data, converting it into an AliGRPObject");
907        m->Print();
908        aGRPData = new AliGRPObject();
909        aGRPData->ReadValuesFromMap(m);
910     }
911
912     else {
913        AliInfo("Found an AliGRPObject in GRP/GRP/Data, reading it");
914        aGRPData = dynamic_cast<AliGRPObject*>(entry->GetObject());  // new GRP entry
915        entry->SetOwner(0);
916     }
917   }
918
919   if (!aGRPData) {
920      AliError("No GRP entry found in OCDB!");
921      return ;
922   }
923   //*** Dealing with the magnetic field map
924     
925   TString lhcState = aGRPData->GetLHCState();
926   if (lhcState==AliGRPObject::GetInvalidString()) {
927     AliError("GRP/GRP/Data entry:  missing value for the LHC state ! Using UNKNOWN");
928     lhcState = "UNKNOWN";
929   }
930
931   TString beamType = aGRPData->GetBeamType();
932   if (beamType==AliGRPObject::GetInvalidString()) {
933     AliError("GRP/GRP/Data entry:  missing value for the beam type ! Using UNKNOWN");
934     beamType = "UNKNOWN";
935   }
936
937   Float_t beamEnergy = aGRPData->GetBeamEnergy();
938   if (beamEnergy==AliGRPObject::GetInvalidFloat()) {
939     AliError("GRP/GRP/Data entry:  missing value for the beam energy ! Using 0");
940     beamEnergy = 0;
941   }
942
943   TString runType = aGRPData->GetRunType();
944   if (runType==AliGRPObject::GetInvalidString()) {
945     AliError("GRP/GRP/Data entry:  missing value for the run type ! Using UNKNOWN");
946     runType = "UNKNOWN";
947   }
948   
949   if ( TGeoGlobalMagField::Instance()->IsLocked() ) {
950     if (TGeoGlobalMagField::Instance()->GetField()->TestBit(AliMagF::kOverrideGRP)) {
951       AliInfo("PHOSEmbedding: MF is locked - GRP information will be ignored !");
952       AliInfo("Running with the externally locked B field !");
953     }
954     else {
955       AliInfo("Destroying existing B field instance!");
956       delete TGeoGlobalMagField::Instance();
957     }   
958   }
959   if ( !TGeoGlobalMagField::Instance()->IsLocked() ) {
960     // Construct the field map out of the information retrieved from GRP.
961     Bool_t ok = kTRUE;
962     // L3
963     Float_t l3Current = aGRPData->GetL3Current((AliGRPObject::Stats)0);
964     if (l3Current == AliGRPObject::GetInvalidFloat()) {
965       AliError("GRP/GRP/Data entry:  missing value for the L3 current !");
966       ok = kFALSE;
967     }
968
969     Char_t l3Polarity = aGRPData->GetL3Polarity();
970     if (l3Polarity == AliGRPObject::GetInvalidChar()) {
971       AliError("GRP/GRP/Data entry:  missing value for the L3 polarity !");
972       ok = kFALSE;
973     }
974
975     // Dipole
976     Float_t diCurrent = aGRPData->GetDipoleCurrent((AliGRPObject::Stats)0);
977     if (diCurrent == AliGRPObject::GetInvalidFloat()) {
978       AliError("GRP/GRP/Data entry:  missing value for the dipole current !");
979       ok = kFALSE;
980     }
981
982     Char_t diPolarity = aGRPData->GetDipolePolarity();
983     if (diPolarity == AliGRPObject::GetInvalidChar()) {
984       AliError("GRP/GRP/Data entry:  missing value for the dipole polarity !");
985       ok = kFALSE;
986     }
987
988     // read special bits for the polarity convention and map type
989     Int_t  polConvention = aGRPData->IsPolarityConventionLHC() ? AliMagF::kConvLHC : AliMagF::kConvDCS2008;
990     Bool_t uniformB = aGRPData->IsUniformBMap();
991
992     if (ok) {
993       AliMagF* fld = AliMagF::CreateFieldMap(TMath::Abs(l3Current) * (l3Polarity ? -1:1),
994                                              TMath::Abs(diCurrent) * (diPolarity ? -1:1),
995                                              polConvention,uniformB,beamEnergy, beamType.Data());
996       if (fld) {
997         TGeoGlobalMagField::Instance()->SetField( fld );
998         TGeoGlobalMagField::Instance()->Lock();
999         AliInfo("Running with the B field constructed out of GRP !");
1000       }
1001       else AliFatal("Failed to create a B field map !");
1002     }
1003     else AliFatal("B field is neither set nor constructed from GRP ! Exitig...");
1004   }
1005   
1006 }
1007 //____________________________________________________________________________
1008 void AliPHOSEmbedding::InitGeometry(){
1009
1010   // Import ideal TGeo geometry and apply misalignment
1011   if (!gGeoManager) {
1012     AliGeomManager::LoadGeometry("geometry.root");
1013     if (!gGeoManager) {
1014       AliFatal("Can not load geometry");
1015     }
1016     if(!AliGeomManager::CheckSymNamesLUT("PHOS")) {
1017       AliFatal("CheckSymNamesLUT");
1018     }
1019   }
1020     
1021
1022   TString detStr = "PHOS";
1023   TString loadAlObjsListOfDets = "PHOS";
1024         
1025   if(AliGeomManager::GetNalignable("GRP") != 0)
1026       loadAlObjsListOfDets.Prepend("GRP "); //add alignment objects for non-sensitive modules
1027 //  AliGeomManager::ApplyAlignObjsFromCDB(loadAlObjsListOfDets.Data());
1028
1029   AliCDBManager::Instance()->UnloadFromCache("*/Align/*");
1030   AliCDBManager::Instance()->UnloadFromCache("GRP/Geometry/Data");
1031  
1032 }
1033 //____________________________________________________________________________
1034 Float_t AliPHOSEmbedding::TestCPV(Double_t dx, Double_t dz, Double_t pt, Int_t charge){
1035   //Parameterization of LHC10h period
1036   //_true if neutral_
1037   
1038   Double_t meanX=0;
1039   Double_t meanZ=0.;
1040   Double_t sx=TMath::Min(5.4,2.59719e+02*TMath::Exp(-pt/1.02053e-01)+
1041               6.58365e-01*5.91917e-01*5.91917e-01/((pt-9.61306e-01)*(pt-9.61306e-01)+5.91917e-01*5.91917e-01)+1.59219);
1042   Double_t sz=TMath::Min(2.75,4.90341e+02*1.91456e-02*1.91456e-02/(pt*pt+1.91456e-02*1.91456e-02)+1.60) ;
1043   
1044   AliESDEvent *event = static_cast<AliESDEvent*>(InputEvent());
1045   Double_t mf = event->GetMagneticField(); //Positive for ++ and negative for --
1046
1047   if(mf<0.){ //field --
1048     meanZ = -0.468318 ;
1049     if(charge>0)
1050       meanX=TMath::Min(7.3, 3.89994*1.20679*1.20679/(pt*pt+1.20679*1.20679)+0.249029+2.49088e+07*TMath::Exp(-pt*3.33650e+01)) ;
1051     else
1052       meanX=-TMath::Min(7.7,3.86040*0.912499*0.912499/(pt*pt+0.912499*0.912499)+1.23114+4.48277e+05*TMath::Exp(-pt*2.57070e+01)) ;
1053   }
1054   else{ //Field ++
1055     meanZ= -0.468318;
1056     if(charge>0)
1057       meanX=-TMath::Min(8.0,3.86040*1.31357*1.31357/(pt*pt+1.31357*1.31357)+0.880579+7.56199e+06*TMath::Exp(-pt*3.08451e+01)) ;
1058     else
1059       meanX= TMath::Min(6.85, 3.89994*1.16240*1.16240/(pt*pt+1.16240*1.16240)-0.120787+2.20275e+05*TMath::Exp(-pt*2.40913e+01)) ;     
1060   }
1061
1062   Double_t rz=(dz-meanZ)/sz ;
1063   Double_t rx=(dx-meanX)/sx ;
1064   return TMath::Sqrt(rx*rx+rz*rz) ;
1065 }
1066