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