Modifications needed to make the frame work in different working modes bu Andrei...
[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
34 //---- AliRoot system ----
35 #include "AliAnaPi0.h"
36 #include "AliCaloTrackReader.h"
37 #include "AliCaloPID.h"
38 #include "AliStack.h"
39
40 #ifdef __PHOSGEO__
41         #include "AliPHOSGeoUtils.h"
42 #endif
43
44 ClassImp(AliAnaPi0)
45
46 //________________________________________________________________________________________________________________________________________________  
47 AliAnaPi0::AliAnaPi0() : AliAnaPartCorrBaseClass(),
48 fNCentrBin(0),fNZvertBin(0),fNrpBin(0),
49 fNPID(0),fNmaxMixEv(0), fZvtxCut(0.),fCalorimeter(""),
50 fEventsList(0x0), fhEtalon(0x0),
51 fhRe1(0x0),fhMi1(0x0),fhRe2(0x0),fhMi2(0x0),fhRe3(0x0),fhMi3(0x0),fhEvents(0x0),
52 fhPrimPt(0x0), fhPrimAccPt(0x0), fhPrimY(0x0), fhPrimAccY(0x0), fhPrimPhi(0x0), fhPrimAccPhi(0x0)
53 {
54 //Default Ctor
55  InitParameters();
56  
57 }
58
59 //________________________________________________________________________________________________________________________________________________
60 AliAnaPi0::AliAnaPi0(const AliAnaPi0 & ex) : AliAnaPartCorrBaseClass(ex),  
61 fNCentrBin(ex.fNCentrBin),fNZvertBin(ex.fNZvertBin),fNrpBin(ex.fNrpBin),
62 fNPID(ex.fNPID),fNmaxMixEv(ex.fNmaxMixEv),fZvtxCut(ex.fZvtxCut), fCalorimeter(ex.fCalorimeter), 
63 fEventsList(ex.fEventsList), fhEtalon(ex.fhEtalon),
64 fhRe1(ex.fhRe1),fhMi1(ex.fhMi1),fhRe2(ex.fhRe2),fhMi2(ex.fhMi2),fhRe3(ex.fhRe3),fhMi3(ex.fhMi3),fhEvents(ex.fhEvents),
65 fhPrimPt(ex.fhPrimPt), fhPrimAccPt(ex.fhPrimAccPt), fhPrimY(ex.fhPrimY), 
66 fhPrimAccY(ex.fhPrimAccY), fhPrimPhi(ex.fhPrimPhi), fhPrimAccPhi(ex.fhPrimAccPhi)
67 {
68         // cpy ctor
69         //Do not need it
70 }
71
72 //________________________________________________________________________________________________________________________________________________
73 AliAnaPi0 & AliAnaPi0::operator = (const AliAnaPi0 & ex)
74 {
75         // assignment operator
76         
77         if(this == &ex)return *this;
78         ((AliAnaPartCorrBaseClass *)this)->operator=(ex);
79         
80         fNCentrBin = ex.fNCentrBin  ; fNZvertBin = ex.fNZvertBin  ; fNrpBin = ex.fNrpBin  ; 
81         fNPID = ex.fNPID  ; fNmaxMixEv = ex.fNmaxMixEv  ; fZvtxCut = ex.fZvtxCut  ;  fCalorimeter = ex.fCalorimeter  ;  
82         fEventsList = ex.fEventsList  ;  fhEtalon = ex.fhEtalon  ; 
83         fhRe1 = ex.fhRe1  ; fhMi1 = ex.fhMi1  ; fhRe2 = ex.fhRe2  ; fhMi2 = ex.fhMi2  ; 
84         fhRe3 = ex.fhRe3  ; fhMi3 = ex.fhMi3  ; fhEvents = ex.fhEvents  ; 
85         fhPrimPt = ex.fhPrimPt  ;  fhPrimAccPt = ex.fhPrimAccPt  ;  fhPrimY = ex.fhPrimY  ;  
86         fhPrimAccY = ex.fhPrimAccY  ;  fhPrimPhi = ex.fhPrimPhi  ;  fhPrimAccPhi = ex.fhPrimAccPhi ;
87         
88         return *this;
89         
90 }
91
92 //________________________________________________________________________________________________________________________________________________
93 AliAnaPi0::~AliAnaPi0() {
94         // Remove event containers
95         if(fEventsList){
96                 for(Int_t ic=0; ic<fNCentrBin; ic++){
97                         for(Int_t iz=0; iz<fNZvertBin; iz++){
98                                 for(Int_t irp=0; irp<fNrpBin; irp++){
99                                         fEventsList[ic*fNZvertBin*fNrpBin+iz*fNrpBin+irp]->Delete() ;
100                                         delete fEventsList[ic*fNZvertBin*fNrpBin+iz*fNrpBin+irp] ;
101                                 }
102                         }
103                 }
104                 delete[] fEventsList; 
105                 fEventsList=0 ;
106         }
107         
108 #ifdef __PHOSGEO__
109     if(fPHOSGeo) delete fPHOSGeo ;
110 #endif  
111 }
112
113 //________________________________________________________________________________________________________________________________________________
114 void AliAnaPi0::InitParameters()
115 {
116 //Init parameters when first called the analysis
117 //Set default parameters
118         SetInputAODName("photons");
119         fNCentrBin = 1;
120         fNZvertBin = 1;
121         fNrpBin    = 1;
122     fNPID      = 9;
123         fNmaxMixEv = 10;
124         fZvtxCut   = 40;
125         fCalorimeter  = "PHOS";
126 }
127 //________________________________________________________________________________________________________________________________________________
128 void AliAnaPi0::Init()
129 {  
130 //Init some data members needed in analysis
131
132 //Histograms binning and range
133         if(!fhEtalon){                                                   //  p_T      alpha   d m_gg    
134                 fhEtalon = new TH3D("hEtalon","Histo with binning parameters",50,0.,25.,10,0.,1.,200,0.,1.) ; 
135                 fhEtalon->SetXTitle("P_{T} (GeV)") ;
136                 fhEtalon->SetYTitle("#alpha") ;
137                 fhEtalon->SetZTitle("m_{#gamma#gamma} (GeV)") ;
138         }
139         
140 }
141
142 //________________________________________________________________________________________________________________________________________________
143 TList * AliAnaPi0::GetCreateOutputObjects()
144 {  
145         // Create histograms to be saved in output file and 
146         // store them in fOutputContainer
147         
148         //create event containers
149         fEventsList = new TList*[fNCentrBin*fNZvertBin*fNrpBin] ;
150         
151         for(Int_t ic=0; ic<fNCentrBin; ic++){
152                 for(Int_t iz=0; iz<fNZvertBin; iz++){
153                         for(Int_t irp=0; irp<fNrpBin; irp++){
154                                 fEventsList[ic*fNZvertBin*fNrpBin+iz*fNrpBin+irp] = new TList() ;
155                         }
156                 }
157         }
158         
159                 //If Geometry library loaded, do geometry selection during analysis.
160 #ifdef __PHOSGEO__
161         printf("PHOS geometry initialized!\n");
162         fPHOSGeo = new AliPHOSGeoUtils("PHOSgeo") ;
163 #endif  
164         
165         TList * outputContainer = new TList() ; 
166         outputContainer->SetName(GetName()); 
167         
168         fhRe1=new TH3D*[fNCentrBin*fNPID] ;
169         fhRe2=new TH3D*[fNCentrBin*fNPID] ;
170         fhRe3=new TH3D*[fNCentrBin*fNPID] ;
171         fhMi1=new TH3D*[fNCentrBin*fNPID] ;
172         fhMi2=new TH3D*[fNCentrBin*fNPID] ;
173         fhMi3=new TH3D*[fNCentrBin*fNPID] ;
174         
175         char key[255] ;
176         char title[255] ;
177         
178         for(Int_t ic=0; ic<fNCentrBin; ic++){
179                 for(Int_t ipid=0; ipid<fNPID; ipid++){
180                                         
181                         //Distance to bad module 1
182                         sprintf(key,"hRe_cen%d_pid%d_dist1",ic,ipid) ;
183                         sprintf(title,"Real m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
184                         
185                         fhEtalon->Clone(key);
186                         fhRe1[ic*fNPID+ipid]=(TH3D*)fhEtalon->Clone(key) ;
187                         fhRe1[ic*fNPID+ipid]->SetName(key) ;
188                         fhRe1[ic*fNPID+ipid]->SetTitle(title) ;
189                         outputContainer->Add(fhRe1[ic*fNPID+ipid]) ;
190                         
191                         sprintf(key,"hMi_cen%d_pid%d_dist1",ic,ipid) ;
192                         sprintf(title,"Mixed m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
193                         fhMi1[ic*fNPID+ipid]=(TH3D*)fhEtalon->Clone(key) ;
194                         fhMi1[ic*fNPID+ipid]->SetName(key) ;
195                         fhMi1[ic*fNPID+ipid]->SetTitle(title) ;
196                         outputContainer->Add(fhMi1[ic*fNPID+ipid]) ;
197                         
198                         //Distance to bad module 2
199                         sprintf(key,"hRe_cen%d_pid%d_dist2",ic,ipid) ;
200                         sprintf(title,"Real m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
201                         fhRe2[ic*fNPID+ipid]=(TH3D*)fhEtalon->Clone(key) ;
202                         fhRe2[ic*fNPID+ipid]->SetName(key) ;
203                         fhRe2[ic*fNPID+ipid]->SetTitle(title) ;
204                         outputContainer->Add(fhRe2[ic*fNPID+ipid]) ;
205                         
206                         sprintf(key,"hMi_cen%d_pid%d_dist2",ic,ipid) ;
207                         sprintf(title,"Mixed m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
208                         fhMi2[ic*fNPID+ipid]=(TH3D*)fhEtalon->Clone(key) ;
209                         fhMi2[ic*fNPID+ipid]->SetName(key) ;
210                         fhMi2[ic*fNPID+ipid]->SetTitle(title) ;
211                         outputContainer->Add(fhMi2[ic*fNPID+ipid]) ;
212                         
213                         //Distance to bad module 3
214                         sprintf(key,"hRe_cen%d_pid%d_dist3",ic,ipid) ;
215                         sprintf(title,"Real m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
216                         fhRe3[ic*fNPID+ipid]=(TH3D*)fhEtalon->Clone(key) ;
217                         fhRe3[ic*fNPID+ipid]->SetName(key) ; 
218                         fhRe3[ic*fNPID+ipid]->SetTitle(title) ;
219                         outputContainer->Add(fhRe3[ic*fNPID+ipid]) ;
220                         
221                         sprintf(key,"hMi_cen%d_pid%d_dist3",ic,ipid) ;
222                         sprintf(title,"Mixed m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
223                         fhMi3[ic*fNPID+ipid]=(TH3D*)fhEtalon->Clone(key) ;
224                         fhMi3[ic*fNPID+ipid]->SetName(key) ;
225                         fhMi3[ic*fNPID+ipid]->SetTitle(title) ;
226                         outputContainer->Add(fhMi3[ic*fNPID+ipid]) ;
227                 }
228     }
229         
230         
231         fhEvents=new TH3D("hEvents","Number of events",fNCentrBin,0.,1.*fNCentrBin,
232                                           fNZvertBin,0.,1.*fNZvertBin,fNrpBin,0.,1.*fNrpBin) ;
233         outputContainer->Add(fhEvents) ;
234
235         //Histograms filled only if MC data is requested        
236         if(IsDataMC() || (GetReader()->GetDataType() == AliCaloTrackReader::kMC) ){
237                 if(fhEtalon->GetXaxis()->GetXbins() && fhEtalon->GetXaxis()->GetXbins()->GetSize()){ //Variable bin size
238                         fhPrimPt = new TH1D("hPrimPt","Primary pi0 pt",fhEtalon->GetXaxis()->GetNbins(),fhEtalon->GetXaxis()->GetXbins()->GetArray()) ;
239                         fhPrimAccPt = new TH1D("hPrimAccPt","Primary pi0 pt with both photons in acceptance",fhEtalon->GetXaxis()->GetNbins(),
240                                                                    fhEtalon->GetXaxis()->GetXbins()->GetArray()) ;
241                 }
242                 else{
243                         fhPrimPt = new TH1D("hPrimPt","Primary pi0 pt",fhEtalon->GetXaxis()->GetNbins(),fhEtalon->GetXaxis()->GetXmin(),fhEtalon->GetXaxis()->GetXmax()) ;
244                         fhPrimAccPt = new TH1D("hPrimAccPt","Primary pi0 pt with both photons in acceptance",
245                                                                    fhEtalon->GetXaxis()->GetNbins(),fhEtalon->GetXaxis()->GetXmin(),fhEtalon->GetXaxis()->GetXmax()) ;
246                 }
247                 outputContainer->Add(fhPrimPt) ;
248                 outputContainer->Add(fhPrimAccPt) ;
249                 
250                 fhPrimY = new TH1D("hPrimaryRapidity","Rapidity of primary pi0",100,-5.,5.) ; 
251                 outputContainer->Add(fhPrimY) ;
252                 
253                 fhPrimAccY = new TH1D("hPrimAccRapidity","Rapidity of primary pi0",100,-5.,5.) ; 
254                 outputContainer->Add(fhPrimAccY) ;
255                 
256                 fhPrimPhi = new TH1D("hPrimaryPhi","Azimithal of primary pi0",180,0.,360.) ; 
257                 outputContainer->Add(fhPrimPhi) ;
258                 
259                 fhPrimAccPhi = new TH1D("hPrimAccPhi","Azimithal of primary pi0 with accepted daughters",180,-0.,360.) ; 
260                 outputContainer->Add(fhPrimAccPhi) ;
261         }
262         
263         //Save parameters used for analysis
264         TString parList ; //this will be list of parameters used for this analysis.
265         char onePar[255] ;
266         sprintf(onePar,"--- AliAnaPi0 ---\n") ;
267         parList+=onePar ;       
268         sprintf(onePar,"Number of bins in Centrality:  %d \n",fNCentrBin) ;
269         parList+=onePar ;
270         sprintf(onePar,"Number of bins in Z vert. pos: %d \n",fNZvertBin) ;
271         parList+=onePar ;
272         sprintf(onePar,"Number of bins in Reac. Plain: %d \n",fNrpBin) ;
273         parList+=onePar ;
274         sprintf(onePar,"Depth of event buffer: %d \n",fNmaxMixEv) ;
275         parList+=onePar ;
276         sprintf(onePar,"Number of different PID used:  %d \n",fNPID) ;
277         parList+=onePar ;
278         sprintf(onePar,"Cuts: \n") ;
279         parList+=onePar ;
280         sprintf(onePar,"Z vertex position: -%f < z < %f \n",fZvtxCut,fZvtxCut) ;
281         parList+=onePar ;
282         sprintf(onePar,"Calorimeter: %s \n",fCalorimeter.Data()) ;
283         parList+=onePar ;
284         
285         TObjString *oString= new TObjString(parList) ;
286         outputContainer->Add(oString);
287         
288         return outputContainer;
289 }
290
291 //_________________________________________________________________________________________________________________________________________________
292 void AliAnaPi0::Print(const Option_t * /*opt*/) const
293 {
294         //Print some relevant parameters set for the analysis
295         AliAnaPartCorrBaseClass::Print(" ");
296         printf("Class AliAnaPi0 for gamma-gamma inv.mass construction \n") ;
297         printf("Number of bins in Centrality:  %d \n",fNCentrBin) ;
298         printf("Number of bins in Z vert. pos: %d \n",fNZvertBin) ;
299         printf("Number of bins in Reac. Plain: %d \n",fNrpBin) ;
300         printf("Depth of event buffer: %d \n",fNmaxMixEv) ;
301         printf("Number of different PID used:  %d \n",fNPID) ;
302         printf("Cuts: \n") ;
303         printf("Z vertex position: -%2.3f < z < %2.3f \n",fZvtxCut,fZvtxCut) ;
304         printf("------------------------------------------------------\n") ;
305
306
307
308 //____________________________________________________________________________________________________________________________________________________
309 void AliAnaPi0::MakeAnalysisFillHistograms() 
310 {
311         //Process one event and extract photons from AOD branch 
312         // filled with AliAnaPhoton and fill histos with invariant mass
313         
314         //Apply some cuts on event: vertex position and centrality range  
315         Int_t iRun=(GetReader()->GetInputEvent())->GetRunNumber() ;
316         if(IsBadRun(iRun)) return ;     
317         
318         Double_t vert[]={0,0,0} ; //vertex ;
319         GetReader()->GetVertex(vert);
320         if(vert[2]<-fZvtxCut || vert[2]> fZvtxCut) return ; //Event can not be used (vertex, centrality,... cuts not fulfilled)
321
322         //Get Centrality and calculate centrality bin
323         //Does not exist in ESD yet???????
324         Int_t curCentrBin=0 ;
325  
326         //Get Reaction Plain position and calculate RP bin
327         //does not exist in ESD yet????
328         Int_t curRPBin=0 ;      
329         
330         Int_t curZvertBin=(Int_t)(0.5*fNZvertBin*(vert[2]+fZvtxCut)/fZvtxCut) ;
331                 
332         fhEvents->Fill(curCentrBin+0.5,curZvertBin+0.5,curRPBin+0.5) ;
333                         
334         Int_t nPhot = GetInputAODBranch()->GetEntriesFast() ;
335         if(GetDebug() > 1) printf("AliAnaPi0::FillHistos: photon entries %d\n", nPhot);
336         
337         for(Int_t i1=0; i1<nPhot-1; i1++){
338                 AliAODPWG4Particle * p1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i1)) ;
339                 TLorentzVector photon1(p1->Px(),p1->Py(),p1->Pz(),p1->E());
340                 for(Int_t i2=i1+1; i2<nPhot; i2++){
341                         AliAODPWG4Particle * p2 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i2)) ;
342                         TLorentzVector photon2(p2->Px(),p2->Py(),p2->Pz(),p2->E());
343                         Double_t m  = (photon1 + photon2).M() ;
344                         Double_t pt = (photon1 + photon2).Pt();
345                         Double_t a  = TMath::Abs(p1->E()-p2->E())/(p1->E()+p2->E()) ;
346                         if(GetDebug() > 2)
347                                 printf("AliAnaPi0::FillHistos: Current Event: pT: photon1 %2.2f, photon2 %2.2f; Pair: pT %2.2f, mass %2.3f, a %f2.3\n",
348                                                 p1->Pt(), p2->Pt(), pt,m,a);
349                         for(Int_t ipid=0; ipid<fNPID; ipid++)
350                         {
351                                 if((p1->IsPIDOK(ipid,AliCaloPID::kPhoton)) && (p2->IsPIDOK(ipid,AliCaloPID::kPhoton))){ 
352                                         fhRe1[curCentrBin*fNPID+ipid]->Fill(pt,a,m) ;
353                                         if(p1->DistToBad()>0 && p2->DistToBad()>0){
354                                                 fhRe2[curCentrBin*fNPID+ipid]->Fill(pt,a,m) ;
355                                                 if(p1->DistToBad()>1 && p2->DistToBad()>1){
356                                                         fhRe3[curCentrBin*fNPID+ipid]->Fill(pt,a,m) ;
357                                                 }
358                                         }
359                                 }
360                         } 
361                 }
362         }
363                 
364         //Fill mixed
365         
366         TList * evMixList=fEventsList[curCentrBin*fNZvertBin*fNrpBin+curZvertBin*fNrpBin+curRPBin] ;
367         Int_t nMixed = evMixList->GetSize() ;
368         for(Int_t ii=0; ii<nMixed; ii++){  
369                 TClonesArray* ev2= (TClonesArray*) (evMixList->At(ii));
370                 Int_t nPhot2=ev2->GetEntriesFast() ;
371                 Double_t m = -999;
372                 if(GetDebug() > 1) printf("AliAnaPi0::FillHistos: Mixed event %d photon entries %d\n", ii, nPhot);
373                 
374                 for(Int_t i1=0; i1<nPhot; i1++){
375                         AliAODPWG4Particle * p1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i1)) ;
376                         TLorentzVector photon1(p1->Px(),p1->Py(),p1->Pz(),p1->E());
377                         for(Int_t i2=0; i2<nPhot2; i2++){
378                                 AliAODPWG4Particle * p2 = (AliAODPWG4Particle*) (ev2->At(i2)) ;
379                                 
380                                 TLorentzVector photon2(p2->Px(),p2->Py(),p2->Pz(),p2->E());
381                                 m =           (photon1+photon2).M() ; 
382                                 Double_t pt = (photon1 + photon2).Pt();
383                                 Double_t a  = TMath::Abs(p1->E()-p2->E())/(p1->E()+p2->E()) ;
384                                 if(GetDebug() > 2)
385                                         printf("AliAnaPi0::FillHistos: Mixed Event: pT: photon1 %2.2f, photon2 %2.2f; Pair: pT %2.2f, mass %2.3f, a %f2.3\n",
386                                                         p1->Pt(), p2->Pt(), pt,m,a);                    
387                                 for(Int_t ipid=0; ipid<fNPID; ipid++){ 
388                                 if((p1->IsPIDOK(ipid,AliCaloPID::kPhoton)) && (p2->IsPIDOK(ipid,AliCaloPID::kPhoton))){ 
389                                                 fhMi1[curCentrBin*fNPID+ipid]->Fill(pt,a,m) ;
390                                                 if(p1->DistToBad()>0 && p2->DistToBad()>0){
391                                                         fhMi2[curCentrBin*fNPID+ipid]->Fill(pt,a,m) ;
392                                                         if(p1->DistToBad()>1 && p2->DistToBad()>1){
393                                                                 fhMi3[curCentrBin*fNPID+ipid]->Fill(pt,a,m) ;
394                                                         }
395                                                         
396                                                 }
397                                         }
398                                 }
399                         }
400                 }
401         }
402         
403         TClonesArray *currentEvent = new TClonesArray(*GetInputAODBranch());
404         //Add current event to buffer and Remove redandant events 
405         if(currentEvent->GetEntriesFast()>0){
406                 evMixList->AddFirst(currentEvent) ;
407                 currentEvent=0 ; //Now list of particles belongs to buffer and it will be deleted with buffer
408                 if(evMixList->GetSize()>=fNmaxMixEv)
409                 {
410                         TClonesArray * tmp = (TClonesArray*) (evMixList->Last()) ;
411                         evMixList->RemoveLast() ;
412                         delete tmp ;
413                 }
414         } 
415         else{ //empty event
416                 delete currentEvent ;
417                 currentEvent=0 ; 
418         }
419                 
420         //Acceptance
421         AliStack * stack = GetMCStack();
422         if(stack && (IsDataMC() || (GetReader()->GetDataType() == AliCaloTrackReader::kMC)) ){
423                 for(Int_t i=0 ; i<stack->GetNprimary(); i++){
424                         TParticle * prim = stack->Particle(i) ;
425                         if(prim->GetPdgCode() == 111){
426                                 Double_t pi0Pt = prim->Pt() ;
427                                 //printf("pi0, pt %2.2f\n",pi0Pt);
428                                 Double_t pi0Y  = 0.5*TMath::Log((prim->Energy()-prim->Pz())/(prim->Energy()+prim->Pz())) ;
429                                 Double_t phi   = TMath::RadToDeg()*prim->Phi() ;
430                                 if(TMath::Abs(pi0Y) < 0.5){
431                                         fhPrimPt->Fill(pi0Pt) ;
432                                 }
433                                 fhPrimY  ->Fill(pi0Y) ;
434                                 fhPrimPhi->Fill(phi) ;
435                                 
436                                 //Check if both photons hit Calorimeter
437                                 Int_t iphot1=prim->GetFirstDaughter() ;
438                                 Int_t iphot2=prim->GetLastDaughter() ;
439                                 if(iphot1>-1 && iphot1<stack->GetNtrack() && iphot2>-1 && iphot2<stack->GetNtrack()){
440                                         TParticle * phot1 = stack->Particle(iphot1) ;
441                                         TParticle * phot2 = stack->Particle(iphot2) ;
442                                         if(phot1 && phot2 && phot1->GetPdgCode()==22 && phot2->GetPdgCode()==22){
443                                                 //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",
444                                                 //      phot1->Pt(), phot1->Phi()*180./3.1415, phot1->Eta(), phot2->Pt(), phot2->Phi()*180./3.1415, phot2->Eta());
445                                                 Bool_t inacceptance = kFALSE;
446 #ifdef __PHOSGEO__
447                                                 Int_t mod ;
448                                                 Double_t x,z ;
449                                                 if(fCalorimeter == "PHOS" && fPHOSGeo->ImpactOnEmc(phot1,mod,z,x) && fPHOSGeo->ImpactOnEmc(phot1,mod,z,x)) 
450                                                         inacceptance = kTRUE;
451                                                 //printf("In REAL PHOS acceptance? %d\n",inacceptance);
452 #else
453                                                 TLorentzVector lv1, lv2;
454                                                 phot1->Momentum(lv1);
455                                                 phot2->Momentum(lv2);
456                                                 if(GetFidutialCut()->IsInFidutialCut(lv1,fCalorimeter) && GetFidutialCut()->IsInFidutialCut(lv2,fCalorimeter)) 
457                                                         inacceptance = kTRUE ;
458                                                 //printf("In %s fidutial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
459 #endif                                                                                                    
460                                                 if(inacceptance){
461                                                         fhPrimAccPt->Fill(pi0Pt) ;
462                                                         fhPrimAccPhi->Fill(phi) ;
463                                                         fhPrimAccY->Fill(pi0Y) ;
464                                                 }//Accepted
465                                         }// 2 photons      
466                                 }//Check daughters exist
467                         }// Primary pi0
468                 }//loop on primaries    
469         }//stack exists and data is MC
470         
471 }       
472
473
474
475
476 //____________________________________________________________________________________________________________________________________________________
477 void AliAnaPi0::Terminate() 
478 {
479   //Do some calculations and plots from the final histograms.
480  
481   printf(" *** %s Terminate:", GetName()) ; 
482   
483   if (!fhRe1) {
484      Error("Terminate", "Remote output histograms not imported in AliAnaPi0 object");
485      return;
486   }
487       
488   printf("        Mgg Real        : %5.3f , RMS : %5.3f \n", fhRe1[0]->GetMean(),   fhRe1[0]->GetRMS() ) ;
489  
490   TCanvas  * cIM = new TCanvas("cIM", "", 400, 10, 600, 700) ;
491   cIM->Divide(2, 2);
492
493   cIM->cd(1) ; 
494   //gPad->SetLogy();
495   TH1D * hIMAllPt = (TH1D*) fhRe1[0]->ProjectionZ();
496   hIMAllPt->SetLineColor(2);
497   hIMAllPt->SetTitle("No cut on  p_{T, #gamma#gamma} ");
498   hIMAllPt->Draw();
499
500   cIM->cd(2) ; 
501   TH3F * hRe1Pt5 = (TH3F*)fhRe1[0]->Clone("IMPt5");
502   hRe1Pt5->GetXaxis()->SetRangeUser(0,5);
503   TH1D * hIMPt5 = (TH1D*) hRe1Pt5->Project3D("z");
504   hIMPt5->SetLineColor(2);  
505   hIMPt5->SetTitle("0 < p_{T, #gamma#gamma} < 5 GeV/c");
506   hIMPt5->Draw();
507   
508   cIM->cd(3) ; 
509   TH3F * hRe1Pt10 =  (TH3F*)fhRe1[0]->Clone("IMPt10");
510   hRe1Pt10->GetXaxis()->SetRangeUser(5,10);
511   TH1D * hIMPt10 = (TH1D*) hRe1Pt10->Project3D("z");
512   hIMPt10->SetLineColor(2);  
513   hIMPt10->SetTitle("5 < p_{T, #gamma#gamma} < 10 GeV/c");
514   hIMPt10->Draw();
515   
516   cIM->cd(4) ; 
517   TH3F * hRe1Pt20 =  (TH3F*)fhRe1[0]->Clone("IMPt20");
518   hRe1Pt20->GetXaxis()->SetRangeUser(10,20);
519   TH1D * hIMPt20 = (TH1D*) hRe1Pt20->Project3D("z");
520   hIMPt20->SetLineColor(2);  
521   hIMPt20->SetTitle("10 < p_{T, #gamma#gamma} < 20 GeV/c");
522   hIMPt20->Draw();
523    
524   cIM->Print("Mgg.eps");
525
526
527  TCanvas  * cPt = new TCanvas("cPt", "", 400, 10, 600, 700) ;
528   cPt->Divide(2, 2);
529
530   cPt->cd(1) ; 
531   //gPad->SetLogy();
532   TH1D * hPt = (TH1D*) fhRe1[0]->Project3D("x");
533   hPt->SetLineColor(2);
534   hPt->SetTitle("No cut on  M_{#gamma#gamma} ");
535   hPt->Draw();
536
537   cPt->cd(2) ; 
538   TH3F * hRe1IM1 = (TH3F*)fhRe1[0]->Clone("PtIM1");
539   hRe1IM1->GetZaxis()->SetRangeUser(0.05,0.21);
540   TH1D * hPtIM1 = (TH1D*) hRe1IM1->Project3D("x");
541   hPtIM1->SetLineColor(2);  
542   hPtIM1->SetTitle("0.05 < M_{#gamma#gamma} < 0.21 GeV/c^{2}");
543   hPtIM1->Draw();
544   
545   cPt->cd(3) ; 
546   TH3F * hRe1IM2 = (TH3F*)fhRe1[0]->Clone("PtIM2");
547   hRe1IM2->GetZaxis()->SetRangeUser(0.09,0.17);
548   TH1D * hPtIM2 = (TH1D*) hRe1IM2->Project3D("x");
549   hPtIM2->SetLineColor(2);  
550   hPtIM2->SetTitle("0.09 < M_{#gamma#gamma} < 0.17 GeV/c^{2}");
551   hPtIM2->Draw();
552
553   cPt->cd(4) ; 
554   TH3F * hRe1IM3 = (TH3F*)fhRe1[0]->Clone("PtIM3");
555   hRe1IM3->GetZaxis()->SetRangeUser(0.11,0.15);
556   TH1D * hPtIM3 = (TH1D*) hRe1IM1->Project3D("x");
557   hPtIM3->SetLineColor(2);  
558   hPtIM3->SetTitle("0.11 < M_{#gamma#gamma} < 0.15 GeV/c^{2}");
559   hPtIM3->Draw();
560    
561   cPt->Print("Pt.eps");
562
563  
564   char line[1024] ; 
565   sprintf(line, ".!tar -zcf %s.tar.gz *.eps", GetName()) ; 
566   gROOT->ProcessLine(line);
567   sprintf(line, ".!rm -fR *.eps"); 
568   gROOT->ProcessLine(line);
569  
570   printf("!! All the eps files are in %s.tar.gz !!!", GetName());
571
572 }
573
574
575
576
577