Protection against empty events.
[u/mrichter/AliRoot.git] / PWG4 / PartCorrDep / AliAnaPi0.cxx
CommitLineData
1c5acb87 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"
6639984f 30#include "TCanvas.h"
31#include "TPad.h"
32#include "TROOT.h"
1c5acb87 33
34//---- AliRoot system ----
35#include "AliAnaPi0.h"
36#include "AliCaloTrackReader.h"
37#include "AliCaloPID.h"
6639984f 38#include "AliStack.h"
1c5acb87 39
40#ifdef __PHOSGEO__
41 #include "AliPHOSGeoUtils.h"
42#endif
43
44ClassImp(AliAnaPi0)
45
46//________________________________________________________________________________________________________________________________________________
47AliAnaPi0::AliAnaPi0() : AliAnaPartCorrBaseClass(),
48fNCentrBin(0),fNZvertBin(0),fNrpBin(0),
49fNPID(0),fNmaxMixEv(0), fZvtxCut(0.),fCalorimeter(""),
50fEventsList(0x0), fhEtalon(0x0),
51fhRe1(0x0),fhMi1(0x0),fhRe2(0x0),fhMi2(0x0),fhRe3(0x0),fhMi3(0x0),fhEvents(0x0),
52fhPrimPt(0x0), fhPrimAccPt(0x0), fhPrimY(0x0), fhPrimAccY(0x0), fhPrimPhi(0x0), fhPrimAccPhi(0x0)
53{
54//Default Ctor
55 InitParameters();
56
57}
58
59//________________________________________________________________________________________________________________________________________________
60AliAnaPi0::AliAnaPi0(const AliAnaPi0 & ex) : AliAnaPartCorrBaseClass(ex),
61fNCentrBin(ex.fNCentrBin),fNZvertBin(ex.fNZvertBin),fNrpBin(ex.fNrpBin),
62fNPID(ex.fNPID),fNmaxMixEv(ex.fNmaxMixEv),fZvtxCut(ex.fZvtxCut), fCalorimeter(ex.fCalorimeter),
63fEventsList(ex.fEventsList), fhEtalon(ex.fhEtalon),
64fhRe1(ex.fhRe1),fhMi1(ex.fhMi1),fhRe2(ex.fhRe2),fhMi2(ex.fhMi2),fhRe3(ex.fhRe3),fhMi3(ex.fhMi3),fhEvents(ex.fhEvents),
65fhPrimPt(ex.fhPrimPt), fhPrimAccPt(ex.fhPrimAccPt), fhPrimY(ex.fhPrimY),
66fhPrimAccY(ex.fhPrimAccY), fhPrimPhi(ex.fhPrimPhi), fhPrimAccPhi(ex.fhPrimAccPhi)
67{
68 // cpy ctor
69 //Do not need it
70}
71
72//________________________________________________________________________________________________________________________________________________
73AliAnaPi0 & 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//________________________________________________________________________________________________________________________________________________
93AliAnaPi0::~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//________________________________________________________________________________________________________________________________________________
114void 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//________________________________________________________________________________________________________________________________________________
128void 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 }
2e557d1c 139
140}
1c5acb87 141
2e557d1c 142//________________________________________________________________________________________________________________________________________________
143TList * AliAnaPi0::GetCreateOutputObjects()
144{
145 // Create histograms to be saved in output file and
146 // store them in fOutputContainer
147
148 //create event containers
1c5acb87 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
1c5acb87 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//_________________________________________________________________________________________________________________________________________________
292void 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//____________________________________________________________________________________________________________________________________________________
6639984f 309void AliAnaPi0::MakeAnalysisFillHistograms()
1c5acb87 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
6639984f 476//____________________________________________________________________________________________________________________________________________________
477void AliAnaPi0::Terminate()
478{
479 //Do some calculations and plots from the final histograms.
480
481 printf(" *** %s Terminate:", GetName()) ;
2e557d1c 482
483 if (!fhRe1) {
484 Error("Terminate", "Remote output histograms not imported in AliAnaPi0 object");
485 return;
486 }
487
6639984f 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");
1c5acb87 562
6639984f 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());
1c5acb87 571
6639984f 572}
1c5acb87 573
574
575
576
577