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