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