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