remove mechanism to avoid dependencies on PHOS and EMCAL Utils
[u/mrichter/AliRoot.git] / PWG4 / PartCorrDep / AliAnaPi0.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15 /* $Id: $ */
16
17 //_________________________________________________________________________
18 // Class to collect two-photon invariant mass distributions for
19 // extractin raw pi0 yield.
20 //
21 //-- Author: Dmitri Peressounko (RRC "KI") 
22 //-- Adapted to PartCorr frame by Lamia Benhabib (SUBATECH)
23 //-- and Gustavo Conesa (INFN-Frascati)
24 //_________________________________________________________________________
25
26
27 // --- ROOT system ---
28 #include "TH3.h"
29 //#include "Riostream.h"
30 #include "TCanvas.h"
31 #include "TPad.h"
32 #include "TROOT.h"
33 #include "TClonesArray.h"
34 #include "TObjString.h"
35
36 //---- AliRoot system ----
37 #include "AliAnaPi0.h"
38 #include "AliCaloTrackReader.h"
39 #include "AliCaloPID.h"
40 #include "AliStack.h"
41 #include "AliFidutialCut.h"
42 #include "TParticle.h"
43 #include "AliAODCaloCluster.h"
44 #include "AliVEvent.h"
45 #include "AliPHOSGeoUtils.h"
46 #include "AliEMCALGeoUtils.h"
47
48 ClassImp(AliAnaPi0)
49
50 //________________________________________________________________________________________________________________________________________________  
51 AliAnaPi0::AliAnaPi0() : AliAnaPartCorrBaseClass(),
52 fNCentrBin(0),fNZvertBin(0),fNrpBin(0),
53 fNPID(0),fNmaxMixEv(0), fZvtxCut(0.),fCalorimeter(""),
54 fEMCALGeoName("EMCAL_COMPLETE"),fEventsList(0x0), fhEtalon(0x0),
55 fhRe1(0x0),fhMi1(0x0),fhRe2(0x0),fhMi2(0x0),fhRe3(0x0),fhMi3(0x0),fhEvents(0x0),
56 fhPrimPt(0x0), fhPrimAccPt(0x0), fhPrimY(0x0), fhPrimAccY(0x0), fhPrimPhi(0x0), fhPrimAccPhi(0x0), 
57 fPHOSGeo(0x0),fEMCALGeo(0x0)
58 {
59 //Default Ctor
60  InitParameters();
61  
62 }
63
64 //________________________________________________________________________________________________________________________________________________
65 AliAnaPi0::AliAnaPi0(const AliAnaPi0 & ex) : AliAnaPartCorrBaseClass(ex),  
66 fNCentrBin(ex.fNCentrBin),fNZvertBin(ex.fNZvertBin),fNrpBin(ex.fNrpBin),
67 fNPID(ex.fNPID),fNmaxMixEv(ex.fNmaxMixEv),fZvtxCut(ex.fZvtxCut), fCalorimeter(ex.fCalorimeter),
68 fEMCALGeoName(ex.fEMCALGeoName), fEventsList(ex.fEventsList), fhEtalon(ex.fhEtalon),
69 fhRe1(ex.fhRe1),fhMi1(ex.fhMi1),fhRe2(ex.fhRe2),fhMi2(ex.fhMi2),fhRe3(ex.fhRe3),fhMi3(ex.fhMi3),fhEvents(ex.fhEvents),
70 fhPrimPt(ex.fhPrimPt), fhPrimAccPt(ex.fhPrimAccPt), fhPrimY(ex.fhPrimY), 
71 fhPrimAccY(ex.fhPrimAccY), fhPrimPhi(ex.fhPrimPhi), fhPrimAccPhi(ex.fhPrimAccPhi),
72 fPHOSGeo(ex.fPHOSGeo),fEMCALGeo(ex.fEMCALGeo)
73 {
74   // cpy ctor
75   //Do not need it
76 }
77
78 //________________________________________________________________________________________________________________________________________________
79 AliAnaPi0 & AliAnaPi0::operator = (const AliAnaPi0 & ex)
80 {
81   // assignment operator
82   
83   if(this == &ex)return *this;
84   ((AliAnaPartCorrBaseClass *)this)->operator=(ex);
85   
86   fNCentrBin = ex.fNCentrBin  ; fNZvertBin = ex.fNZvertBin  ; fNrpBin = ex.fNrpBin  ; 
87   fNPID = ex.fNPID  ; fNmaxMixEv = ex.fNmaxMixEv  ; fZvtxCut = ex.fZvtxCut  ;  fCalorimeter = ex.fCalorimeter  ;  
88   fEMCALGeoName = ex.fEMCALGeoName ; fEventsList = ex.fEventsList  ;  fhEtalon = ex.fhEtalon  ; 
89   fhRe1 = ex.fhRe1  ; fhMi1 = ex.fhMi1  ; fhRe2 = ex.fhRe2  ; fhMi2 = ex.fhMi2  ; 
90   fhRe3 = ex.fhRe3  ; fhMi3 = ex.fhMi3  ; fhEvents = ex.fhEvents  ; 
91   fhPrimPt = ex.fhPrimPt  ;  fhPrimAccPt = ex.fhPrimAccPt  ;  fhPrimY = ex.fhPrimY  ;  
92   fhPrimAccY = ex.fhPrimAccY  ;  fhPrimPhi = ex.fhPrimPhi  ;  fhPrimAccPhi = ex.fhPrimAccPhi ;
93   
94   return *this;
95   
96 }
97
98 //________________________________________________________________________________________________________________________________________________
99 AliAnaPi0::~AliAnaPi0() {
100   // Remove event containers
101   if(fEventsList){
102     for(Int_t ic=0; ic<fNCentrBin; ic++){
103       for(Int_t iz=0; iz<fNZvertBin; iz++){
104         for(Int_t irp=0; irp<fNrpBin; irp++){
105           fEventsList[ic*fNZvertBin*fNrpBin+iz*fNrpBin+irp]->Delete() ;
106           delete fEventsList[ic*fNZvertBin*fNrpBin+iz*fNrpBin+irp] ;
107         }
108       }
109     }
110     delete[] fEventsList; 
111     fEventsList=0 ;
112   }
113   
114   if(fPHOSGeo)  delete fPHOSGeo  ;
115   if(fEMCALGeo) delete fEMCALGeo ;
116         
117 }
118
119 //________________________________________________________________________________________________________________________________________________
120 void AliAnaPi0::InitParameters()
121 {
122 //Init parameters when first called the analysis
123 //Set default parameters
124   SetInputAODName("PWG4Particle");
125   
126   AddToHistogramsName("AnaPi0_");
127
128   fNCentrBin = 1;
129   fNZvertBin = 1;
130   fNrpBin    = 1;
131   fNPID      = 9;
132   fNmaxMixEv = 10;
133   fZvtxCut   = 40;
134   fCalorimeter  = "PHOS";
135   fEMCALGeoName = "EMCAL_COMPLETE";
136 }
137 //________________________________________________________________________________________________________________________________________________
138 void AliAnaPi0::Init()
139 {  
140   //Init some data members needed in analysis
141   
142   //Histograms binning and range
143   if(!fhEtalon){                                                   //  p_T      alpha   d m_gg    
144     fhEtalon = new TH3D("hEtalon","Histo with binning parameters",50,0.,25.,10,0.,1.,200,0.,1.) ; 
145     fhEtalon->SetXTitle("P_{T} (GeV)") ;
146     fhEtalon->SetYTitle("#alpha") ;
147     fhEtalon->SetZTitle("m_{#gamma#gamma} (GeV)") ;
148   }
149   
150 }
151
152 //________________________________________________________________________________________________________________________________________________
153 TList * AliAnaPi0::GetCreateOutputObjects()
154 {  
155   // Create histograms to be saved in output file and 
156   // store them in fOutputContainer
157   
158   //create event containers
159   fEventsList = new TList*[fNCentrBin*fNZvertBin*fNrpBin] ;
160         
161   for(Int_t ic=0; ic<fNCentrBin; ic++){
162     for(Int_t iz=0; iz<fNZvertBin; iz++){
163       for(Int_t irp=0; irp<fNrpBin; irp++){
164         fEventsList[ic*fNZvertBin*fNrpBin+iz*fNrpBin+irp] = new TList() ;
165       }
166     }
167   }
168   
169   //If Geometry library loaded, do geometry selection during analysis.
170   printf("AliAnaPi0::GetCreateOutputObjects() - PHOS geometry initialized!\n");
171   fPHOSGeo = new AliPHOSGeoUtils("PHOSgeo") ;
172   
173   printf("AliAnaPi0::GetCreateOutputObjects() - EMCAL geometry initialized!\n");
174   fEMCALGeo = new AliEMCALGeoUtils(fEMCALGeoName) ;
175
176   TList * outputContainer = new TList() ; 
177   outputContainer->SetName(GetName()); 
178   
179   fhRe1 = new TH3D*[fNCentrBin*fNPID] ;
180   fhRe2 = new TH3D*[fNCentrBin*fNPID] ;
181   fhRe3 = new TH3D*[fNCentrBin*fNPID] ;
182   fhMi1 = new TH3D*[fNCentrBin*fNPID] ;
183   fhMi2 = new TH3D*[fNCentrBin*fNPID] ;
184   fhMi3 = new TH3D*[fNCentrBin*fNPID] ;
185   
186   char key[255] ;
187   char title[255] ;
188   
189   for(Int_t ic=0; ic<fNCentrBin; ic++){
190     for(Int_t ipid=0; ipid<fNPID; ipid++){
191       
192       //Distance to bad module 1
193       sprintf(key,"hRe_cen%d_pid%d_dist1",ic,ipid) ;
194       sprintf(title,"Real m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
195       
196       fhEtalon->Clone(key);
197       fhRe1[ic*fNPID+ipid]=(TH3D*)fhEtalon->Clone(key) ;
198       fhRe1[ic*fNPID+ipid]->SetName(key) ;
199       fhRe1[ic*fNPID+ipid]->SetTitle(title) ;
200       outputContainer->Add(fhRe1[ic*fNPID+ipid]) ;
201       
202       sprintf(key,"hMi_cen%d_pid%d_dist1",ic,ipid) ;
203       sprintf(title,"Mixed m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
204       fhMi1[ic*fNPID+ipid]=(TH3D*)fhEtalon->Clone(key) ;
205       fhMi1[ic*fNPID+ipid]->SetName(key) ;
206       fhMi1[ic*fNPID+ipid]->SetTitle(title) ;
207       outputContainer->Add(fhMi1[ic*fNPID+ipid]) ;
208       
209       //Distance to bad module 2
210       sprintf(key,"hRe_cen%d_pid%d_dist2",ic,ipid) ;
211       sprintf(title,"Real m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
212       fhRe2[ic*fNPID+ipid]=(TH3D*)fhEtalon->Clone(key) ;
213       fhRe2[ic*fNPID+ipid]->SetName(key) ;
214       fhRe2[ic*fNPID+ipid]->SetTitle(title) ;
215       outputContainer->Add(fhRe2[ic*fNPID+ipid]) ;
216       
217       sprintf(key,"hMi_cen%d_pid%d_dist2",ic,ipid) ;
218       sprintf(title,"Mixed m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
219       fhMi2[ic*fNPID+ipid]=(TH3D*)fhEtalon->Clone(key) ;
220       fhMi2[ic*fNPID+ipid]->SetName(key) ;
221       fhMi2[ic*fNPID+ipid]->SetTitle(title) ;
222       outputContainer->Add(fhMi2[ic*fNPID+ipid]) ;
223       
224       //Distance to bad module 3
225       sprintf(key,"hRe_cen%d_pid%d_dist3",ic,ipid) ;
226       sprintf(title,"Real m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
227       fhRe3[ic*fNPID+ipid]=(TH3D*)fhEtalon->Clone(key) ;
228       fhRe3[ic*fNPID+ipid]->SetName(key) ; 
229       fhRe3[ic*fNPID+ipid]->SetTitle(title) ;
230       outputContainer->Add(fhRe3[ic*fNPID+ipid]) ;
231       
232       sprintf(key,"hMi_cen%d_pid%d_dist3",ic,ipid) ;
233       sprintf(title,"Mixed m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
234       fhMi3[ic*fNPID+ipid]=(TH3D*)fhEtalon->Clone(key) ;
235       fhMi3[ic*fNPID+ipid]->SetName(key) ;
236       fhMi3[ic*fNPID+ipid]->SetTitle(title) ;
237       outputContainer->Add(fhMi3[ic*fNPID+ipid]) ;
238     }
239   }
240   
241   
242   fhEvents=new TH3D("hEvents","Number of events",fNCentrBin,0.,1.*fNCentrBin,
243                     fNZvertBin,0.,1.*fNZvertBin,fNrpBin,0.,1.*fNrpBin) ;
244   outputContainer->Add(fhEvents) ;
245   
246   //Histograms filled only if MC data is requested      
247   if(IsDataMC() || (GetReader()->GetDataType() == AliCaloTrackReader::kMC) ){
248     if(fhEtalon->GetXaxis()->GetXbins() && fhEtalon->GetXaxis()->GetXbins()->GetSize()){ //Variable bin size
249       fhPrimPt = new TH1D("hPrimPt","Primary pi0 pt",fhEtalon->GetXaxis()->GetNbins(),fhEtalon->GetXaxis()->GetXbins()->GetArray()) ;
250       fhPrimAccPt = new TH1D("hPrimAccPt","Primary pi0 pt with both photons in acceptance",fhEtalon->GetXaxis()->GetNbins(),
251                              fhEtalon->GetXaxis()->GetXbins()->GetArray()) ;
252     }
253     else{
254       fhPrimPt = new TH1D("hPrimPt","Primary pi0 pt",fhEtalon->GetXaxis()->GetNbins(),fhEtalon->GetXaxis()->GetXmin(),fhEtalon->GetXaxis()->GetXmax()) ;
255       fhPrimAccPt = new TH1D("hPrimAccPt","Primary pi0 pt with both photons in acceptance",
256                              fhEtalon->GetXaxis()->GetNbins(),fhEtalon->GetXaxis()->GetXmin(),fhEtalon->GetXaxis()->GetXmax()) ;
257     }
258     outputContainer->Add(fhPrimPt) ;
259     outputContainer->Add(fhPrimAccPt) ;
260     
261     fhPrimY = new TH1D("hPrimaryRapidity","Rapidity of primary pi0",100,-5.,5.) ; 
262     outputContainer->Add(fhPrimY) ;
263     
264     fhPrimAccY = new TH1D("hPrimAccRapidity","Rapidity of primary pi0",100,-5.,5.) ; 
265     outputContainer->Add(fhPrimAccY) ;
266     
267     fhPrimPhi = new TH1D("hPrimaryPhi","Azimithal of primary pi0",180,0.,360.) ; 
268     outputContainer->Add(fhPrimPhi) ;
269     
270     fhPrimAccPhi = new TH1D("hPrimAccPhi","Azimithal of primary pi0 with accepted daughters",180,-0.,360.) ; 
271     outputContainer->Add(fhPrimAccPhi) ;
272   }
273
274   //Save parameters used for analysis
275   TString parList ; //this will be list of parameters used for this analysis.
276   char onePar[255] ;
277   sprintf(onePar,"--- AliAnaPi0 ---\n") ;
278   parList+=onePar ;     
279   sprintf(onePar,"Number of bins in Centrality:  %d \n",fNCentrBin) ;
280   parList+=onePar ;
281   sprintf(onePar,"Number of bins in Z vert. pos: %d \n",fNZvertBin) ;
282   parList+=onePar ;
283   sprintf(onePar,"Number of bins in Reac. Plain: %d \n",fNrpBin) ;
284   parList+=onePar ;
285   sprintf(onePar,"Depth of event buffer: %d \n",fNmaxMixEv) ;
286   parList+=onePar ;
287   sprintf(onePar,"Number of different PID used:  %d \n",fNPID) ;
288   parList+=onePar ;
289   sprintf(onePar,"Cuts: \n") ;
290   parList+=onePar ;
291   sprintf(onePar,"Z vertex position: -%f < z < %f \n",fZvtxCut,fZvtxCut) ;
292   parList+=onePar ;
293   sprintf(onePar,"Calorimeter: %s \n",fCalorimeter.Data()) ;
294   parList+=onePar ;
295   
296   TObjString *oString= new TObjString(parList) ;
297   outputContainer->Add(oString);
298   
299   return outputContainer;
300 }
301
302 //_________________________________________________________________________________________________________________________________________________
303 void AliAnaPi0::Print(const Option_t * /*opt*/) const
304 {
305   //Print some relevant parameters set for the analysis
306   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
307   AliAnaPartCorrBaseClass::Print(" ");
308
309   printf("Number of bins in Centrality:  %d \n",fNCentrBin) ;
310   printf("Number of bins in Z vert. pos: %d \n",fNZvertBin) ;
311   printf("Number of bins in Reac. Plain: %d \n",fNrpBin) ;
312   printf("Depth of event buffer: %d \n",fNmaxMixEv) ;
313   printf("Number of different PID used:  %d \n",fNPID) ;
314   printf("Cuts: \n") ;
315   printf("Z vertex position: -%2.3f < z < %2.3f \n",fZvtxCut,fZvtxCut) ;
316   printf("------------------------------------------------------\n") ;
317
318
319
320 //____________________________________________________________________________________________________________________________________________________
321 void AliAnaPi0::MakeAnalysisFillHistograms() 
322 {
323   //Process one event and extract photons from AOD branch 
324   // filled with AliAnaPhoton and fill histos with invariant mass
325   
326   //Apply some cuts on event: vertex position and centrality range  
327   Int_t iRun=(GetReader()->GetInputEvent())->GetRunNumber() ;
328   if(IsBadRun(iRun)) return ;   
329   
330   Double_t vert[]={0,0,0} ; //vertex ;
331   GetReader()->GetVertex(vert);
332   if(vert[2]<-fZvtxCut || vert[2]> fZvtxCut) return ; //Event can not be used (vertex, centrality,... cuts not fulfilled)
333   
334   //Get Centrality and calculate centrality bin
335   //Does not exist in ESD yet???????
336   Int_t curCentrBin=0 ;
337   
338   //Get Reaction Plain position and calculate RP bin
339   //does not exist in ESD yet????
340   Int_t curRPBin=0 ;    
341   
342   Int_t curZvertBin=(Int_t)(0.5*fNZvertBin*(vert[2]+fZvtxCut)/fZvtxCut) ;
343   
344   fhEvents->Fill(curCentrBin+0.5,curZvertBin+0.5,curRPBin+0.5) ;
345   
346   Int_t nPhot = GetInputAODBranch()->GetEntriesFast() ;
347   if(GetDebug() > 1) printf("AliAnaPi0::MakeAnalysisFillHistograms() - Photon entries %d\n", nPhot);
348   
349   for(Int_t i1=0; i1<nPhot-1; i1++){
350     AliAODPWG4Particle * p1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i1)) ;
351     TLorentzVector photon1(p1->Px(),p1->Py(),p1->Pz(),p1->E());
352     for(Int_t i2=i1+1; i2<nPhot; i2++){
353       AliAODPWG4Particle * p2 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i2)) ;
354       TLorentzVector photon2(p2->Px(),p2->Py(),p2->Pz(),p2->E());
355       Double_t m  = (photon1 + photon2).M() ;
356       Double_t pt = (photon1 + photon2).Pt();
357       Double_t a  = TMath::Abs(p1->E()-p2->E())/(p1->E()+p2->E()) ;
358       if(GetDebug() > 2)
359         printf("AliAnaPi0::MakeAnalysisFillHistograms() - Current Event: pT: photon1 %2.2f, photon2 %2.2f; Pair: pT %2.2f, mass %2.3f, a %f2.3\n",
360                p1->Pt(), p2->Pt(), pt,m,a);
361       for(Int_t ipid=0; ipid<fNPID; ipid++)
362         {
363           if((p1->IsPIDOK(ipid,AliCaloPID::kPhoton)) && (p2->IsPIDOK(ipid,AliCaloPID::kPhoton))){ 
364             fhRe1[curCentrBin*fNPID+ipid]->Fill(pt,a,m) ;
365             if(p1->DistToBad()>0 && p2->DistToBad()>0){
366               fhRe2[curCentrBin*fNPID+ipid]->Fill(pt,a,m) ;
367               if(p1->DistToBad()>1 && p2->DistToBad()>1){
368                 fhRe3[curCentrBin*fNPID+ipid]->Fill(pt,a,m) ;
369               }
370             }
371           }
372         } 
373     }
374   }
375   
376   //Fill mixed
377   
378   TList * evMixList=fEventsList[curCentrBin*fNZvertBin*fNrpBin+curZvertBin*fNrpBin+curRPBin] ;
379   Int_t nMixed = evMixList->GetSize() ;
380   for(Int_t ii=0; ii<nMixed; ii++){  
381     TClonesArray* ev2= (TClonesArray*) (evMixList->At(ii));
382     Int_t nPhot2=ev2->GetEntriesFast() ;
383     Double_t m = -999;
384     if(GetDebug() > 1) printf("AliAnaPi0::MakeAnalysisFillHistograms() - Mixed event %d photon entries %d\n", ii, nPhot);
385     
386     for(Int_t i1=0; i1<nPhot; i1++){
387       AliAODPWG4Particle * p1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i1)) ;
388       TLorentzVector photon1(p1->Px(),p1->Py(),p1->Pz(),p1->E());
389       for(Int_t i2=0; i2<nPhot2; i2++){
390         AliAODPWG4Particle * p2 = (AliAODPWG4Particle*) (ev2->At(i2)) ;
391         
392         TLorentzVector photon2(p2->Px(),p2->Py(),p2->Pz(),p2->E());
393         m =           (photon1+photon2).M() ; 
394         Double_t pt = (photon1 + photon2).Pt();
395         Double_t a  = TMath::Abs(p1->E()-p2->E())/(p1->E()+p2->E()) ;
396         if(GetDebug() > 2)
397           printf("AliAnaPi0::MakeAnalysisFillHistograms() - Mixed Event: pT: photon1 %2.2f, photon2 %2.2f; Pair: pT %2.2f, mass %2.3f, a %f2.3\n",
398                  p1->Pt(), p2->Pt(), pt,m,a);                   
399         for(Int_t ipid=0; ipid<fNPID; ipid++){ 
400           if((p1->IsPIDOK(ipid,AliCaloPID::kPhoton)) && (p2->IsPIDOK(ipid,AliCaloPID::kPhoton))){ 
401             fhMi1[curCentrBin*fNPID+ipid]->Fill(pt,a,m) ;
402             if(p1->DistToBad()>0 && p2->DistToBad()>0){
403               fhMi2[curCentrBin*fNPID+ipid]->Fill(pt,a,m) ;
404               if(p1->DistToBad()>1 && p2->DistToBad()>1){
405                 fhMi3[curCentrBin*fNPID+ipid]->Fill(pt,a,m) ;
406               }
407               
408             }
409           }
410         }
411       }
412     }
413   }
414   
415   TClonesArray *currentEvent = new TClonesArray(*GetInputAODBranch());
416   //Add current event to buffer and Remove redandant events 
417   if(currentEvent->GetEntriesFast()>0){
418     evMixList->AddFirst(currentEvent) ;
419     currentEvent=0 ; //Now list of particles belongs to buffer and it will be deleted with buffer
420     if(evMixList->GetSize()>=fNmaxMixEv)
421       {
422         TClonesArray * tmp = (TClonesArray*) (evMixList->Last()) ;
423         evMixList->RemoveLast() ;
424         delete tmp ;
425       }
426   } 
427   else{ //empty event
428     delete currentEvent ;
429     currentEvent=0 ; 
430   }
431   
432   //Acceptance
433   if(GetReader()->ReadStack()){ 
434           AliStack * stack = GetMCStack();
435           if(stack && (IsDataMC() || (GetReader()->GetDataType() == AliCaloTrackReader::kMC)) ){
436              for(Int_t i=0 ; i<stack->GetNprimary(); i++){
437                  TParticle * prim = stack->Particle(i) ;
438                          if(prim->GetPdgCode() == 111){
439                      Double_t pi0Pt = prim->Pt() ;
440                      //printf("pi0, pt %2.2f\n",pi0Pt);
441                     if(prim->Energy() == TMath::Abs(prim->Pz()))  continue ; //Protection against floating point exception        
442                     Double_t pi0Y  = 0.5*TMath::Log((prim->Energy()-prim->Pz())/(prim->Energy()+prim->Pz())) ;
443                     Double_t phi   = TMath::RadToDeg()*prim->Phi() ;
444                 if(TMath::Abs(pi0Y) < 0.5){
445                    fhPrimPt->Fill(pi0Pt) ;
446                 }
447                 fhPrimY  ->Fill(pi0Y) ;
448                     fhPrimPhi->Fill(phi) ;
449         
450                 //Check if both photons hit Calorimeter
451                 Int_t iphot1=prim->GetFirstDaughter() ;
452                 Int_t iphot2=prim->GetLastDaughter() ;
453                 if(iphot1>-1 && iphot1<stack->GetNtrack() && iphot2>-1 && iphot2<stack->GetNtrack()){
454                    TParticle * phot1 = stack->Particle(iphot1) ;
455                    TParticle * phot2 = stack->Particle(iphot2) ;
456                    if(phot1 && phot2 && phot1->GetPdgCode()==22 && phot2->GetPdgCode()==22){
457                    //printf("2 photons: photon 1: pt %2.2f, phi %3.2f, eta %1.2f; photon 2: pt %2.2f, phi %3.2f, eta %1.2f\n",
458                            //   phot1->Pt(), phot1->Phi()*180./3.1415, phot1->Eta(), phot2->Pt(), phot2->Phi()*180./3.1415, phot2->Eta());
459                    Bool_t inacceptance = kFALSE;
460
461                            if(fCalorimeter == "PHOS"){
462                                    if(fPHOSGeo){
463                                            Int_t mod ;
464                                            Double_t x,z ;
465                                            if(fPHOSGeo->ImpactOnEmc(phot1,mod,z,x) && fPHOSGeo->ImpactOnEmc(phot2,mod,z,x)) 
466                                                    inacceptance = kTRUE;
467                                            if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
468                                    }
469                                    else{
470                                            TLorentzVector lv1, lv2;
471                                            phot1->Momentum(lv1);
472                                            phot2->Momentum(lv2);
473                                            if(GetFidutialCut()->IsInFidutialCut(lv1,fCalorimeter) && GetFidutialCut()->IsInFidutialCut(lv2,fCalorimeter)) 
474                                            inacceptance = kTRUE ;
475                                            if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
476                                    }
477                                    
478                            }       
479                            else if(fCalorimeter == "EMCAL"){
480                                    if(fEMCALGeo){
481                                            if(fEMCALGeo->Impact(phot1) && fEMCALGeo->Impact(phot2)) 
482                                                    inacceptance = kTRUE;
483                                            if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
484                                         }
485                                    else{
486                                            TLorentzVector lv1, lv2;
487                                            phot1->Momentum(lv1);
488                                            phot2->Momentum(lv2);
489                                            if(GetFidutialCut()->IsInFidutialCut(lv1,fCalorimeter) && GetFidutialCut()->IsInFidutialCut(lv2,fCalorimeter)) 
490                                                    inacceptance = kTRUE ;
491                                            if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
492                                    }
493                    }      
494                 
495                    if(inacceptance){
496                       fhPrimAccPt->Fill(pi0Pt) ;
497                       fhPrimAccPhi->Fill(phi) ;
498                       fhPrimAccY->Fill(pi0Y) ;
499                    }//Accepted
500                }// 2 photons      
501              }//Check daughters exist
502        }// Primary pi0
503      }//loop on primaries       
504    }//stack exists and data is MC
505   }//read stack
506   else if(GetReader()->ReadAODMCParticles()){
507                 if(GetDebug() >= 0)  printf("AliAnaPi0::MakeAnalysisFillHistograms() - Acceptance calculation with MCParticles not implemented yet\n");
508   }     
509   
510 }       
511
512 //________________________________________________________________________
513 void AliAnaPi0::ReadHistograms(TList* outputList)
514 {
515         // Needed when Terminate is executed in distributed environment
516         // Refill analysis histograms of this class with corresponding histograms in output list. 
517         
518         // Histograms of this analsys are kept in the same list as other analysis, recover the position of
519         // the first one and then add the next.
520         Int_t index = outputList->IndexOf(outputList->FindObject(GetAddedHistogramsStringToName()+"hRe_cen0_pid0_dist1"));
521         
522         if(!fhRe1) fhRe1 = new TH3D*[fNCentrBin*fNPID] ;
523         if(!fhRe2) fhRe2 = new TH3D*[fNCentrBin*fNPID] ;
524         if(!fhRe3) fhRe3 = new TH3D*[fNCentrBin*fNPID] ;
525         if(!fhMi1) fhMi1 = new TH3D*[fNCentrBin*fNPID] ;
526         if(!fhMi2) fhMi2 = new TH3D*[fNCentrBin*fNPID] ;
527         if(!fhMi3) fhMi3 = new TH3D*[fNCentrBin*fNPID] ;        
528         
529     for(Int_t ic=0; ic<fNCentrBin; ic++){
530         for(Int_t ipid=0; ipid<fNPID; ipid++){
531                         fhRe1[ic*fNPID+ipid] = (TH3D*) outputList->At(index++);
532                         fhMi1[ic*fNPID+ipid] = (TH3D*) outputList->At(index++);
533                         fhRe2[ic*fNPID+ipid] = (TH3D*) outputList->At(index++);
534                         fhMi2[ic*fNPID+ipid] = (TH3D*) outputList->At(index++);
535                         fhRe3[ic*fNPID+ipid] = (TH3D*) outputList->At(index++);
536                         fhMi3[ic*fNPID+ipid] = (TH3D*) outputList->At(index++);
537                 }
538         }
539         
540         fhEvents = (TH3D *) outputList->At(index++); 
541         
542         //Histograms filled only if MC data is requested        
543         if(IsDataMC() || (GetReader()->GetDataType() == AliCaloTrackReader::kMC) ){
544                  fhPrimPt     = (TH1D*)  outputList->At(index++);
545                  fhPrimAccPt  = (TH1D*)  outputList->At(index++);
546                  fhPrimY      = (TH1D*)  outputList->At(index++);
547                  fhPrimAccY   = (TH1D*)  outputList->At(index++);
548                  fhPrimPhi    = (TH1D*)  outputList->At(index++);
549                  fhPrimAccPhi = (TH1D*)  outputList->At(index);
550         }
551         
552 }
553
554
555 //____________________________________________________________________________________________________________________________________________________
556 void AliAnaPi0::Terminate(TList* outputList) 
557 {
558   //Do some calculations and plots from the final histograms.
559   
560   printf(" *** %s Terminate:\n", GetName()) ; 
561  
562   //Recover histograms from output histograms list, needed for distributed analysis.    
563   ReadHistograms(outputList);
564         
565   if (!fhRe1) {
566      printf("AliAnaPi0::Terminate() - Error: Remote output histograms not imported in AliAnaPi0 object");
567      return;
568   }
569         
570   printf("AliAnaPi0::Terminate()         Mgg Real        : %5.3f , RMS : %5.3f \n", fhRe1[0]->GetMean(),   fhRe1[0]->GetRMS() ) ;
571  
572   char nameIM[128];
573   sprintf(nameIM,"AliAnaPi0_%s_cPt",fCalorimeter.Data());
574   TCanvas  * cIM = new TCanvas(nameIM, "", 400, 10, 600, 700) ;
575   cIM->Divide(2, 2);
576
577   cIM->cd(1) ; 
578   //gPad->SetLogy();
579   TH1D * hIMAllPt = (TH1D*) fhRe1[0]->ProjectionZ();
580   hIMAllPt->SetLineColor(2);
581   hIMAllPt->SetTitle("No cut on  p_{T, #gamma#gamma} ");
582   hIMAllPt->Draw();
583
584   cIM->cd(2) ; 
585   TH3F * hRe1Pt5 = (TH3F*)fhRe1[0]->Clone("IMPt5");
586   hRe1Pt5->GetXaxis()->SetRangeUser(0,5);
587   TH1D * hIMPt5 = (TH1D*) hRe1Pt5->Project3D("z");
588   hIMPt5->SetLineColor(2);  
589   hIMPt5->SetTitle("0 < p_{T, #gamma#gamma} < 5 GeV/c");
590   hIMPt5->Draw();
591   
592   cIM->cd(3) ; 
593   TH3F * hRe1Pt10 =  (TH3F*)fhRe1[0]->Clone("IMPt10");
594   hRe1Pt10->GetXaxis()->SetRangeUser(5,10);
595   TH1D * hIMPt10 = (TH1D*) hRe1Pt10->Project3D("z");
596   hIMPt10->SetLineColor(2);  
597   hIMPt10->SetTitle("5 < p_{T, #gamma#gamma} < 10 GeV/c");
598   hIMPt10->Draw();
599   
600   cIM->cd(4) ; 
601   TH3F * hRe1Pt20 =  (TH3F*)fhRe1[0]->Clone("IMPt20");
602   hRe1Pt20->GetXaxis()->SetRangeUser(10,20);
603   TH1D * hIMPt20 = (TH1D*) hRe1Pt20->Project3D("z");
604   hIMPt20->SetLineColor(2);  
605   hIMPt20->SetTitle("10 < p_{T, #gamma#gamma} < 20 GeV/c");
606   hIMPt20->Draw();
607    
608   char nameIMF[128];
609   sprintf(nameIMF,"AliAnaPi0_%s_Mgg.eps",fCalorimeter.Data());
610   cIM->Print(nameIMF);
611
612   char namePt[128];
613   sprintf(namePt,"AliAnaPi0_%s_cPt",fCalorimeter.Data());
614   TCanvas  * cPt = new TCanvas(namePt, "", 400, 10, 600, 700) ;
615   cPt->Divide(2, 2);
616
617   cPt->cd(1) ; 
618   //gPad->SetLogy();
619   TH1D * hPt = (TH1D*) fhRe1[0]->Project3D("x");
620   hPt->SetLineColor(2);
621   hPt->SetTitle("No cut on  M_{#gamma#gamma} ");
622   hPt->Draw();
623
624   cPt->cd(2) ; 
625   TH3F * hRe1IM1 = (TH3F*)fhRe1[0]->Clone("PtIM1");
626   hRe1IM1->GetZaxis()->SetRangeUser(0.05,0.21);
627   TH1D * hPtIM1 = (TH1D*) hRe1IM1->Project3D("x");
628   hPtIM1->SetLineColor(2);  
629   hPtIM1->SetTitle("0.05 < M_{#gamma#gamma} < 0.21 GeV/c^{2}");
630   hPtIM1->Draw();
631   
632   cPt->cd(3) ; 
633   TH3F * hRe1IM2 = (TH3F*)fhRe1[0]->Clone("PtIM2");
634   hRe1IM2->GetZaxis()->SetRangeUser(0.09,0.17);
635   TH1D * hPtIM2 = (TH1D*) hRe1IM2->Project3D("x");
636   hPtIM2->SetLineColor(2);  
637   hPtIM2->SetTitle("0.09 < M_{#gamma#gamma} < 0.17 GeV/c^{2}");
638   hPtIM2->Draw();
639
640   cPt->cd(4) ; 
641   TH3F * hRe1IM3 = (TH3F*)fhRe1[0]->Clone("PtIM3");
642   hRe1IM3->GetZaxis()->SetRangeUser(0.11,0.15);
643   TH1D * hPtIM3 = (TH1D*) hRe1IM1->Project3D("x");
644   hPtIM3->SetLineColor(2);  
645   hPtIM3->SetTitle("0.11 < M_{#gamma#gamma} < 0.15 GeV/c^{2}");
646   hPtIM3->Draw();
647    
648   char namePtF[128];
649   sprintf(namePtF,"AliAnaPi0_%s_Pt.eps",fCalorimeter.Data());
650   cPt->Print(namePtF);
651
652  
653   char line[1024] ; 
654   sprintf(line, ".!tar -zcf %s_%s.tar.gz *.eps", GetName(),fCalorimeter.Data()) ; 
655   gROOT->ProcessLine(line);
656   sprintf(line, ".!rm -fR AliAnaPi0_%s*.eps",fCalorimeter.Data()); 
657   gROOT->ProcessLine(line);
658  
659   printf(" AliAnaPi0::Terminate() - !! All the eps files are in %s_%s.tar.gz !!!\n", GetName(), fCalorimeter.Data());
660
661 }
662
663
664
665
666