Event embedding tasks for PHOS are added (D.Peressounko)
[u/mrichter/AliRoot.git] / PWG4 / UserTasks / PHOS_embedding / AliAnalysisTaskPi0Efficiency.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 "AliAnalysisTaskPi0Efficiency.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(AliAnalysisTaskPi0Efficiency)
40
41 //________________________________________________________________________
42 AliAnalysisTaskPi0Efficiency::AliAnalysisTaskPi0Efficiency(const char *name) 
43 : AliAnalysisTaskSE(name),
44   fStack(0),
45   fOutputContainer(0),
46   fPHOSEvent(0),
47   fPHOSCalibData(0),
48   fNonLinCorr(0),
49   fRPfull(0),
50   fRPA(0),
51   fRPC(0),
52   fRPFar(0),
53   fRPAFar(0),
54   fRPCFar(0),
55   fCentrality(0),
56   fCenBin(0),
57   fPHOSGeo(0),
58   fEventCounter(0)
59 {
60   // Constructor
61   for(Int_t i=0;i<1;i++){
62     for(Int_t j=0;j<10;j++)
63       for(Int_t k=0;k<11;k++)
64         fPHOSEvents[i][j][k]=0 ;
65   }
66   
67   // Output slots #0 write into a TH1 container
68   DefineOutput(1,TList::Class());
69
70   // Set bad channel map
71   char key[55] ;
72   for(Int_t i=0; i<6; i++){
73     sprintf(key,"PHOS_BadMap_mod%d",i) ;
74     fPHOSBadMap[i]=new TH2I(key,"Bad Modules map",64,0.,64.,56,0.,56.) ;
75   }
76   // Initialize the PHOS geometry
77   fPHOSGeo = AliPHOSGeometry::GetInstance("IHEP") ;
78
79
80 }
81
82 //________________________________________________________________________
83 void AliAnalysisTaskPi0Efficiency::UserCreateOutputObjects()
84 {
85   // Create histograms
86   // Called once
87
88   // ESD histograms
89   if(fOutputContainer != NULL){
90     delete fOutputContainer;
91   }
92   fOutputContainer = new TList();
93   fOutputContainer->SetOwner(kTRUE);
94
95   //Event selection
96   fOutputContainer->Add(new TH1F("hSelEvents","Event celection", 10,0.,10.)) ;
97
98   //vertex distribution
99   fOutputContainer->Add(new TH1F("hZvertex","Z vertex position", 50,-25.,25.)) ;
100
101   //Centrality
102   fOutputContainer->Add(new TH1F("hCentrality","Event centrality", 100,0.,100.)) ;
103
104   //QA histograms                       
105   fOutputContainer->Add(new TH2F("hCluM1","Cell (X,Z), M1" ,64,0.5,64.5, 56,0.5,56.5));
106   fOutputContainer->Add(new TH2F("hCluM2","Cell (X,Z), M2" ,64,0.5,64.5, 56,0.5,56.5));
107   fOutputContainer->Add(new TH2F("hCluM3","Cell (X,Z), M3" ,64,0.5,64.5, 56,0.5,56.5));
108
109   Int_t nM       = 500;
110   Double_t mMin  = 0.0;
111   Double_t mMax  = 1.0;
112   Int_t nPt      = 200;
113   Double_t ptMin = 0;
114   Double_t ptMax = 20;
115
116   char key[55] ;
117   for(Int_t cent=0; cent<6; cent++){
118     //Single photon
119     snprintf(key,55,"hPhotAll_cen%d",cent) ;
120     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
121     snprintf(key,55,"hPhotCPV_cen%d",cent) ;
122     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
123     snprintf(key,55,"hPhotDisp_cen%d",cent) ;
124     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
125     snprintf(key,55,"hPhotBoth_cen%d",cent) ;
126     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
127     
128     snprintf(key,55,"hMassPtAll_cen%d",cent) ;
129     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
130     snprintf(key,55,"hMassPtCPV_cen%d",cent) ;
131     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
132     snprintf(key,55,"hMassPtDisp_cen%d",cent) ;
133     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
134     snprintf(key,55,"hMassPtBoth_cen%d",cent) ;
135     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
136     
137     //Mixed
138     snprintf(key,55,"hMiMassPtAll_cen%d",cent) ;
139     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
140     snprintf(key,55,"hMiMassPtCPV_cen%d",cent) ;
141     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
142     snprintf(key,55,"hMiMassPtDisp_cen%d",cent) ;
143     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
144     snprintf(key,55,"hMiMassPtBoth_cen%d",cent) ;
145     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
146
147     snprintf(key,55,"hMCMassPtAll_cen%d",cent) ;
148     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
149     snprintf(key,55,"hMCMassPtCPV_cen%d",cent) ;
150     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
151     snprintf(key,55,"hMCMassPtDisp_cen%d",cent) ;
152     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
153     snprintf(key,55,"hMCMassPtBoth_cen%d",cent) ;
154     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
155
156     //Single photon
157     snprintf(key,55,"hMCPhotAll_cen%d",cent) ;
158     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
159     snprintf(key,55,"hMCPhotCPV_cen%d",cent) ;
160     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
161     snprintf(key,55,"hMCPhotDisp_cen%d",cent) ;
162     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
163     snprintf(key,55,"hMCPhotBoth_cen%d",cent) ;
164     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
165
166     
167     
168   }
169
170
171   //MC
172   for(Int_t cent=0; cent<6; cent++){
173     snprintf(key,55,"hMC_rap_gamma_cen%d",cent) ;
174     fOutputContainer->Add(new TH1F(key,"Rapidity pi0",200,-1.,1.)) ;
175     snprintf(key,55,"hMC_rap_pi0_cen%d",cent) ;
176     fOutputContainer->Add(new TH1F(key,"Rapidity pi0",200,-1.,1.)) ;
177     snprintf(key,55,"hMC_rap_eta_cen%d",cent) ;
178     fOutputContainer->Add(new TH1F("hMC_rap_eta","Rapidity eta",200,-1.,1.)) ;
179     snprintf(key,55,"hMC_phi_gamma_cen%d",cent) ;
180     fOutputContainer->Add(new TH1F(key,"Phi pi0",200,0.,TMath::TwoPi())) ;
181     snprintf(key,55,"hMC_phi_pi0_cen%d",cent) ;
182     fOutputContainer->Add(new TH1F(key,"Phi pi0",200,0.,TMath::TwoPi())) ;
183     snprintf(key,55,"hMC_phi_eta_cen%d",cent) ;
184     fOutputContainer->Add(new TH1F(key,"Phi eta",200,0.,TMath::TwoPi())) ;
185     snprintf(key,55,"hMC_all_gamma_cen%d",cent) ;
186     fOutputContainer->Add(new TH1F(key,"Rapidity photon",250,0.,25.)) ;
187     snprintf(key,55,"hMC_all_pi0_cen%d",cent) ;
188     fOutputContainer->Add(new TH1F(key,"Rapidity pi0",250,0.,25.)) ;
189     snprintf(key,55,"hMC_all_eta_cen%d",cent) ;
190     fOutputContainer->Add(new TH1F(key,"Rapidity eta",250,0.,25.)) ;
191     snprintf(key,55,"hMC_unitEta_gamma_cen%d",cent) ;
192     fOutputContainer->Add(new TH1F(key,"Pt photon",250,0.,25.)) ;
193     snprintf(key,55,"hMC_unitEta_pi0_cen%d",cent) ;
194     fOutputContainer->Add(new TH1F(key,"Rapidity eta",250,0.,25.)) ;
195     snprintf(key,55,"hMC_unitEta_eta_cen%d",cent) ;
196     fOutputContainer->Add(new TH1F(key,"Rapidity eta",250,0.,25.)) ;
197   }
198   
199   PostData(1, fOutputContainer);
200
201 }
202
203 //________________________________________________________________________
204 void AliAnalysisTaskPi0Efficiency::UserExec(Option_t *) 
205 {
206   // Main loop, called for each event
207   // Analyze ESD/AOD
208
209   FillHistogram("hSelEvents",0.5) ;  
210   
211   AliAODEvent *event = dynamic_cast<AliAODEvent*>(InputEvent());
212   if (!event) {
213     Printf("ERROR: Could not retrieve event");
214     PostData(1, fOutputContainer);
215     return;
216   }
217   
218   FillHistogram("hSelEvents",1.5) ;
219   AliAODHeader *header = event->GetHeader() ;
220   
221   // Checks if we have a primary vertex
222   // Get primary vertices form ESD
223   const AliAODVertex *esdVertex5 = event->GetPrimaryVertex();
224
225  // don't rely on ESD vertex, assume (0,0,0)
226   Double_t vtx0[3] ={0.,0.,0.};
227   Double_t vtx5[3] ={0.,0.,0.};
228   
229   vtx5[0] = esdVertex5->GetX();
230   vtx5[1] = esdVertex5->GetY();
231   vtx5[2] = esdVertex5->GetZ();
232   
233   
234   FillHistogram("hZvertex",esdVertex5->GetZ());
235   if (TMath::Abs(esdVertex5->GetZ()) > 10. ){
236     PostData(1, fOutputContainer);
237     return;
238   }
239   FillHistogram("hSelEvents",2.5) ;
240
241   //Vtx class z-bin
242   //  Int_t zvtx = (Int_t)((vtx5[2]+10.)/2.) ;
243   //  if(zvtx<0)zvtx=0 ;
244   //  if(zvtx>9)zvtx=9 ;
245   Int_t zvtx=0 ;
246
247 //  fCentrality=header->GetCentralityP()->GetCentralityPercentile("V0M"); // returns the centrality percentile, 
248 //                                                          //a float from 0 to 100 (or to the trigger efficiency)
249    fCentrality=header->GetZDCN2Energy() ;
250
251   if( fCentrality < 0. || fCentrality>80.){
252     PostData(1, fOutputContainer);
253     return;
254   }
255   FillHistogram("hSelEvents",3.5) ;
256   Float_t bins[7]={0.,5.,10.,20.,40.,60.,80.} ;
257   fCenBin=0 ;
258   while(fCenBin<6 && fCentrality > bins[fCenBin+1])
259     fCenBin++ ; 
260
261  
262   //reaction plain
263   fRPfull= header->GetZDCN1Energy() ;
264   if(fRPfull==999){ //reaction plain was not defined
265     PostData(1, fOutputContainer);
266     return;
267   } 
268
269   FillHistogram("hSelEvents",4.5) ;
270   //All event selections done
271   FillHistogram("hCentrality",fCentrality) ;
272   //Reaction plain is defined in the range (-pi/2;pi/2)
273   //We have 10 bins
274   Int_t irp=Int_t(10.*(fRPfull+TMath::PiOver2())/TMath::Pi());
275   if(irp>9)irp=9 ;
276
277   if(!fPHOSEvents[zvtx][fCenBin][irp]) 
278     fPHOSEvents[zvtx][fCenBin][irp]=new TList() ;
279   TList * prevPHOS = fPHOSEvents[zvtx][fCenBin][irp] ;
280
281   // Get PHOS rotation matrices from ESD and set them to the PHOS geometry
282   if(fEventCounter == 0) {
283     for(Int_t mod=0; mod<5; mod++) {
284       const TGeoHMatrix* m =header->GetPHOSMatrix(mod) ;
285       fPHOSGeo->SetMisalMatrix(m,mod) ;
286       Printf("PHOS geo matrix for module # %d is set: %p\n", mod,m);
287     }
288     fEventCounter++ ;
289   }
290
291   ProcessMC() ;
292
293   if(fPHOSEvent)
294     fPHOSEvent->Clear() ;
295   else
296     fPHOSEvent = new TClonesArray("AliCaloPhoton",200) ;
297
298
299   char key[55] ;
300   Int_t inPHOS=0 ;
301   TVector3 vertex(vtx0);
302   TClonesArray * clusters = (TClonesArray*)event->FindListObject("EmbeddedCaloClusters") ;
303   Int_t multClust = clusters->GetEntriesFast();
304   for (Int_t i=0; i<multClust; i++) {
305     AliVCluster *clu = (AliVCluster*) clusters->At(i);
306     if ( !clu->IsPHOS() || clu->E()<0.3) continue;
307
308     Float_t  position[3];
309     clu->GetPosition(position);
310     TVector3 global(position) ;
311     Int_t relId[4] ;
312     fPHOSGeo->GlobalPos2RelId(global,relId) ;
313     Int_t mod  = relId[0] ;
314     Int_t cellX = relId[2];
315     Int_t cellZ = relId[3] ;
316     if ( !IsGoodChannel("PHOS",mod,cellX,cellZ) ) 
317       continue ;
318     if(clu->GetNCells()<3)
319       continue ;
320
321     sprintf(key,"hCluM%d",mod) ;
322     FillHistogram(key,cellX,cellZ,1.);
323
324     TLorentzVector pv1 ;
325     clu->GetMomentum(pv1 ,vtx0);
326     
327     if(inPHOS>=fPHOSEvent->GetSize()){
328       fPHOSEvent->Expand(inPHOS+50) ;
329     }
330     new((*fPHOSEvent)[inPHOS]) AliCaloPhoton(pv1.X(),pv1.Py(),pv1.Z(),pv1.E()) ;
331     AliCaloPhoton * ph = (AliCaloPhoton*)fPHOSEvent->At(inPHOS) ;
332     ph->SetModule(mod) ;
333     ph->SetMomV2(&pv1) ;
334     ph->SetNCells(clu->GetNCells());
335     ph->SetDispBit(TestLambda(clu->GetM20(),clu->GetM02())) ;
336     ph->SetCPVBit(clu->GetEmcCpvDistance()>1.) ;
337
338     inPHOS++ ;
339   }
340   //Single photon
341   for (Int_t i1=0; i1<inPHOS; i1++) {
342     AliCaloPhoton * ph1=(AliCaloPhoton*)fPHOSEvent->At(i1) ;
343     snprintf(key,55,"hPhotAll_cen%d",fCenBin) ;
344     FillHistogram(key,ph1->Pt()) ;
345     if(ph1->IsCPVOK() ){
346       snprintf(key,55,"hPhotCPV_cen%d",fCenBin) ;
347       FillHistogram(key,ph1->Pt()) ;
348     }
349     if(ph1->IsDispOK()){
350       snprintf(key,55,"hPhotDisp_cen%d",fCenBin) ;
351       FillHistogram(key,ph1->Pt()) ;
352       if(ph1->IsCPVOK()){
353         snprintf(key,55,"hPhotBoth_cen%d",fCenBin) ;
354         FillHistogram(key,ph1->Pt()) ;
355       }
356     } // end of loop i2
357   } // end of loop i1 
358
359   // Fill Real disribution
360   for (Int_t i1=0; i1<inPHOS-1; i1++) {
361     AliCaloPhoton * ph1=(AliCaloPhoton*)fPHOSEvent->At(i1) ;
362     for (Int_t i2=i1+1; i2<inPHOS; i2++) {
363       AliCaloPhoton * ph2=(AliCaloPhoton*)fPHOSEvent->At(i2) ;
364       TLorentzVector p12  = *ph1  + *ph2;
365       TLorentzVector pv12 = *(ph1->GetMomV2()) + *(ph2->GetMomV2());      
366
367       snprintf(key,55,"hMassPtAll_cen%d",fCenBin) ;
368       FillHistogram(key,p12.M() ,p12.Pt()) ;
369       if(ph1->IsCPVOK() && ph2->IsCPVOK()){
370         snprintf(key,55,"hMassPtCPV_cen%d",fCenBin) ;
371         FillHistogram(key,p12.M() ,p12.Pt()) ;
372       }
373       if(ph1->IsDispOK() && ph2->IsDispOK()){
374         snprintf(key,55,"hMassPtDisp_cen%d",fCenBin) ;
375         FillHistogram(key,p12.M() ,p12.Pt()) ;
376
377         if(ph1->IsCPVOK() && ph2->IsCPVOK()){
378           snprintf(key,55,"hMassPtBoth_cen%d",fCenBin) ;
379           FillHistogram(key,p12.M() ,p12.Pt()) ;
380         }
381       }
382     } // end of loop i2
383   } // end of loop i1 
384
385   //now mixed
386   for (Int_t i1=0; i1<inPHOS; i1++) {
387     AliCaloPhoton * ph1=(AliCaloPhoton*)fPHOSEvent->At(i1) ;
388     for(Int_t ev=0; ev<prevPHOS->GetSize();ev++){
389       TClonesArray * mixPHOS = static_cast<TClonesArray*>(prevPHOS->At(ev)) ;
390       for(Int_t i2=0; i2<mixPHOS->GetEntriesFast();i2++){
391         AliCaloPhoton * ph2=(AliCaloPhoton*)mixPHOS->At(i2) ;
392         TLorentzVector p12  = *ph1  + *ph2;
393         TLorentzVector pv12 = *(ph1->GetMomV2()) + *(ph2->GetMomV2());
394         
395         snprintf(key,55,"hMiMassPtAll_cen%d",fCenBin) ;
396         FillHistogram(key,p12.M() ,p12.Pt()) ;
397         if(ph1->IsCPVOK() && ph2->IsCPVOK()){
398           snprintf(key,55,"hMiMassPtCPV_cen%d",fCenBin) ;
399           FillHistogram(key,p12.M() ,p12.Pt()) ;
400         }
401         if(ph1->IsDispOK() && ph2->IsDispOK()){
402           snprintf(key,55,"hMiMassPtDisp_cen%d",fCenBin) ;
403           FillHistogram(key,p12.M() ,p12.Pt()) ;
404           
405           if(ph1->IsCPVOK() && ph2->IsCPVOK()){
406             snprintf(key,55,"hMiMassPtBoth_cen%d",fCenBin) ;
407             FillHistogram(key,p12.M() ,p12.Pt()) ;
408           }
409         }
410       } // end of loop i2
411     }
412   } // end of loop i1
413   
414   
415   //Now we either add current events to stack or remove
416   //If no photons in current event - no need to add it to mixed
417   if(fPHOSEvent->GetEntriesFast()>0){
418     prevPHOS->AddFirst(fPHOSEvent) ;
419     fPHOSEvent=0;
420     if(prevPHOS->GetSize()>100){//Remove redundant events
421       TClonesArray * tmp = static_cast<TClonesArray*>(prevPHOS->Last()) ;
422       prevPHOS->RemoveLast() ;
423       delete tmp ;
424     }
425   }
426   // Post output data.
427   PostData(1, fOutputContainer);
428   fEventCounter++;
429 }
430
431 //________________________________________________________________________
432 void AliAnalysisTaskPi0Efficiency::Terminate(Option_t *)
433 {
434   // Draw result to the screen
435   // Called once at the end of the query
436   
437 }
438
439 //________________________________________________________________________
440 Bool_t AliAnalysisTaskPi0Efficiency::IsGoodChannel(const char * det, Int_t mod, Int_t ix, Int_t iz)
441 {
442   //Check if this channel belogs to the good ones
443
444   if(strcmp(det,"PHOS")==0){
445     if(mod>5 || mod<1){
446       AliError(Form("No bad map for PHOS module %d ",mod)) ;
447       return kTRUE ;
448     }
449     if(!fPHOSBadMap[mod]){
450       AliError(Form("No Bad map for PHOS module %d",mod)) ;
451       return kTRUE ;
452     }
453     if(fPHOSBadMap[mod]->GetBinContent(ix,iz)>0)
454       return kFALSE ;
455     else
456       return kTRUE ;
457   }
458   else{
459     AliError(Form("Can not find bad channels for detector %s ",det)) ;
460   }
461   return kTRUE ;
462 }
463 //_____________________________________________________________________________
464 void AliAnalysisTaskPi0Efficiency::FillHistogram(const char * key,Double_t x)const{
465   //FillHistogram
466   TH1I * tmpI = dynamic_cast<TH1I*>(fOutputContainer->FindObject(key)) ;
467   if(tmpI){
468     tmpI->Fill(x) ;
469     return ;
470   }
471   TH1F * tmpF = dynamic_cast<TH1F*>(fOutputContainer->FindObject(key)) ;
472   if(tmpF){
473     tmpF->Fill(x) ;
474     return ;
475   }
476   TH1D * tmpD = dynamic_cast<TH1D*>(fOutputContainer->FindObject(key)) ;
477   if(tmpD){
478     tmpD->Fill(x) ;
479     return ;
480   }
481   AliInfo(Form("can not find histogram <%s> ",key)) ;
482 }
483 //_____________________________________________________________________________
484 void AliAnalysisTaskPi0Efficiency::FillHistogram(const char * key,Double_t x,Double_t y)const{
485   //FillHistogram
486   TObject * tmp = fOutputContainer->FindObject(key) ;
487   if(!tmp){
488     AliInfo(Form("can not find histogram <%s> ",key)) ;
489     return ;
490   }
491   if(tmp->IsA() == TClass::GetClass("TH1F")){
492     ((TH1F*)tmp)->Fill(x,y) ;
493     return ;
494   }
495   if(tmp->IsA() == TClass::GetClass("TH2F")){
496     ((TH2F*)tmp)->Fill(x,y) ;
497     return ;
498   }
499   AliError(Form("Calling FillHistogram with 2 parameters for histo <%s> of type %s",key,tmp->IsA()->GetName())) ;
500 }
501
502 //_____________________________________________________________________________
503 void AliAnalysisTaskPi0Efficiency::FillHistogram(const char * key,Double_t x,Double_t y, Double_t z) const{
504   //Fills 1D histograms with key
505   TObject * tmp = fOutputContainer->FindObject(key) ;
506   if(!tmp){
507     AliInfo(Form("can not find histogram <%s> ",key)) ;
508     return ;
509   }
510   if(tmp->IsA() == TClass::GetClass("TH2F")){
511     ((TH2F*)tmp)->Fill(x,y,z) ;
512     return ;
513   }
514   if(tmp->IsA() == TClass::GetClass("TH3F")){
515     ((TH3F*)tmp)->Fill(x,y,z) ;
516     return ;
517   }
518 }
519 //_____________________________________________________________________________
520 Bool_t AliAnalysisTaskPi0Efficiency::TestLambda(Double_t l1,Double_t l2){
521   Double_t l1Mean=1.22 ;
522   Double_t l2Mean=2.0 ;
523   Double_t l1Sigma=0.42 ;
524   Double_t l2Sigma=0.71 ;
525   Double_t c=-0.59 ;
526   Double_t R2=(l1-l1Mean)*(l1-l1Mean)/l1Sigma/l1Sigma+(l2-l2Mean)*(l2-l2Mean)/l2Sigma/l2Sigma-c*(l1-l1Mean)*(l2-l2Mean)/l1Sigma/l2Sigma ;
527   return (R2<9.) ;
528   
529 }
530 //___________________________________________________________________________
531 void AliAnalysisTaskPi0Efficiency::ProcessMC(){
532   //fill histograms for efficiensy etc. calculation
533   const Double_t rcut = 1. ; //cut for primary particles
534   //---------First pi0/eta-----------------------------
535   char partName[10] ;
536   char hkey[55] ;
537
538   AliAODEvent *event = dynamic_cast<AliAODEvent*>(InputEvent());
539   TClonesArray *mcArray = (TClonesArray*)event->FindListObject(AliAODMCParticle::StdBranchName());
540   for(Int_t i=0;i<mcArray->GetEntriesFast();i++){
541      AliAODMCParticle* particle =  (AliAODMCParticle*) mcArray->At(i);
542     if(particle->GetPdgCode() == 111)
543       sprintf(partName,"pi0") ;
544     else
545       if(particle->GetPdgCode() == 221)
546         sprintf(partName,"eta") ;
547       else
548         if(particle->GetPdgCode() == 22)
549            sprintf(partName,"gamma") ;
550         else
551            continue ;
552
553     //Primary particle
554     Double_t r=TMath::Sqrt(particle->Xv()*particle->Xv()+particle->Yv()*particle->Yv());
555     if(r >rcut)
556       continue ;
557
558     Double_t pt = particle->Pt() ;
559     //Total number of pi0 with creation radius <1 cm
560     sprintf(hkey,"hMC_all_%s_cen%d",partName,fCenBin) ;
561     FillHistogram(hkey,pt) ;
562     if(TMath::Abs(particle->Y())<0.12){
563       sprintf(hkey,"hMC_unitEta_%s_cen%d",partName,fCenBin) ;
564       FillHistogram(hkey,pt) ;
565     }
566
567     sprintf(hkey,"hMC_rap_%s_cen%d",partName,fCenBin) ;
568     FillHistogram(hkey,particle->Y()) ;
569     
570     Double_t phi=particle->Phi() ;
571     while(phi<0.)phi+=TMath::TwoPi() ;
572     while(phi>TMath::TwoPi())phi-=TMath::TwoPi() ;
573     sprintf(hkey,"hMC_phi_%s_cen%d",partName,fCenBin) ;
574     FillHistogram(hkey,phi) ;
575    
576
577     //Check if one of photons converted
578     if(particle->GetNDaughters()!=2)
579      continue ; //Do not account Dalitz decays
580
581 /*
582     TParticle * gamma1 = fStack->Particle(particle->GetFirstDaughter());
583     TParticle * gamma2 = fStack->Particle(particle->GetLastDaughter());
584     //Number of pi0s decayed into acceptance
585     Int_t mod1,mod2 ;
586     Double_t x=0.,z=0. ;
587     Bool_t hitPHOS1 = fPHOSGeo->ImpactOnEmc(gamma1, mod1, z,x) ;
588     Bool_t hitPHOS2 = fPHOSGeo->ImpactOnEmc(gamma2, mod2, z,x) ;
589
590     Bool_t goodPair=kFALSE ;
591     if( hitPHOS1 && hitPHOS2){
592       sprintf(hkey,"hMC_PHOSacc_%s",partName) ;
593       FillHistogram(hkey,pt) ;
594       goodPair=kTRUE ;
595     }
596
597 */
598   }
599  
600   //Now calculate "Real" distribution of clusters with primary
601   TClonesArray cluPrim("AliCaloPhoton",200) ; //clusters with primary
602   Int_t multClust = event->GetNumberOfCaloClusters();
603   Int_t inPHOS=0 ;
604   Double_t vtx0[3] = {0,0,0}; 
605   for (Int_t i=0; i<multClust; i++) {
606     AliVCluster *clu = event->GetCaloCluster(i);
607     if ( !clu->IsPHOS() || clu->E()<0.3) continue;
608     if(clu->GetLabel()<0) continue ;
609
610     Float_t  position[3];
611     clu->GetPosition(position);
612     TVector3 global(position) ;
613     Int_t relId[4] ;
614     fPHOSGeo->GlobalPos2RelId(global,relId) ;
615     Int_t mod  = relId[0] ;
616     Int_t cellX = relId[2];
617     Int_t cellZ = relId[3] ;
618     if ( !IsGoodChannel("PHOS",mod,cellX,cellZ) ) 
619       continue ;
620     if(clu->GetNCells()<3)
621       continue ;
622
623     TLorentzVector pv1 ;
624     clu->GetMomentum(pv1 ,vtx0);
625     
626     if(inPHOS>=cluPrim.GetSize()){
627       cluPrim.Expand(inPHOS+50) ;
628     }
629     AliCaloPhoton * ph = new(cluPrim[inPHOS]) AliCaloPhoton(pv1.X(),pv1.Py(),pv1.Z(),pv1.E()) ;
630     //AliCaloPhoton * ph = (AliCaloPhoton*)fPHOSEvent->At(inPHOS) ;
631     ph->SetModule(mod) ;
632     ph->SetMomV2(&pv1) ;
633     ph->SetNCells(clu->GetNCells());
634     ph->SetDispBit(TestLambda(clu->GetM20(),clu->GetM02())) ;
635     ph->SetCPVBit(clu->GetEmcCpvDistance()>1.) ;
636
637     inPHOS++ ;
638
639   }
640   
641   //Single photon
642   char key[55] ;
643   for (Int_t i1=0; i1<inPHOS; i1++) {
644     AliCaloPhoton * ph1=(AliCaloPhoton*)cluPrim.At(i1) ;
645     snprintf(key,55,"hMCPhotAll_cen%d",fCenBin) ;
646     FillHistogram(key,ph1->Pt()) ;
647     if(ph1->IsCPVOK() ){
648       snprintf(key,55,"hMCPhotCPV_cen%d",fCenBin) ;
649       FillHistogram(key,ph1->Pt()) ;
650       
651     }
652     if(ph1->IsDispOK()){
653       snprintf(key,55,"hMCPhotDisp_cen%d",fCenBin) ;
654       FillHistogram(key,ph1->Pt()) ;
655       if(ph1->IsCPVOK()){
656         snprintf(key,55,"hMCPhotBoth_cen%d",fCenBin) ;
657         FillHistogram(key,ph1->Pt()) ;
658       }
659     } // end of loop i2
660   } // end of loop i1 
661
662   // Fill Real disribution
663   for (Int_t i1=0; i1<inPHOS-1; i1++) {
664     AliCaloPhoton * ph1=(AliCaloPhoton*)cluPrim.At(i1) ;
665     for (Int_t i2=i1+1; i2<inPHOS; i2++) {
666       AliCaloPhoton * ph2=(AliCaloPhoton*)cluPrim.At(i2) ;
667       TLorentzVector p12  = *ph1  + *ph2;
668       TLorentzVector pv12 = *(ph1->GetMomV2()) + *(ph2->GetMomV2());      
669   
670       snprintf(key,55,"hMCMassPtAll_cen%d",fCenBin) ;
671       FillHistogram(key,p12.M() ,p12.Pt()) ;
672       if(ph1->IsCPVOK() && ph2->IsCPVOK()){
673         snprintf(key,55,"hMCMassPtCPV_cen%d",fCenBin) ;
674         FillHistogram(key,p12.M() ,p12.Pt()) ;
675       }
676       if(ph1->IsDispOK() && ph2->IsDispOK()){
677         snprintf(key,55,"hMCMassPtDisp_cen%d",fCenBin) ;
678         FillHistogram(key,p12.M() ,p12.Pt()) ;
679
680         if(ph1->IsCPVOK() && ph2->IsCPVOK()){
681           snprintf(key,55,"hMCMassPtBoth_cen%d",fCenBin) ;
682           FillHistogram(key,p12.M() ,p12.Pt()) ;
683         }
684       }
685     } // end of loop i2
686   } // end of loop i1
687
688   
689 }
690
691