Event embedding tasks for PHOS are added (D.Peressounko)
[u/mrichter/AliRoot.git] / PWG4 / UserTasks / PHOS_embedding / AliAnalysisTaskPi0DiffEfficiency.cxx
1 #include "TChain.h"
2 #include "TTree.h"
3 #include "TObjArray.h"
4 #include "TF1.h"
5 #include "TFile.h"
6 #include "TH1F.h"
7 #include "TH2F.h"
8 #include "TH2I.h"
9 #include "TH3F.h"
10 #include "TParticle.h"
11 #include "TCanvas.h"
12 #include "TStyle.h"
13 #include "TRandom.h"
14
15 #include "AliAODMCParticle.h"
16 #include "AliAnalysisManager.h"
17 #include "AliMCEventHandler.h"
18 #include "AliMCEvent.h"
19 #include "AliStack.h"
20 #include "AliAnalysisTaskSE.h"
21 #include "AliAnalysisTaskPi0DiffEfficiency.h"
22 #include "AliCaloPhoton.h"
23 #include "AliPHOSGeometry.h"
24 #include "AliPHOSEsdCluster.h"
25 #include "AliPHOSCalibData.h"
26 #include "AliAODEvent.h"
27 #include "AliAODCaloCluster.h"
28 #include "AliAODVertex.h"
29 #include "AliESDtrackCuts.h"
30 #include "AliLog.h"
31 #include "AliPID.h"
32 #include "AliCDBManager.h"
33 #include "AliCentrality.h" 
34
35 // Analysis task to fill histograms with PHOS ESD clusters and cells
36 // Authors: Yuri Kharlov
37 // Date   : 28.05.2009
38
39 ClassImp(AliAnalysisTaskPi0DiffEfficiency)
40
41 //________________________________________________________________________
42 AliAnalysisTaskPi0DiffEfficiency::AliAnalysisTaskPi0DiffEfficiency(const char *name) 
43 : AliAnalysisTaskSE(name),
44   fStack(0),
45   fOutputContainer(0),
46   fPHOSEvent1(0),
47   fPHOSEvent2(0),
48   fPHOSCalibData(0),
49   fNonLinCorr(0),
50   fRPfull(0),
51   fRPA(0),
52   fRPC(0),
53   fRPFar(0),
54   fRPAFar(0),
55   fRPCFar(0),
56   fCentrality(0),
57   fCenBin(0),
58   fPHOSGeo(0),
59   fEventCounter(0)
60 {
61   // Constructor
62   for(Int_t i=0;i<1;i++){
63     for(Int_t j=0;j<10;j++)
64       for(Int_t k=0;k<11;k++)
65         fPHOSEvents[i][j][k]=0 ;
66   }
67   
68   // Output slots #0 write into a TH1 container
69   DefineOutput(1,TList::Class());
70
71   // Set bad channel map
72   char key[55] ;
73   for(Int_t i=0; i<6; i++){
74     sprintf(key,"PHOS_BadMap_mod%d",i) ;
75     fPHOSBadMap[i]=new TH2I(key,"Bad Modules map",64,0.,64.,56,0.,56.) ;
76   }
77   // Initialize the PHOS geometry
78   fPHOSGeo = AliPHOSGeometry::GetInstance("IHEP") ;
79
80
81 }
82
83 //________________________________________________________________________
84 void AliAnalysisTaskPi0DiffEfficiency::UserCreateOutputObjects()
85 {
86   // Create histograms
87   // Called once
88
89   // ESD histograms
90   if(fOutputContainer != NULL){
91     delete fOutputContainer;
92   }
93   fOutputContainer = new TList();
94   fOutputContainer->SetOwner(kTRUE);
95
96   //Event selection
97   fOutputContainer->Add(new TH1F("hSelEvents","Event celection", 10,0.,10.)) ;
98
99   //vertex distribution
100   fOutputContainer->Add(new TH1F("hZvertex","Z vertex position", 50,-25.,25.)) ;
101
102   //Centrality
103   fOutputContainer->Add(new TH1F("hCentrality","Event centrality", 100,0.,100.)) ;
104
105   //QA histograms                       
106   fOutputContainer->Add(new TH2F("hCluM1","Cell (X,Z), M1" ,64,0.5,64.5, 56,0.5,56.5));
107   fOutputContainer->Add(new TH2F("hCluM2","Cell (X,Z), M2" ,64,0.5,64.5, 56,0.5,56.5));
108   fOutputContainer->Add(new TH2F("hCluM3","Cell (X,Z), M3" ,64,0.5,64.5, 56,0.5,56.5));
109
110   Int_t nM       = 500;
111   Double_t mMin  = 0.0;
112   Double_t mMax  = 1.0;
113   Int_t nPt      = 200;
114   Double_t ptMin = 0;
115   Double_t ptMax = 20;
116
117   char key[55] ;
118   for(Int_t cent=0; cent<6; cent++){
119     //Single photon
120     snprintf(key,55,"hPhotAll_cen%d",cent) ;
121     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
122     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
123     snprintf(key,55,"hPhotCPV_cen%d",cent) ;
124     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
125     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
126     snprintf(key,55,"hPhotDisp_cen%d",cent) ;
127     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
128     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
129     snprintf(key,55,"hPhotBoth_cen%d",cent) ;
130     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
131     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
132
133     snprintf(key,55,"hOldMassPtAll_cen%d",cent) ;
134     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
135     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
136     snprintf(key,55,"hOldMassPtCPV_cen%d",cent) ;
137     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
138     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
139     snprintf(key,55,"hOldMassPtDisp_cen%d",cent) ;
140     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
141     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
142     snprintf(key,55,"hOldMassPtBoth_cen%d",cent) ;
143     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
144     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
145     
146     snprintf(key,55,"hNewMassPtAll_cen%d",cent) ;
147     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
148     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
149     snprintf(key,55,"hNewMassPtCPV_cen%d",cent) ;
150     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
151     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
152     snprintf(key,55,"hNewMassPtDisp_cen%d",cent) ;
153     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
154     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
155     snprintf(key,55,"hNewMassPtBoth_cen%d",cent) ;
156     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
157     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
158     
159     snprintf(key,55,"hMassPtAll_cen%d",cent) ;
160     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
161     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
162     snprintf(key,55,"hMassPtCPV_cen%d",cent) ;
163     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
164     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
165     snprintf(key,55,"hMassPtDisp_cen%d",cent) ;
166     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
167     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
168     snprintf(key,55,"hMassPtBoth_cen%d",cent) ;
169     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
170     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
171     
172     //Mixed
173     snprintf(key,55,"hMiMassPtAll_cen%d",cent) ;
174     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
175     snprintf(key,55,"hMiMassPtCPV_cen%d",cent) ;
176     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
177     snprintf(key,55,"hMiMassPtDisp_cen%d",cent) ;
178     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
179     snprintf(key,55,"hMiMassPtBoth_cen%d",cent) ;
180     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
181
182     snprintf(key,55,"hMCMassPtAll_cen%d",cent) ;
183     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
184     snprintf(key,55,"hMCMassPtCPV_cen%d",cent) ;
185     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
186     snprintf(key,55,"hMCMassPtDisp_cen%d",cent) ;
187     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
188     snprintf(key,55,"hMCMassPtBoth_cen%d",cent) ;
189     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
190
191     //Single photon
192     snprintf(key,55,"hMCPhotAll_cen%d",cent) ;
193     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
194     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
195     snprintf(key,55,"hMCPhotCPV_cen%d",cent) ;
196     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
197     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
198     snprintf(key,55,"hMCPhotDisp_cen%d",cent) ;
199     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
200     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
201     snprintf(key,55,"hMCPhotBoth_cen%d",cent) ;
202     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
203     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
204
205     
206     
207   }
208
209
210   //MC
211   for(Int_t cent=0; cent<6; cent++){
212     snprintf(key,55,"hMC_rap_gamma_cen%d",cent) ;
213     fOutputContainer->Add(new TH1F(key,"Rapidity pi0",200,-1.,1.)) ;
214     snprintf(key,55,"hMC_rap_pi0_cen%d",cent) ;
215     fOutputContainer->Add(new TH1F(key,"Rapidity pi0",200,-1.,1.)) ;
216     snprintf(key,55,"hMC_rap_eta_cen%d",cent) ;
217     fOutputContainer->Add(new TH1F("hMC_rap_eta","Rapidity eta",200,-1.,1.)) ;
218     snprintf(key,55,"hMC_phi_gamma_cen%d",cent) ;
219     fOutputContainer->Add(new TH1F(key,"Phi pi0",200,0.,TMath::TwoPi())) ;
220     snprintf(key,55,"hMC_phi_pi0_cen%d",cent) ;
221     fOutputContainer->Add(new TH1F(key,"Phi pi0",200,0.,TMath::TwoPi())) ;
222     snprintf(key,55,"hMC_phi_eta_cen%d",cent) ;
223     fOutputContainer->Add(new TH1F(key,"Phi eta",200,0.,TMath::TwoPi())) ;
224     snprintf(key,55,"hMC_all_gamma_cen%d",cent) ;
225     fOutputContainer->Add(new TH1F(key,"Rapidity photon",250,0.,25.)) ;
226     snprintf(key,55,"hMC_all_pi0_cen%d",cent) ;
227     fOutputContainer->Add(new TH1F(key,"Rapidity pi0",250,0.,25.)) ;
228     snprintf(key,55,"hMC_all_eta_cen%d",cent) ;
229     fOutputContainer->Add(new TH1F(key,"Rapidity eta",250,0.,25.)) ;
230     snprintf(key,55,"hMC_unitEta_gamma_cen%d",cent) ;
231     fOutputContainer->Add(new TH1F(key,"Pt photon",250,0.,25.)) ;
232     snprintf(key,55,"hMC_unitEta_pi0_cen%d",cent) ;
233     fOutputContainer->Add(new TH1F(key,"Rapidity eta",250,0.,25.)) ;
234     snprintf(key,55,"hMC_unitEta_eta_cen%d",cent) ;
235     fOutputContainer->Add(new TH1F(key,"Rapidity eta",250,0.,25.)) ;
236   }
237   
238   PostData(1, fOutputContainer);
239
240 }
241
242 //________________________________________________________________________
243 void AliAnalysisTaskPi0DiffEfficiency::UserExec(Option_t *) 
244 {
245   // Main loop, called for each event
246   // Analyze ESD/AOD
247
248   FillHistogram("hSelEvents",0.5) ;  
249   
250   AliAODEvent *event = dynamic_cast<AliAODEvent*>(InputEvent());
251   if (!event) {
252     Printf("ERROR: Could not retrieve event");
253     PostData(1, fOutputContainer);
254     return;
255   }
256   
257   FillHistogram("hSelEvents",1.5) ;
258   AliAODHeader *header = event->GetHeader() ;
259   
260   // Checks if we have a primary vertex
261   // Get primary vertices form ESD
262   const AliAODVertex *esdVertex5 = event->GetPrimaryVertex();
263
264  // don't rely on ESD vertex, assume (0,0,0)
265   Double_t vtx0[3] ={0.,0.,0.};
266   Double_t vtx5[3] ={0.,0.,0.};
267   
268   vtx5[0] = esdVertex5->GetX();
269   vtx5[1] = esdVertex5->GetY();
270   vtx5[2] = esdVertex5->GetZ();
271   
272   
273   FillHistogram("hZvertex",esdVertex5->GetZ());
274   if (TMath::Abs(esdVertex5->GetZ()) > 10. ){
275     PostData(1, fOutputContainer);
276     return;
277   }
278   FillHistogram("hSelEvents",2.5) ;
279
280   //Vtx class z-bin
281   //  Int_t zvtx = (Int_t)((vtx5[2]+10.)/2.) ;
282   //  if(zvtx<0)zvtx=0 ;
283   //  if(zvtx>9)zvtx=9 ;
284   Int_t zvtx=0 ;
285
286 //  fCentrality=header->GetCentralityP()->GetCentralityPercentile("V0M"); // returns the centrality percentile, 
287 //                                                          //a float from 0 to 100 (or to the trigger efficiency)
288    fCentrality=header->GetZDCN2Energy() ;
289
290   if( fCentrality < 0. || fCentrality>80.){
291     PostData(1, fOutputContainer);
292     return;
293   }
294   FillHistogram("hSelEvents",3.5) ;
295   Float_t bins[7]={0.,5.,10.,20.,40.,60.,80.} ;
296   fCenBin=0 ;
297   while(fCenBin<6 && fCentrality > bins[fCenBin+1])
298     fCenBin++ ; 
299
300  
301   //reaction plain
302   fRPfull= header->GetZDCN1Energy() ;
303   if(fRPfull==999){ //reaction plain was not defined
304     PostData(1, fOutputContainer);
305     return;
306   } 
307
308   FillHistogram("hSelEvents",4.5) ;
309   //All event selections done
310   FillHistogram("hCentrality",fCentrality) ;
311   //Reaction plain is defined in the range (-pi/2;pi/2)
312   //We have 10 bins
313   Int_t irp=Int_t(10.*(fRPfull+TMath::PiOver2())/TMath::Pi());
314   if(irp>9)irp=9 ;
315
316   if(!fPHOSEvents[zvtx][fCenBin][irp]) 
317     fPHOSEvents[zvtx][fCenBin][irp]=new TList() ;
318   TList * prevPHOS = fPHOSEvents[zvtx][fCenBin][irp] ;
319
320   // Get PHOS rotation matrices from ESD and set them to the PHOS geometry
321   if(fEventCounter == 0) {
322     for(Int_t mod=0; mod<5; mod++) {
323       const TGeoHMatrix* m =header->GetPHOSMatrix(mod) ;
324       fPHOSGeo->SetMisalMatrix(m,mod) ;
325       Printf("PHOS geo matrix for module # %d is set: %p\n", mod,m);
326     }
327     fEventCounter++ ;
328   }
329
330   ProcessMC() ;
331
332   if(fPHOSEvent1){
333     fPHOSEvent1->Clear() ;
334     fPHOSEvent2->Clear() ;
335   }
336   else{
337     fPHOSEvent1 = new TClonesArray("AliCaloPhoton",200) ;
338     fPHOSEvent2 = new TClonesArray("AliCaloPhoton",200) ;
339   }
340
341   TClonesArray * clustersEmb = (TClonesArray*)event->FindListObject("EmbeddedCaloClusters") ;
342   TClonesArray * clustersOld = event->GetCaloClusters() ;
343   TVector3 vertex(vtx0);
344   char key[55] ;
345   //Before Embedding
346   Int_t multClustOld = clustersOld->GetEntriesFast();
347   Int_t multClustEmb = clustersEmb->GetEntriesFast();
348   Int_t inPHOSold=0 ;
349   Int_t inPHOSemb=0 ;
350   for (Int_t i=0; i<multClustOld; i++) {
351     AliAODCaloCluster *clu = (AliAODCaloCluster*)clustersOld->At(i);
352     if ( !clu->IsPHOS() || clu->E()<0.3) continue;
353
354     Bool_t survive=kFALSE ;
355     for(Int_t ii=0;(ii<multClustEmb)&&(!survive);ii++){
356        AliAODCaloCluster *clu2 = (AliAODCaloCluster*)clustersEmb->At(ii);
357        survive=IsSameCluster(clu,clu2);
358     }
359
360
361     Float_t  position[3];
362     clu->GetPosition(position);
363     TVector3 global(position) ;
364     Int_t relId[4] ;
365     fPHOSGeo->GlobalPos2RelId(global,relId) ;
366     Int_t mod  = relId[0] ;
367     Int_t cellX = relId[2];
368     Int_t cellZ = relId[3] ;
369     if ( !IsGoodChannel("PHOS",mod,cellX,cellZ) ) 
370       continue ;
371     if(clu->GetNCells()<3)
372       continue ;
373
374     sprintf(key,"hCluM%d",mod) ;
375     FillHistogram(key,cellX,cellZ,1.);
376
377     TLorentzVector pv1 ;
378     clu->GetMomentum(pv1 ,vtx0);
379     
380     if(inPHOSold>=fPHOSEvent1->GetSize()){
381       fPHOSEvent1->Expand(inPHOSold+50) ;
382     }
383     AliCaloPhoton * ph = new((*fPHOSEvent1)[inPHOSold]) AliCaloPhoton(pv1.X(),pv1.Py(),pv1.Z(),pv1.E()) ;
384     ph->SetModule(mod) ;
385     ph->SetMomV2(&pv1) ;
386     ph->SetNCells(clu->GetNCells());
387     ph->SetDispBit(TestLambda(clu->GetM20(),clu->GetM02())) ;
388     ph->SetCPVBit(clu->GetEmcCpvDistance()>1.) ;
389     if(!survive) //this cluster found in list after embedding, skipping it
390      ph->SetTagged(1) ;
391
392     inPHOSold++ ;
393   }
394
395   for (Int_t i=0; i<multClustEmb; i++) {
396     AliAODCaloCluster *clu = (AliAODCaloCluster*)clustersEmb->At(i);
397     if ( !clu->IsPHOS() || clu->E()<0.3) continue;
398
399     Bool_t survive=kFALSE ;
400     for(Int_t ii=0;(ii<multClustOld)&&(!survive);ii++){
401        AliAODCaloCluster *clu2 = (AliAODCaloCluster*)clustersOld->At(ii);
402        survive=IsSameCluster(clu,clu2);
403     }
404     
405     Float_t  position[3];
406     clu->GetPosition(position);
407     TVector3 global(position) ;
408     Int_t relId[4] ;
409     fPHOSGeo->GlobalPos2RelId(global,relId) ;
410     Int_t mod  = relId[0] ;
411     Int_t cellX = relId[2];
412     Int_t cellZ = relId[3] ;
413     if ( !IsGoodChannel("PHOS",mod,cellX,cellZ) ) 
414       continue ;
415     if(clu->GetNCells()<3)
416       continue ;
417
418     sprintf(key,"hCluM%d",mod) ;
419     FillHistogram(key,cellX,cellZ,1.);
420
421     TLorentzVector pv1 ;
422     clu->GetMomentum(pv1 ,vtx0);
423     
424     if(inPHOSemb>=fPHOSEvent2->GetSize()){
425       fPHOSEvent2->Expand(inPHOSemb+50) ;
426     }
427     AliCaloPhoton * ph = new((*fPHOSEvent2)[inPHOSemb]) AliCaloPhoton(pv1.X(),pv1.Py(),pv1.Z(),pv1.E()) ;
428     ph->SetModule(mod) ;
429     ph->SetMomV2(&pv1) ;
430     ph->SetNCells(clu->GetNCells());
431     ph->SetDispBit(TestLambda(clu->GetM20(),clu->GetM02())) ;
432     ph->SetCPVBit(clu->GetEmcCpvDistance()>1.) ;
433     if(!survive) //this cluster found in list after embedding, skipping it
434      ph->SetTagged(1) ;
435
436     inPHOSemb++ ;
437   }
438   
439   //Single photon
440   for (Int_t i1=0; i1<inPHOSold; i1++) {
441     AliCaloPhoton * ph1=(AliCaloPhoton*)fPHOSEvent1->At(i1) ;
442     if(!ph1->IsTagged())
443       continue ;
444     snprintf(key,55,"hPhotAll_cen%d",fCenBin) ;
445     FillHistogram(key,ph1->Pt(),-1.) ;
446     if(ph1->IsCPVOK() ){
447       snprintf(key,55,"hPhotCPV_cen%d",fCenBin) ;
448       FillHistogram(key,ph1->Pt(),-1.) ;
449     }
450     if(ph1->IsDispOK()){
451       snprintf(key,55,"hPhotDisp_cen%d",fCenBin) ;
452       FillHistogram(key,ph1->Pt(),-1.) ;
453       if(ph1->IsCPVOK()){
454         snprintf(key,55,"hPhotBoth_cen%d",fCenBin) ;
455         FillHistogram(key,ph1->Pt(),-1.) ;
456       }
457     } // end of loop i2
458   } // end of loop i1 
459
460   for (Int_t i1=0; i1<inPHOSemb; i1++) {
461     AliCaloPhoton * ph1=(AliCaloPhoton*)fPHOSEvent2->At(i1) ;
462     if(!ph1->IsTagged())
463       continue ;
464     snprintf(key,55,"hPhotAll_cen%d",fCenBin) ;
465     FillHistogram(key,ph1->Pt(),1.) ;
466     if(ph1->IsCPVOK() ){
467       snprintf(key,55,"hPhotCPV_cen%d",fCenBin) ;
468       FillHistogram(key,ph1->Pt(),1.) ;
469     }
470     if(ph1->IsDispOK()){
471       snprintf(key,55,"hPhotDisp_cen%d",fCenBin) ;
472       FillHistogram(key,ph1->Pt(),1.) ;
473       if(ph1->IsCPVOK()){
474         snprintf(key,55,"hPhotBoth_cen%d",fCenBin) ;
475         FillHistogram(key,ph1->Pt(),1.) ;
476       }
477     } // end of loop i2
478   } // end of loop i1 
479
480
481
482   // Fill Real disribution:
483   // Disappeared clusters enter with negative contribution
484   // In addition fill control histogam with Real before embedding
485   for (Int_t i1=0; i1<inPHOSold-1; i1++) {
486     AliCaloPhoton * ph1=(AliCaloPhoton*)fPHOSEvent1->At(i1) ;
487     for (Int_t i2=i1+1; i2<inPHOSold; i2++) {
488       AliCaloPhoton * ph2=(AliCaloPhoton*)fPHOSEvent1->At(i2) ;
489
490       TLorentzVector p12  = *ph1  + *ph2;
491       TLorentzVector pv12 = *(ph1->GetMomV2()) + *(ph2->GetMomV2());      
492       //Fill Controll histogram: Real before embedding
493       snprintf(key,55,"hOldMassPtAll_cen%d",fCenBin) ;
494       FillHistogram(key,p12.M() ,p12.Pt(),-1.) ;
495       if(ph1->IsCPVOK() && ph2->IsCPVOK()){
496         snprintf(key,55,"hOldMassPtCPV_cen%d",fCenBin) ;
497         FillHistogram(key,p12.M() ,p12.Pt(),-1) ;
498       }
499       if(ph1->IsDispOK() && ph2->IsDispOK()){
500         snprintf(key,55,"hOldMassPtDisp_cen%d",fCenBin) ;
501         FillHistogram(key,p12.M() ,p12.Pt(),-1) ;
502
503         if(ph1->IsCPVOK() && ph2->IsCPVOK()){
504           snprintf(key,55,"hOldMassPtBoth_cen%d",fCenBin) ;
505           FillHistogram(key,p12.M() ,p12.Pt(),-1) ;
506         }
507       }
508       
509       //Now fill mail histograms with negative contributions
510       if(!(ph1->IsTagged() || ph2->IsTagged()) )
511         continue ;
512       snprintf(key,55,"hMassPtAll_cen%d",fCenBin) ;
513       FillHistogram(key,p12.M() ,p12.Pt(),-1.) ;
514       if(ph1->IsCPVOK() && ph2->IsCPVOK()){
515         snprintf(key,55,"hMassPtCPV_cen%d",fCenBin) ;
516         FillHistogram(key,p12.M() ,p12.Pt(),-1) ;
517       }
518       if(ph1->IsDispOK() && ph2->IsDispOK()){
519         snprintf(key,55,"hMassPtDisp_cen%d",fCenBin) ;
520         FillHistogram(key,p12.M() ,p12.Pt(),-1) ;
521
522         if(ph1->IsCPVOK() && ph2->IsCPVOK()){
523           snprintf(key,55,"hMassPtBoth_cen%d",fCenBin) ;
524           FillHistogram(key,p12.M() ,p12.Pt(),-1) ;
525         }
526       }
527     } // end of loop i2
528   } // end of loop i1 
529
530
531   // Further fill Real disribution
532   // now with positive contribution from new clusters
533   // ass well fill controll histogram
534   for (Int_t i1=0; i1<inPHOSemb-1; i1++) {
535     AliCaloPhoton * ph1=(AliCaloPhoton*)fPHOSEvent2->At(i1) ;
536     for (Int_t i2=i1+1; i2<inPHOSemb; i2++) {
537       AliCaloPhoton * ph2=(AliCaloPhoton*)fPHOSEvent2->At(i2) ;
538
539       TLorentzVector p12  = *ph1  + *ph2;
540       TLorentzVector pv12 = *(ph1->GetMomV2()) + *(ph2->GetMomV2());      
541
542       // Controll histogram: Real after embedding
543       snprintf(key,55,"hNewMassPtAll_cen%d",fCenBin) ;
544       FillHistogram(key,p12.M() ,p12.Pt(),1.) ;
545       if(ph1->IsCPVOK() && ph2->IsCPVOK()){
546         snprintf(key,55,"hNewMassPtCPV_cen%d",fCenBin) ;
547         FillHistogram(key,p12.M() ,p12.Pt(),1) ;
548       }
549       if(ph1->IsDispOK() && ph2->IsDispOK()){
550         snprintf(key,55,"hNewMassPtDisp_cen%d",fCenBin) ;
551         FillHistogram(key,p12.M() ,p12.Pt(),1) ;
552
553         if(ph1->IsCPVOK() && ph2->IsCPVOK()){
554           snprintf(key,55,"hNewMassPtBoth_cen%d",fCenBin) ;
555           FillHistogram(key,p12.M() ,p12.Pt(),1) ;
556         }
557       }
558      
559       //Now fill main histogamm
560       //new clusters with positive contribution
561       if(!(ph1->IsTagged() || ph2->IsTagged()) )
562         continue ;
563       snprintf(key,55,"hMassPtAll_cen%d",fCenBin) ;
564       FillHistogram(key,p12.M() ,p12.Pt(),1.) ;
565       if(ph1->IsCPVOK() && ph2->IsCPVOK()){
566         snprintf(key,55,"hMassPtCPV_cen%d",fCenBin) ;
567         FillHistogram(key,p12.M() ,p12.Pt(),1) ;
568       }
569       if(ph1->IsDispOK() && ph2->IsDispOK()){
570         snprintf(key,55,"hMassPtDisp_cen%d",fCenBin) ;
571         FillHistogram(key,p12.M() ,p12.Pt(),1) ;
572
573         if(ph1->IsCPVOK() && ph2->IsCPVOK()){
574           snprintf(key,55,"hMassPtBoth_cen%d",fCenBin) ;
575           FillHistogram(key,p12.M() ,p12.Pt(),1) ;
576         }
577       }
578     } // end of loop i2
579   } // end of loop i1 
580
581
582   //now mixed, does not really matter old or new list of clusters
583   for (Int_t i1=0; i1<inPHOSemb; i1++) {
584     AliCaloPhoton * ph1=(AliCaloPhoton*)fPHOSEvent2->At(i1) ;
585     for(Int_t ev=0; ev<prevPHOS->GetSize();ev++){
586       TClonesArray * mixPHOS = static_cast<TClonesArray*>(prevPHOS->At(ev)) ;
587       for(Int_t i2=0; i2<mixPHOS->GetEntriesFast();i2++){
588         AliCaloPhoton * ph2=(AliCaloPhoton*)mixPHOS->At(i2) ;
589         TLorentzVector p12  = *ph1  + *ph2;
590         TLorentzVector pv12 = *(ph1->GetMomV2()) + *(ph2->GetMomV2());
591         
592         snprintf(key,55,"hMiMassPtAll_cen%d",fCenBin) ;
593         FillHistogram(key,p12.M() ,p12.Pt(),1.) ;
594         if(ph1->IsCPVOK() && ph2->IsCPVOK()){
595           snprintf(key,55,"hMiMassPtCPV_cen%d",fCenBin) ;
596           FillHistogram(key,p12.M() ,p12.Pt(),1.) ;
597         }
598         if(ph1->IsDispOK() && ph2->IsDispOK()){
599           snprintf(key,55,"hMiMassPtDisp_cen%d",fCenBin) ;
600           FillHistogram(key,p12.M() ,p12.Pt(),1.) ;
601           
602           if(ph1->IsCPVOK() && ph2->IsCPVOK()){
603             snprintf(key,55,"hMiMassPtBoth_cen%d",fCenBin) ;
604             FillHistogram(key,p12.M() ,p12.Pt(),1.) ;
605           }
606         }
607       } // end of loop i2
608     }
609   } // end of loop i1
610   
611   
612   //Now we either add current events to stack or remove
613   //If no photons in current event - no need to add it to mixed
614   if(fPHOSEvent2->GetEntriesFast()>0){
615     prevPHOS->AddFirst(fPHOSEvent2) ;
616     fPHOSEvent2=0;
617     delete fPHOSEvent1;
618     fPHOSEvent1=0;
619     if(prevPHOS->GetSize()>100){//Remove redundant events
620       TClonesArray * tmp = static_cast<TClonesArray*>(prevPHOS->Last()) ;
621       prevPHOS->RemoveLast() ;
622       delete tmp ;
623     }
624   }
625   // Post output data.
626   PostData(1, fOutputContainer);
627   fEventCounter++;
628 }
629
630 //________________________________________________________________________
631 void AliAnalysisTaskPi0DiffEfficiency::Terminate(Option_t *)
632 {
633   // Draw result to the screen
634   // Called once at the end of the query
635   
636 }
637
638 //________________________________________________________________________
639 Bool_t AliAnalysisTaskPi0DiffEfficiency::IsGoodChannel(const char * det, Int_t mod, Int_t ix, Int_t iz)
640 {
641   //Check if this channel belogs to the good ones
642
643   if(strcmp(det,"PHOS")==0){
644     if(mod>5 || mod<1){
645       AliError(Form("No bad map for PHOS module %d ",mod)) ;
646       return kTRUE ;
647     }
648     if(!fPHOSBadMap[mod]){
649       AliError(Form("No Bad map for PHOS module %d",mod)) ;
650       return kTRUE ;
651     }
652     if(fPHOSBadMap[mod]->GetBinContent(ix,iz)>0)
653       return kFALSE ;
654     else
655       return kTRUE ;
656   }
657   else{
658     AliError(Form("Can not find bad channels for detector %s ",det)) ;
659   }
660   return kTRUE ;
661 }
662 //_____________________________________________________________________________
663 void AliAnalysisTaskPi0DiffEfficiency::FillHistogram(const char * key,Double_t x)const{
664   //FillHistogram
665   TH1I * tmpI = dynamic_cast<TH1I*>(fOutputContainer->FindObject(key)) ;
666   if(tmpI){
667     tmpI->Fill(x) ;
668     return ;
669   }
670   TH1F * tmpF = dynamic_cast<TH1F*>(fOutputContainer->FindObject(key)) ;
671   if(tmpF){
672     tmpF->Fill(x) ;
673     return ;
674   }
675   TH1D * tmpD = dynamic_cast<TH1D*>(fOutputContainer->FindObject(key)) ;
676   if(tmpD){
677     tmpD->Fill(x) ;
678     return ;
679   }
680   AliInfo(Form("can not find histogram <%s> ",key)) ;
681 }
682 //_____________________________________________________________________________
683 void AliAnalysisTaskPi0DiffEfficiency::FillHistogram(const char * key,Double_t x,Double_t y)const{
684   //FillHistogram
685   TObject * tmp = fOutputContainer->FindObject(key) ;
686   if(!tmp){
687     AliInfo(Form("can not find histogram <%s> ",key)) ;
688     return ;
689   }
690   if(tmp->IsA() == TClass::GetClass("TH1F")){
691     ((TH1F*)tmp)->Fill(x,y) ;
692     return ;
693   }
694   if(tmp->IsA() == TClass::GetClass("TH2F")){
695     ((TH2F*)tmp)->Fill(x,y) ;
696     return ;
697   }
698   AliError(Form("Calling FillHistogram with 2 parameters for histo <%s> of type %s",key,tmp->IsA()->GetName())) ;
699 }
700
701 //_____________________________________________________________________________
702 void AliAnalysisTaskPi0DiffEfficiency::FillHistogram(const char * key,Double_t x,Double_t y, Double_t z) const{
703   //Fills 1D histograms with key
704   TObject * tmp = fOutputContainer->FindObject(key) ;
705   if(!tmp){
706     AliInfo(Form("can not find histogram <%s> ",key)) ;
707     return ;
708   }
709   if(tmp->IsA() == TClass::GetClass("TH2F")){
710     ((TH2F*)tmp)->Fill(x,y,z) ;
711     return ;
712   }
713   if(tmp->IsA() == TClass::GetClass("TH3F")){
714     ((TH3F*)tmp)->Fill(x,y,z) ;
715     return ;
716   }
717 }
718 //_____________________________________________________________________________
719 Bool_t AliAnalysisTaskPi0DiffEfficiency::TestLambda(Double_t l1,Double_t l2){
720   Double_t l1Mean=1.22 ;
721   Double_t l2Mean=2.0 ;
722   Double_t l1Sigma=0.42 ;
723   Double_t l2Sigma=0.71 ;
724   Double_t c=-0.59 ;
725   Double_t R2=(l1-l1Mean)*(l1-l1Mean)/l1Sigma/l1Sigma+(l2-l2Mean)*(l2-l2Mean)/l2Sigma/l2Sigma-c*(l1-l1Mean)*(l2-l2Mean)/l1Sigma/l2Sigma ;
726   return (R2<9.) ;
727   
728 }
729 //___________________________________________________________________________
730 void AliAnalysisTaskPi0DiffEfficiency::ProcessMC(){
731   //fill histograms for efficiensy etc. calculation
732   const Double_t rcut = 1. ; //cut for primary particles
733   //---------First pi0/eta-----------------------------
734   char partName[10] ;
735   char hkey[55] ;
736
737   AliAODEvent *event = dynamic_cast<AliAODEvent*>(InputEvent());
738   TClonesArray *mcArray = (TClonesArray*)event->FindListObject(AliAODMCParticle::StdBranchName());
739   for(Int_t i=0;i<mcArray->GetEntriesFast();i++){
740      AliAODMCParticle* particle =  (AliAODMCParticle*) mcArray->At(i);
741     if(particle->GetPdgCode() == 111)
742       sprintf(partName,"pi0") ;
743     else
744       if(particle->GetPdgCode() == 221)
745         sprintf(partName,"eta") ;
746       else
747         if(particle->GetPdgCode() == 22)
748            sprintf(partName,"gamma") ;
749         else
750            continue ;
751
752     //Primary particle
753     Double_t r=TMath::Sqrt(particle->Xv()*particle->Xv()+particle->Yv()*particle->Yv());
754     if(r >rcut)
755       continue ;
756
757     Double_t pt = particle->Pt() ;
758     //Total number of pi0 with creation radius <1 cm
759     sprintf(hkey,"hMC_all_%s_cen%d",partName,fCenBin) ;
760     FillHistogram(hkey,pt) ;
761     if(TMath::Abs(particle->Y())<0.12){
762       sprintf(hkey,"hMC_unitEta_%s_cen%d",partName,fCenBin) ;
763       FillHistogram(hkey,pt) ;
764     }
765
766     sprintf(hkey,"hMC_rap_%s_cen%d",partName,fCenBin) ;
767     FillHistogram(hkey,particle->Y()) ;
768     
769     Double_t phi=particle->Phi() ;
770     while(phi<0.)phi+=TMath::TwoPi() ;
771     while(phi>TMath::TwoPi())phi-=TMath::TwoPi() ;
772     sprintf(hkey,"hMC_phi_%s_cen%d",partName,fCenBin) ;
773     FillHistogram(hkey,phi) ;
774    
775
776     //Check if one of photons converted
777     if(particle->GetNDaughters()!=2)
778      continue ; //Do not account Dalitz decays
779
780 /*
781     TParticle * gamma1 = fStack->Particle(particle->GetFirstDaughter());
782     TParticle * gamma2 = fStack->Particle(particle->GetLastDaughter());
783     //Number of pi0s decayed into acceptance
784     Int_t mod1,mod2 ;
785     Double_t x=0.,z=0. ;
786     Bool_t hitPHOS1 = fPHOSGeo->ImpactOnEmc(gamma1, mod1, z,x) ;
787     Bool_t hitPHOS2 = fPHOSGeo->ImpactOnEmc(gamma2, mod2, z,x) ;
788
789     Bool_t goodPair=kFALSE ;
790     if( hitPHOS1 && hitPHOS2){
791       sprintf(hkey,"hMC_PHOSacc_%s",partName) ;
792       FillHistogram(hkey,pt) ;
793       goodPair=kTRUE ;
794     }
795
796 */
797   }
798  
799   //Now calculate "Real" distribution of clusters with primary
800   TClonesArray cluPrim("AliCaloPhoton",200) ; //clusters with primary
801   Int_t multClust = event->GetNumberOfCaloClusters();
802   Int_t inPHOS=0 ;
803   Double_t vtx0[3] = {0,0,0}; 
804   for (Int_t i=0; i<multClust; i++) {
805     AliVCluster *clu = event->GetCaloCluster(i);
806     if ( !clu->IsPHOS() || clu->E()<0.3) continue;
807     if(clu->GetLabel()<0) continue ;
808
809     Float_t  position[3];
810     clu->GetPosition(position);
811     TVector3 global(position) ;
812     Int_t relId[4] ;
813     fPHOSGeo->GlobalPos2RelId(global,relId) ;
814     Int_t mod  = relId[0] ;
815     Int_t cellX = relId[2];
816     Int_t cellZ = relId[3] ;
817     if ( !IsGoodChannel("PHOS",mod,cellX,cellZ) ) 
818       continue ;
819     if(clu->GetNCells()<3)
820       continue ;
821
822     TLorentzVector pv1 ;
823     clu->GetMomentum(pv1 ,vtx0);
824     
825     if(inPHOS>=cluPrim.GetSize()){
826       cluPrim.Expand(inPHOS+50) ;
827     }
828     AliCaloPhoton * ph = new(cluPrim[inPHOS]) AliCaloPhoton(pv1.X(),pv1.Py(),pv1.Z(),pv1.E()) ;
829     //AliCaloPhoton * ph = (AliCaloPhoton*)fPHOSEvent->At(inPHOS) ;
830     ph->SetModule(mod) ;
831     ph->SetMomV2(&pv1) ;
832     ph->SetNCells(clu->GetNCells());
833     ph->SetDispBit(TestLambda(clu->GetM20(),clu->GetM02())) ;
834     ph->SetCPVBit(clu->GetEmcCpvDistance()>1.) ;
835
836     inPHOS++ ;
837
838   }
839   
840   //Single photon
841   char key[55] ;
842   for (Int_t i1=0; i1<inPHOS; i1++) {
843     AliCaloPhoton * ph1=(AliCaloPhoton*)cluPrim.At(i1) ;
844     snprintf(key,55,"hMCPhotAll_cen%d",fCenBin) ;
845     FillHistogram(key,ph1->Pt()) ;
846     if(ph1->IsCPVOK() ){
847       snprintf(key,55,"hMCPhotCPV_cen%d",fCenBin) ;
848       FillHistogram(key,ph1->Pt()) ;
849       
850     }
851     if(ph1->IsDispOK()){
852       snprintf(key,55,"hMCPhotDisp_cen%d",fCenBin) ;
853       FillHistogram(key,ph1->Pt()) ;
854       if(ph1->IsCPVOK()){
855         snprintf(key,55,"hMCPhotBoth_cen%d",fCenBin) ;
856         FillHistogram(key,ph1->Pt()) ;
857       }
858     } // end of loop i2
859   } // end of loop i1 
860
861   // Fill Real disribution
862   for (Int_t i1=0; i1<inPHOS-1; i1++) {
863     AliCaloPhoton * ph1=(AliCaloPhoton*)cluPrim.At(i1) ;
864     for (Int_t i2=i1+1; i2<inPHOS; i2++) {
865       AliCaloPhoton * ph2=(AliCaloPhoton*)cluPrim.At(i2) ;
866       TLorentzVector p12  = *ph1  + *ph2;
867       TLorentzVector pv12 = *(ph1->GetMomV2()) + *(ph2->GetMomV2());      
868   
869       snprintf(key,55,"hMCMassPtAll_cen%d",fCenBin) ;
870       FillHistogram(key,p12.M() ,p12.Pt()) ;
871       if(ph1->IsCPVOK() && ph2->IsCPVOK()){
872         snprintf(key,55,"hMCMassPtCPV_cen%d",fCenBin) ;
873         FillHistogram(key,p12.M() ,p12.Pt()) ;
874       }
875       if(ph1->IsDispOK() && ph2->IsDispOK()){
876         snprintf(key,55,"hMCMassPtDisp_cen%d",fCenBin) ;
877         FillHistogram(key,p12.M() ,p12.Pt()) ;
878
879         if(ph1->IsCPVOK() && ph2->IsCPVOK()){
880           snprintf(key,55,"hMCMassPtBoth_cen%d",fCenBin) ;
881           FillHistogram(key,p12.M() ,p12.Pt()) ;
882         }
883       }
884     } // end of loop i2
885   } // end of loop i1
886
887   
888 }
889 //___________________________________________________________________________
890 Bool_t AliAnalysisTaskPi0DiffEfficiency::IsSameCluster(AliAODCaloCluster * c1, AliAODCaloCluster * c2)const{
891  //Compare clusters before and after embedding
892  //clusters are the same if 
893  // - Energy changed less than 0.1%  (numerical accuracy in reconstruction)
894  // - lists of digits are the same
895   
896  if(c1->GetNCells() != c2->GetNCells())
897    return kFALSE ;
898  
899  if(TMath::Abs(c1->E()-c2->E())>0.001*c1->E())
900    return kFALSE ;
901
902  UShort_t *list1 = c1->GetCellsAbsId() ; 
903  UShort_t *list2 = c2->GetCellsAbsId() ; 
904  for(Int_t i=0; i<c1->GetNCells(); i++){
905   if(list1[i] != list2[i])
906     return kFALSE ;
907  }
908  return kTRUE ; 
909   
910 }
911
912
913