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