1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
17 //_________________________________________________________________________
18 // Class to collect two-photon invariant mass distributions for
19 // extractin raw pi0 yield.
21 //-- Author: Dmitri Peressounko (RRC "KI")
22 //-- Adapted to PartCorr frame by Lamia Benhabib (SUBATECH)
23 //-- and Gustavo Conesa (INFN-Frascati)
24 //_________________________________________________________________________
27 // --- ROOT system ---
30 //#include "Riostream.h"
34 #include "TClonesArray.h"
35 #include "TObjString.h"
37 //---- AliRoot system ----
38 #include "AliAnaPi0.h"
39 #include "AliCaloTrackReader.h"
40 #include "AliCaloPID.h"
42 #include "AliFiducialCut.h"
43 #include "TParticle.h"
44 #include "AliVEvent.h"
45 #include "AliESDCaloCluster.h"
46 #include "AliESDEvent.h"
47 #include "AliAODEvent.h"
48 #include "AliNeutralMesonSelection.h"
49 #include "AliMixedEvent.h"
54 //________________________________________________________________________________________________________________________________________________
55 AliAnaPi0::AliAnaPi0() : AliAnaPartCorrBaseClass(),
56 fDoOwnMix(kFALSE),fNCentrBin(0),//fNZvertBin(0),fNrpBin(0),
57 fNPID(0),fNmaxMixEv(0), fCalorimeter(""),
58 fNModules(12), fUseAngleCut(kFALSE), fEventsList(0x0), fMultiCutAna(kFALSE),
59 fNPtCuts(0),fNAsymCuts(0), fNCellNCuts(0),fNPIDBits(0),
60 fhReMod(0x0),fhReDiffMod(0x0),
61 fhRe1(0x0), fhMi1(0x0), fhRe2(0x0), fhMi2(0x0), fhRe3(0x0), fhMi3(0x0),
62 fhReInvPt1(0x0), fhMiInvPt1(0x0), fhReInvPt2(0x0), fhMiInvPt2(0x0), fhReInvPt3(0x0), fhMiInvPt3(0x0),
63 fhRePtNCellAsymCuts(0x0), fhRePIDBits(0x0),
64 fhRePtMult(0x0), fhRePtMultAsy07(0x0) ,
65 fhEvents(0x0), fhRealOpeningAngle(0x0),fhRealCosOpeningAngle(0x0),
66 fhPrimPt(0x0), fhPrimAccPt(0x0), fhPrimY(0x0), fhPrimAccY(0x0), fhPrimPhi(0x0), fhPrimAccPhi(0x0),
67 fhPrimOpeningAngle(0x0),fhPrimCosOpeningAngle(0x0)
74 //________________________________________________________________________________________________________________________________________________
75 AliAnaPi0::~AliAnaPi0() {
76 // Remove event containers
78 if(fDoOwnMix && fEventsList){
79 for(Int_t ic=0; ic<fNCentrBin; ic++){
80 for(Int_t iz=0; iz<GetNZvertBin(); iz++){
81 for(Int_t irp=0; irp<GetNRPBin(); irp++){
82 fEventsList[ic*GetNZvertBin()*GetNRPBin()+iz*GetNRPBin()+irp]->Delete() ;
83 delete fEventsList[ic*GetNZvertBin()*GetNRPBin()+iz*GetNRPBin()+irp] ;
93 //________________________________________________________________________________________________________________________________________________
94 void AliAnaPi0::InitParameters()
96 //Init parameters when first called the analysis
97 //Set default parameters
98 SetInputAODName("PWG4Particle");
100 AddToHistogramsName("AnaPi0_");
101 fNModules = 12; // set maximum to maximum number of EMCAL modules
108 fCalorimeter = "PHOS";
109 fUseAngleCut = kFALSE;
111 fMultiCutAna = kFALSE;
114 fPtCuts[0] = 0.; fPtCuts[1] = 0.2; fPtCuts[2] = 0.3;fPtCuts[3] = 0.; fPtCuts[4] = 0.;
117 fAsymCuts[0] = 0.7; fAsymCuts[1] = 0.8; fAsymCuts[2] = 1.; fAsymCuts[3] = 0.; fAsymCuts[4] = 0.;
120 fCellNCuts[0] = 1; fCellNCuts[1] = 2; fCellNCuts[2] = 3; fCellNCuts[3] = 0; fCellNCuts[4] = 0;
123 fPIDBits[0] = 2; fPIDBits[1] = 4; fPIDBits[2] = 6; // check dispersion, neutral, dispersion&&neutral
124 fPIDBits[3] = 0; fPIDBits[4] = 0;
128 //________________________________________________________________________________________________________________________________________________
129 TObjString * AliAnaPi0::GetAnalysisCuts()
131 //Save parameters used for analysis
132 TString parList ; //this will be list of parameters used for this analysis.
133 const Int_t buffersize = 255;
134 char onePar[buffersize] ;
135 snprintf(onePar,buffersize,"--- AliAnaPi0 ---\n") ;
137 snprintf(onePar,buffersize,"Number of bins in Centrality: %d \n",fNCentrBin) ;
139 snprintf(onePar,buffersize,"Number of bins in Z vert. pos: %d \n",GetNZvertBin()) ;
141 snprintf(onePar,buffersize,"Number of bins in Reac. Plain: %d \n",GetNRPBin()) ;
143 snprintf(onePar,buffersize,"Depth of event buffer: %d \n",fNmaxMixEv) ;
145 snprintf(onePar,buffersize,"Number of different PID used: %d \n",fNPID) ;
147 snprintf(onePar,buffersize,"Cuts: \n") ;
149 snprintf(onePar,buffersize,"Z vertex position: -%f < z < %f \n",GetZvertexCut(),GetZvertexCut()) ;
151 snprintf(onePar,buffersize,"Calorimeter: %s \n",fCalorimeter.Data()) ;
153 snprintf(onePar,buffersize,"Number of modules: %d \n",fNModules) ;
156 snprintf(onePar, buffersize," pT cuts: n = %d, pt > ",fNPtCuts) ;
157 for(Int_t i = 0; i < fNPtCuts; i++) snprintf(onePar,buffersize,"%s %2.2f;",onePar,fPtCuts[i]);
160 snprintf(onePar,buffersize, " N cell in cluster cuts: n = %d, nCell > ",fNCellNCuts) ;
161 for(Int_t i = 0; i < fNCellNCuts; i++) snprintf(onePar,buffersize,"%s %d;",onePar,fCellNCuts[i]);
164 snprintf(onePar,buffersize," Asymmetry cuts: n = %d, asymmetry < ",fNAsymCuts) ;
165 for(Int_t i = 0; i < fNAsymCuts; i++) snprintf(onePar,buffersize,"%s %2.2f;",onePar,fAsymCuts[i]);
168 snprintf(onePar,buffersize," PID selection bits: n = %d, PID bit =\n",fNPIDBits) ;
169 for(Int_t i = 0; i < fNPIDBits; i++) snprintf(onePar,buffersize,"%s %d;",onePar,fPIDBits[i]);
173 return new TObjString(parList) ;
176 //________________________________________________________________________________________________________________________________________________
177 TList * AliAnaPi0::GetCreateOutputObjects()
179 // Create histograms to be saved in output file and
180 // store them in fOutputContainer
182 //create event containers
183 fEventsList = new TList*[fNCentrBin*GetNZvertBin()*GetNRPBin()] ;
185 for(Int_t ic=0; ic<fNCentrBin; ic++){
186 for(Int_t iz=0; iz<GetNZvertBin(); iz++){
187 for(Int_t irp=0; irp<GetNRPBin(); irp++){
188 fEventsList[ic*GetNZvertBin()*GetNRPBin()+iz*GetNRPBin()+irp] = new TList() ;
193 TList * outputContainer = new TList() ;
194 outputContainer->SetName(GetName());
196 fhReMod = new TH3D*[fNModules] ;
197 fhReDiffMod = new TH3D*[fNModules+1] ;
199 fhRe1 = new TH3D*[fNCentrBin*fNPID] ;
200 fhRe2 = new TH3D*[fNCentrBin*fNPID] ;
201 fhRe3 = new TH3D*[fNCentrBin*fNPID] ;
202 fhMi1 = new TH3D*[fNCentrBin*fNPID] ;
203 fhMi2 = new TH3D*[fNCentrBin*fNPID] ;
204 fhMi3 = new TH3D*[fNCentrBin*fNPID] ;
206 fhReInvPt1 = new TH3D*[fNCentrBin*fNPID] ;
207 fhReInvPt2 = new TH3D*[fNCentrBin*fNPID] ;
208 fhReInvPt3 = new TH3D*[fNCentrBin*fNPID] ;
209 fhMiInvPt1 = new TH3D*[fNCentrBin*fNPID] ;
210 fhMiInvPt2 = new TH3D*[fNCentrBin*fNPID] ;
211 fhMiInvPt3 = new TH3D*[fNCentrBin*fNPID] ;
213 const Int_t buffersize = 255;
214 char key[buffersize] ;
215 char title[buffersize] ;
217 Int_t nptbins = GetHistoPtBins();
218 Int_t nphibins = GetHistoPhiBins();
219 Int_t netabins = GetHistoEtaBins();
220 Float_t ptmax = GetHistoPtMax();
221 Float_t phimax = GetHistoPhiMax();
222 Float_t etamax = GetHistoEtaMax();
223 Float_t ptmin = GetHistoPtMin();
224 Float_t phimin = GetHistoPhiMin();
225 Float_t etamin = GetHistoEtaMin();
227 Int_t nmassbins = GetHistoMassBins();
228 Int_t nasymbins = GetHistoAsymmetryBins();
229 Float_t massmax = GetHistoMassMax();
230 Float_t asymmax = GetHistoAsymmetryMax();
231 Float_t massmin = GetHistoMassMin();
232 Float_t asymmin = GetHistoAsymmetryMin();
234 for(Int_t ic=0; ic<fNCentrBin; ic++){
235 for(Int_t ipid=0; ipid<fNPID; ipid++){
237 //Distance to bad module 1
238 snprintf(key, buffersize,"hRe_cen%d_pid%d_dist1",ic,ipid) ;
239 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
240 fhRe1[ic*fNPID+ipid] = new TH3D(key,title,nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax,nmassbins,massmin,massmax) ;
241 outputContainer->Add(fhRe1[ic*fNPID+ipid]) ;
243 //Distance to bad module 2
244 snprintf(key, buffersize,"hRe_cen%d_pid%d_dist2",ic,ipid) ;
245 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
246 fhRe2[ic*fNPID+ipid] = new TH3D(key,title,nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax,nmassbins,massmin,massmax) ;
247 outputContainer->Add(fhRe2[ic*fNPID+ipid]) ;
249 //Distance to bad module 3
250 snprintf(key, buffersize,"hRe_cen%d_pid%d_dist3",ic,ipid) ;
251 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
252 fhRe3[ic*fNPID+ipid] = new TH3D(key,title,nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax,nmassbins,massmin,massmax) ;
253 outputContainer->Add(fhRe3[ic*fNPID+ipid]) ;
256 //Distance to bad module 1
257 snprintf(key, buffersize,"hReInvPt_cen%d_pid%d_dist1",ic,ipid) ;
258 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
259 fhReInvPt1[ic*fNPID+ipid] = new TH3D(key,title,nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax,nmassbins,massmin,massmax) ;
260 outputContainer->Add(fhReInvPt1[ic*fNPID+ipid]) ;
262 //Distance to bad module 2
263 snprintf(key, buffersize,"hReInvPt_cen%d_pid%d_dist2",ic,ipid) ;
264 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
265 fhReInvPt2[ic*fNPID+ipid] = new TH3D(key,title,nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax,nmassbins,massmin,massmax) ;
266 outputContainer->Add(fhReInvPt2[ic*fNPID+ipid]) ;
268 //Distance to bad module 3
269 snprintf(key, buffersize,"hReInvPt_cen%d_pid%d_dist3",ic,ipid) ;
270 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
271 fhReInvPt3[ic*fNPID+ipid] = new TH3D(key,title,nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax,nmassbins,massmin,massmax) ;
272 outputContainer->Add(fhReInvPt3[ic*fNPID+ipid]) ;
275 //Distance to bad module 1
276 snprintf(key, buffersize,"hMi_cen%d_pid%d_dist1",ic,ipid) ;
277 snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
278 fhMi1[ic*fNPID+ipid] = new TH3D(key,title,nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax,nmassbins,massmin,massmax) ;
279 outputContainer->Add(fhMi1[ic*fNPID+ipid]) ;
281 //Distance to bad module 2
282 snprintf(key, buffersize,"hMi_cen%d_pid%d_dist2",ic,ipid) ;
283 snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
284 fhMi2[ic*fNPID+ipid] = new TH3D(key,title,nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax,nmassbins,massmin,massmax) ;
285 outputContainer->Add(fhMi2[ic*fNPID+ipid]) ;
287 //Distance to bad module 3
288 snprintf(key, buffersize,"hMi_cen%d_pid%d_dist3",ic,ipid) ;
289 snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
290 fhMi3[ic*fNPID+ipid] = new TH3D(key,title,nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax,nmassbins,massmin,massmax) ;
291 outputContainer->Add(fhMi3[ic*fNPID+ipid]) ;
294 //Distance to bad module 1
295 snprintf(key, buffersize,"hMiInvPt_cen%d_pid%d_dist1",ic,ipid) ;
296 snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
297 fhMiInvPt1[ic*fNPID+ipid] = new TH3D(key,title,nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax,nmassbins,massmin,massmax) ;
298 outputContainer->Add(fhMiInvPt1[ic*fNPID+ipid]) ;
300 //Distance to bad module 2
301 snprintf(key, buffersize,"hMiInvPt_cen%d_pid%d_dist2",ic,ipid) ;
302 snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
303 fhMiInvPt2[ic*fNPID+ipid] = new TH3D(key,title,nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax,nmassbins,massmin,massmax) ;
304 outputContainer->Add(fhMiInvPt2[ic*fNPID+ipid]) ;
306 //Distance to bad module 3
307 snprintf(key, buffersize,"hMiInvPt_cen%d_pid%d_dist3",ic,ipid) ;
308 snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
309 fhMiInvPt3[ic*fNPID+ipid] = new TH3D(key,title,nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax,nmassbins,massmin,massmax) ;
310 outputContainer->Add(fhMiInvPt3[ic*fNPID+ipid]) ;
319 fhRePIDBits = new TH2D*[fNPIDBits];
320 for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
321 snprintf(key, buffersize,"hRe_pidbit%d",ipid) ;
322 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for PIDBit=%d",fPIDBits[ipid]) ;
323 fhRePIDBits[ipid] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
324 outputContainer->Add(fhRePIDBits[ipid]) ;
327 fhRePtNCellAsymCuts = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
328 for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
329 for(Int_t icell=0; icell<fNCellNCuts; icell++){
330 for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
331 snprintf(key, buffersize,"hRe_pt%d_cell%d_asym%d",ipt,icell,iasym) ;
332 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for pt >%2.2f, ncell>%d and asym >%2.2f ",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]) ;
333 Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
334 //printf("ipt %d, icell %d, iassym %d, index %d\n",ipt, icell, iasym, index);
335 fhRePtNCellAsymCuts[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
336 outputContainer->Add(fhRePtNCellAsymCuts[index]) ;
341 fhRePtMult = new TH3D("hRePtMult","(p_{T},C,M)_{#gamma#gamma}, 0<A<1.0", nptbins,ptmin,ptmax,40,0.,160.,nmassbins,massmin,massmax);
342 fhRePtMultAsy07 = new TH3D("hRePtMultAsy07","(p_{T},C,M)_{#gamma#gamma}, 0<A<0.7",nptbins,ptmin,ptmax,40,0.,160.,nmassbins,massmin,massmax);
343 outputContainer->Add(fhRePtMult) ;
344 outputContainer->Add(fhRePtMultAsy07) ;
346 }// multi cuts analysis
348 fhEvents=new TH3D("hEvents","Number of events",fNCentrBin,0.,1.*fNCentrBin,
349 GetNZvertBin(),0.,1.*GetNZvertBin(),GetNRPBin(),0.,1.*GetNRPBin()) ;
350 outputContainer->Add(fhEvents) ;
352 fhRealOpeningAngle = new TH2D
353 ("hRealOpeningAngle","Angle between all #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,200,0,0.5);
354 fhRealOpeningAngle->SetYTitle("#theta(rad)");
355 fhRealOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
356 outputContainer->Add(fhRealOpeningAngle) ;
358 fhRealCosOpeningAngle = new TH2D
359 ("hRealCosOpeningAngle","Cosinus of angle between all #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,200,-1,1);
360 fhRealCosOpeningAngle->SetYTitle("cos (#theta) ");
361 fhRealCosOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
362 outputContainer->Add(fhRealCosOpeningAngle) ;
365 //Histograms filled only if MC data is requested
368 fhPrimPt = new TH1D("hPrimPt","Primary pi0 pt",nptbins,ptmin,ptmax) ;
369 fhPrimAccPt = new TH1D("hPrimAccPt","Primary pi0 pt with both photons in acceptance",nptbins,ptmin,ptmax) ;
370 outputContainer->Add(fhPrimPt) ;
371 outputContainer->Add(fhPrimAccPt) ;
373 fhPrimY = new TH1D("hPrimaryRapidity","Rapidity of primary pi0",netabins,etamin,etamax) ;
374 outputContainer->Add(fhPrimY) ;
376 fhPrimAccY = new TH1D("hPrimAccRapidity","Rapidity of primary pi0",netabins,etamin,etamax) ;
377 outputContainer->Add(fhPrimAccY) ;
379 fhPrimPhi = new TH1D("hPrimaryPhi","Azimithal of primary pi0",nphibins,phimin*TMath::RadToDeg(),phimax*TMath::RadToDeg()) ;
380 outputContainer->Add(fhPrimPhi) ;
382 fhPrimAccPhi = new TH1D("hPrimAccPhi","Azimithal of primary pi0 with accepted daughters",nphibins,phimin*TMath::RadToDeg(),phimax*TMath::RadToDeg()) ;
383 outputContainer->Add(fhPrimAccPhi) ;
386 fhPrimOpeningAngle = new TH2D
387 ("hPrimOpeningAngle","Angle between all primary #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,100,0,0.5);
388 fhPrimOpeningAngle->SetYTitle("#theta(rad)");
389 fhPrimOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
390 outputContainer->Add(fhPrimOpeningAngle) ;
392 fhPrimCosOpeningAngle = new TH2D
393 ("hPrimCosOpeningAngle","Cosinus of angle between all primary #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,100,-1,1);
394 fhPrimCosOpeningAngle->SetYTitle("cos (#theta) ");
395 fhPrimCosOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
396 outputContainer->Add(fhPrimCosOpeningAngle) ;
400 TString * pairname = new TString[fNModules];
401 if(fCalorimeter=="EMCAL"){
402 pairname[0]="A side (0-2)";
403 pairname[1]="C side (1-3)";
404 pairname[2]="Sector 0 (0-1)";
405 pairname[3]="Sector 1 (2-3)";
406 for(Int_t i = 4 ; i < fNModules ; i++) pairname[i]="";}
407 if(fCalorimeter=="PHOS") {
411 for(Int_t i = 3 ; i < fNModules ; i++) pairname[i]="";}
413 for(Int_t imod=0; imod<fNModules; imod++){
414 //Module dependent invariant mass
415 snprintf(key, buffersize,"hReMod_%d",imod) ;
416 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for Module %d",imod) ;
417 fhReMod[imod] = new TH3D(key,title,nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax,nmassbins,massmin,massmax) ;
418 outputContainer->Add(fhReMod[imod]) ;
420 snprintf(key, buffersize,"hReDiffMod_%d",imod) ;
421 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for Different Modules: %s",(pairname[imod]).Data()) ;
422 fhReDiffMod[imod] = new TH3D(key,title,nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax,nmassbins,massmin,massmax) ;
423 outputContainer->Add(fhReDiffMod[imod]) ;
428 snprintf(key, buffersize,"hReDiffMod_%d",fNModules) ;
429 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for all Modules Combination") ;
430 fhReDiffMod[fNModules] = new TH3D(key,title,nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax,nmassbins,massmin,massmax) ;
431 outputContainer->Add(fhReDiffMod[fNModules]) ;
434 // for(Int_t i = 0; i < outputContainer->GetEntries() ; i++){
436 // printf("Histogram %d, name: %s\n ",i, outputContainer->At(i)->GetName());
440 return outputContainer;
443 //_________________________________________________________________________________________________________________________________________________
444 void AliAnaPi0::Print(const Option_t * /*opt*/) const
446 //Print some relevant parameters set for the analysis
447 printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
448 AliAnaPartCorrBaseClass::Print(" ");
450 printf("Number of bins in Centrality: %d \n",fNCentrBin) ;
451 printf("Number of bins in Z vert. pos: %d \n",GetNZvertBin()) ;
452 printf("Number of bins in Reac. Plain: %d \n",GetNRPBin()) ;
453 printf("Depth of event buffer: %d \n",fNmaxMixEv) ;
454 printf("Number of different PID used: %d \n",fNPID) ;
456 printf("Z vertex position: -%2.3f < z < %2.3f \n",GetZvertexCut(),GetZvertexCut()) ;
457 printf("Number of modules: %d \n",fNModules) ;
458 printf("Select pairs with their angle: %d \n",fUseAngleCut) ;
460 printf("pT cuts: n = %d, \n",fNPtCuts) ;
462 for(Int_t i = 0; i < fNPtCuts; i++) printf("%2.2f ",fPtCuts[i]);
465 printf("N cell in cluster cuts: n = %d, \n",fNCellNCuts) ;
466 printf("\tnCell > ");
467 for(Int_t i = 0; i < fNCellNCuts; i++) printf("%d ",fCellNCuts[i]);
470 printf("Asymmetry cuts: n = %d, \n",fNAsymCuts) ;
471 printf("\tasymmetry < ");
472 for(Int_t i = 0; i < fNAsymCuts; i++) printf("%2.2f ",fAsymCuts[i]);
475 printf("PID selection bits: n = %d, \n",fNPIDBits) ;
476 printf("\tPID bit = ");
477 for(Int_t i = 0; i < fNPIDBits; i++) printf("%d ",fPIDBits[i]);
481 printf("------------------------------------------------------\n") ;
484 //_____________________________________________________________
485 void AliAnaPi0::FillAcceptanceHistograms(){
486 //Fill acceptance histograms if MC data is available
488 if(IsDataMC() && GetReader()->ReadStack()){
489 AliStack * stack = GetMCStack();
490 if(stack && (IsDataMC() || (GetReader()->GetDataType() == AliCaloTrackReader::kMC)) ){
491 for(Int_t i=0 ; i<stack->GetNprimary(); i++){
492 TParticle * prim = stack->Particle(i) ;
493 if(prim->GetPdgCode() == 111){
494 Double_t pi0Pt = prim->Pt() ;
495 //printf("pi0, pt %2.2f\n",pi0Pt);
496 if(prim->Energy() == TMath::Abs(prim->Pz())) continue ; //Protection against floating point exception
497 Double_t pi0Y = 0.5*TMath::Log((prim->Energy()-prim->Pz())/(prim->Energy()+prim->Pz())) ;
498 Double_t phi = TMath::RadToDeg()*prim->Phi() ;
499 if(TMath::Abs(pi0Y) < 0.5){
500 fhPrimPt->Fill(pi0Pt) ;
502 fhPrimY ->Fill(pi0Y) ;
503 fhPrimPhi->Fill(phi) ;
505 //Check if both photons hit Calorimeter
506 Int_t iphot1=prim->GetFirstDaughter() ;
507 Int_t iphot2=prim->GetLastDaughter() ;
508 if(iphot1>-1 && iphot1<stack->GetNtrack() && iphot2>-1 && iphot2<stack->GetNtrack()){
509 TParticle * phot1 = stack->Particle(iphot1) ;
510 TParticle * phot2 = stack->Particle(iphot2) ;
511 if(phot1 && phot2 && phot1->GetPdgCode()==22 && phot2->GetPdgCode()==22){
512 //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",
513 // phot1->Pt(), phot1->Phi()*180./3.1415, phot1->Eta(), phot2->Pt(), phot2->Phi()*180./3.1415, phot2->Eta());
515 TLorentzVector lv1, lv2;
516 phot1->Momentum(lv1);
517 phot2->Momentum(lv2);
519 Bool_t inacceptance = kFALSE;
520 if(fCalorimeter == "PHOS"){
521 if(GetPHOSGeometry() && GetCaloUtils()->IsPHOSGeoMatrixSet()){
524 if(GetPHOSGeometry()->ImpactOnEmc(phot1,mod,z,x) && GetPHOSGeometry()->ImpactOnEmc(phot2,mod,z,x))
525 inacceptance = kTRUE;
526 if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
530 if(GetFiducialCut()->IsInFiducialCut(lv1,fCalorimeter) && GetFiducialCut()->IsInFiducialCut(lv2,fCalorimeter))
531 inacceptance = kTRUE ;
532 if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
536 else if(fCalorimeter == "EMCAL" && GetCaloUtils()->IsEMCALGeoMatrixSet()){
537 if(GetEMCALGeometry()){
538 if(GetEMCALGeometry()->Impact(phot1) && GetEMCALGeometry()->Impact(phot2))
539 inacceptance = kTRUE;
540 if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
543 if(GetFiducialCut()->IsInFiducialCut(lv1,fCalorimeter) && GetFiducialCut()->IsInFiducialCut(lv2,fCalorimeter))
544 inacceptance = kTRUE ;
545 if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
551 fhPrimAccPt->Fill(pi0Pt) ;
552 fhPrimAccPhi->Fill(phi) ;
553 fhPrimAccY->Fill(pi0Y) ;
554 Double_t angle = lv1.Angle(lv2.Vect());
555 fhPrimOpeningAngle ->Fill(pi0Pt,angle);
556 fhPrimCosOpeningAngle->Fill(pi0Pt,TMath::Cos(angle));
560 }//Check daughters exist
563 }//stack exists and data is MC
565 else if(GetReader()->ReadAODMCParticles()){
566 if(GetDebug() >= 0) printf("AliAnaPi0::FillAcceptanceHistograms() - Acceptance calculation with MCParticles not implemented yet\n");
570 //____________________________________________________________________________________________________________________________________________________
571 void AliAnaPi0::MakeAnalysisFillHistograms()
573 //Process one event and extract photons from AOD branch
574 // filled with AliAnaPhoton and fill histos with invariant mass
576 //In case of MC data, fill acceptance histograms
577 FillAcceptanceHistograms();
579 //Apply some cuts on event: vertex position and centrality range
580 Int_t iRun=(GetReader()->GetInputEvent())->GetRunNumber() ;
581 if(IsBadRun(iRun)) return ;
583 Int_t nPhot = GetInputAODBranch()->GetEntriesFast() ;
585 printf("AliAnaPi0::MakeAnalysisFillHistograms() - Photon entries %d\n", nPhot);
590 Double_t vert[] = {0.0, 0.0, 0.0} ; //vertex
591 Int_t evtIndex1 = 0 ;
592 Int_t currentEvtIndex = -1 ;
593 Int_t curCentrBin = 0 ;
595 Int_t curZvertBin = 0 ;
597 for(Int_t i1=0; i1<nPhot-1; i1++){
598 AliAODPWG4Particle * p1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i1)) ;
599 // get the event index in the mixed buffer where the photon comes from
600 // in case of mixing with analysis frame, not own mixing
601 evtIndex1 = GetEventIndex(p1, vert) ;
602 //printf("charge = %d\n", track->Charge());
603 if ( evtIndex1 == -1 )
605 if ( evtIndex1 == -2 )
607 if(TMath::Abs(vert[2]) > GetZvertexCut()) continue ; //vertex cut
608 if (evtIndex1 != currentEvtIndex) {
609 //Get Reaction Plan position and calculate RP bin
610 //does not exist in ESD yet????
613 curZvertBin = (Int_t)(0.5*GetNZvertBin()*(vert[2]+GetZvertexCut())/GetZvertexCut()) ;
614 fhEvents->Fill(curCentrBin+0.5,curZvertBin+0.5,curRPBin+0.5) ;
615 currentEvtIndex = evtIndex1 ;
618 //printf("AliAnaPi0::MakeAnalysisFillHistograms(): Photon 1 Evt %d Vertex : %f,%f,%f\n",evtIndex1, GetVertex(evtIndex1)[0] ,GetVertex(evtIndex1)[1],GetVertex(evtIndex1)[2]);
620 TLorentzVector photon1(p1->Px(),p1->Py(),p1->Pz(),p1->E());
622 module1 = GetModuleNumber(p1);
623 for(Int_t i2=i1+1; i2<nPhot; i2++){
624 AliAODPWG4Particle * p2 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i2)) ;
625 Int_t evtIndex2 = GetEventIndex(p2, vert) ;
626 if ( evtIndex2 == -1 )
628 if ( evtIndex2 == -2 )
630 if (GetMixedEvent() && (evtIndex1 == evtIndex2))
632 //printf("AliAnaPi0::MakeAnalysisFillHistograms(): Photon 2 Evt %d Vertex : %f,%f,%f\n",evtIndex2, GetVertex(evtIndex2)[0] ,GetVertex(evtIndex2)[1],GetVertex(evtIndex2)[2]);
633 TLorentzVector photon2(p2->Px(),p2->Py(),p2->Pz(),p2->E());
635 module2 = GetModuleNumber(p2);
636 Double_t m = (photon1 + photon2).M() ;
637 Double_t pt = (photon1 + photon2).Pt();
638 Double_t a = TMath::Abs(p1->E()-p2->E())/(p1->E()+p2->E()) ;
640 printf("AliAnaPi0::MakeAnalysisFillHistograms() - Current Event: pT: photon1 %2.2f, photon2 %2.2f; Pair: pT %2.2f, mass %2.3f, a %f2.3\n",
641 p1->Pt(), p2->Pt(), pt,m,a);
642 //Check if opening angle is too large or too small compared to what is expected
643 Double_t angle = photon1.Angle(photon2.Vect());
644 //if(fUseAngleCut && !GetNeutralMesonSelection()->IsAngleInWindow((photon1+photon2).E(),angle)) continue;
645 //printf("angle %f\n",angle);
646 if(fUseAngleCut && angle < 0.1)
648 fhRealOpeningAngle ->Fill(pt,angle);
649 fhRealCosOpeningAngle->Fill(pt,TMath::Cos(angle));
651 //Fill module dependent histograms
652 //if(module1==module2) printf("mod1 %d\n",module1);
653 if(module1==module2 && module1 >=0 && module1<fNModules)
654 fhReMod[module1]->Fill(pt,a,m) ;
656 fhReDiffMod[fNModules]->Fill(pt,a,m) ;
658 if(fCalorimeter=="EMCAL"){
659 if((module1==0 && module2==2) || (module1==2 && module2==0)) fhReDiffMod[0]->Fill(pt,a,m) ;
660 if((module1==1 && module2==3) || (module1==3 && module2==1)) fhReDiffMod[1]->Fill(pt,a,m) ;
661 if((module1==0 && module2==1) || (module1==1 && module2==0)) fhReDiffMod[2]->Fill(pt,a,m) ;
662 if((module1==2 && module2==3) || (module1==3 && module2==2)) fhReDiffMod[3]->Fill(pt,a,m) ;
665 if((module1==0 && module2==1) || (module1==1 && module2==0)) fhReDiffMod[0]->Fill(pt,a,m) ;
666 if((module1==0 && module2==2) || (module1==2 && module2==0)) fhReDiffMod[1]->Fill(pt,a,m) ;
667 if((module1==1 && module2==2) || (module1==2 && module2==1)) fhReDiffMod[2]->Fill(pt,a,m) ;
670 for(Int_t ipid=0; ipid<fNPID; ipid++){
671 if((p1->IsPIDOK(ipid,AliCaloPID::kPhoton)) && (p2->IsPIDOK(ipid,AliCaloPID::kPhoton))){
672 fhRe1 [curCentrBin*fNPID+ipid]->Fill(pt,a,m) ;
673 fhReInvPt1[curCentrBin*fNPID+ipid]->Fill(pt,a,m,1./pt) ;
674 if(p1->DistToBad()>0 && p2->DistToBad()>0){
675 fhRe2 [curCentrBin*fNPID+ipid]->Fill(pt,a,m) ;
676 fhReInvPt2[curCentrBin*fNPID+ipid]->Fill(pt,a,m,1./pt) ;
677 if(p1->DistToBad()>1 && p2->DistToBad()>1){
678 fhRe3 [curCentrBin*fNPID+ipid]->Fill(pt,a,m) ;
679 fhReInvPt3[curCentrBin*fNPID+ipid]->Fill(pt,a,m,1./pt) ;
685 //Multi cuts analysis
687 //Histograms for different PID bits selection
688 for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
690 if(p1->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton) &&
691 p2->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton)) fhRePIDBits[ipid]->Fill(pt,m) ;
693 //printf("ipt %d, ipid%d, name %s\n",ipt, ipid, fhRePtPIDCuts[ipt*fNPIDBits+ipid]->GetName());
694 } // pid bit cut loop
696 //Several pt,ncell and asymetry cuts
697 //Get the number of cells
700 AliVEvent * event = GetReader()->GetInputEvent();
702 for(Int_t iclus = 0; iclus < event->GetNumberOfCaloClusters(); iclus++){
703 AliVCluster *cluster = event->GetCaloCluster(iclus);
706 if (fCalorimeter == "EMCAL" && GetReader()->IsEMCALCluster(cluster)) is = kTRUE;
707 else if(fCalorimeter == "PHOS" && GetReader()->IsPHOSCluster (cluster)) is = kTRUE;
710 if (p1->GetCaloLabel(0) == cluster->GetID()) ncell1 = cluster->GetNCells();
711 else if (p2->GetCaloLabel(0) == cluster->GetID()) ncell2 = cluster->GetNCells();
712 } // PHOS or EMCAL cluster as requested in analysis
714 if(ncell2 > 0 && ncell1 > 0) break; // No need to continue the iteration
717 //printf("e 1: %2.2f, e 2: %2.2f, ncells: n1 %d, n2 %d\n", p1->E(), p2->E(),ncell1,ncell2);
719 for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
720 for(Int_t icell=0; icell<fNCellNCuts; icell++){
721 for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
723 if(p1->Pt() > fPtCuts[ipt] && p2->Pt() > fPtCuts[ipt] &&
724 a < fAsymCuts[iasym] &&
725 ncell1 >= fCellNCuts[icell] && ncell2 >= fCellNCuts[icell]) fhRePtNCellAsymCuts[((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym]->Fill(pt,m) ;
727 //printf("ipt %d, icell%d, iasym %d, name %s\n",ipt, icell, iasym, fhRePtNCellAsymCuts[((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym]->GetName());
732 fhRePtMult->Fill(pt,GetTrackMultiplicity(),m) ;
734 fhRePtMultAsy07->Fill(pt,GetTrackMultiplicity(),m) ;
736 }// multiple cuts analysis
737 }// second same event particle
742 TList * evMixList=fEventsList[curCentrBin*GetNZvertBin()*GetNRPBin()+curZvertBin*GetNRPBin()+curRPBin] ;
743 Int_t nMixed = evMixList->GetSize() ;
744 for(Int_t ii=0; ii<nMixed; ii++){
745 TClonesArray* ev2= (TClonesArray*) (evMixList->At(ii));
746 Int_t nPhot2=ev2->GetEntriesFast() ;
748 if(GetDebug() > 1) printf("AliAnaPi0::MakeAnalysisFillHistograms() - Mixed event %d photon entries %d\n", ii, nPhot);
750 for(Int_t i1=0; i1<nPhot; i1++){
751 AliAODPWG4Particle * p1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i1)) ;
752 TLorentzVector photon1(p1->Px(),p1->Py(),p1->Pz(),p1->E());
753 for(Int_t i2=0; i2<nPhot2; i2++){
754 AliAODPWG4Particle * p2 = (AliAODPWG4Particle*) (ev2->At(i2)) ;
756 TLorentzVector photon2(p2->Px(),p2->Py(),p2->Pz(),p2->E());
757 m = (photon1+photon2).M() ;
758 Double_t pt = (photon1 + photon2).Pt();
759 Double_t a = TMath::Abs(p1->E()-p2->E())/(p1->E()+p2->E()) ;
761 //Check if opening angle is too large or too small compared to what is expected
762 Double_t angle = photon1.Angle(photon2.Vect());
763 //if(fUseAngleCut && !GetNeutralMesonSelection()->IsAngleInWindow((photon1+photon2).E(),angle)) continue;
764 if(fUseAngleCut && angle < 0.1) continue;
767 printf("AliAnaPi0::MakeAnalysisFillHistograms() - Mixed Event: pT: photon1 %2.2f, photon2 %2.2f; Pair: pT %2.2f, mass %2.3f, a %f2.3\n",
768 p1->Pt(), p2->Pt(), pt,m,a);
769 for(Int_t ipid=0; ipid<fNPID; ipid++){
770 if((p1->IsPIDOK(ipid,AliCaloPID::kPhoton)) && (p2->IsPIDOK(ipid,AliCaloPID::kPhoton))){
771 fhMi1 [curCentrBin*fNPID+ipid]->Fill(pt, a,m) ;
772 fhMiInvPt1[curCentrBin*fNPID+ipid]->Fill(1./pt,a,m) ;
773 if(p1->DistToBad()>0 && p2->DistToBad()>0){
774 fhMi2 [curCentrBin*fNPID+ipid]->Fill(pt, a,m) ;
775 fhMiInvPt2[curCentrBin*fNPID+ipid]->Fill(1./pt,a,m) ;
776 if(p1->DistToBad()>1 && p2->DistToBad()>1){
777 fhMi3 [curCentrBin*fNPID+ipid]->Fill(pt, a,m) ;
778 fhMiInvPt3[curCentrBin*fNPID+ipid]->Fill(1./pt,a,m) ;
783 }//loop for histograms
784 }// second cluster loop
785 }//first cluster loop
786 }//loop on mixed events
788 TClonesArray *currentEvent = new TClonesArray(*GetInputAODBranch());
789 //Ad d current event to buffer and Remove redundant events
790 if(currentEvent->GetEntriesFast()>0){
791 evMixList->AddFirst(currentEvent) ;
792 currentEvent=0 ; //Now list of particles belongs to buffer and it will be deleted with buffer
793 if(evMixList->GetSize()>=fNmaxMixEv)
795 TClonesArray * tmp = (TClonesArray*) (evMixList->Last()) ;
796 evMixList->RemoveLast() ;
801 delete currentEvent ;
808 //________________________________________________________________________
809 void AliAnaPi0::ReadHistograms(TList* outputList)
811 // Needed when Terminate is executed in distributed environment
812 // Refill analysis histograms of this class with corresponding histograms in output list.
814 // Histograms of this analsys are kept in the same list as other analysis, recover the position of
815 // the first one and then add the next.
816 Int_t index = outputList->IndexOf(outputList->FindObject(GetAddedHistogramsStringToName()+"hRe_cen0_pid0_dist1"));
818 if(!fhRe1) fhRe1 = new TH3D*[fNCentrBin*fNPID] ;
819 if(!fhRe2) fhRe2 = new TH3D*[fNCentrBin*fNPID] ;
820 if(!fhRe3) fhRe3 = new TH3D*[fNCentrBin*fNPID] ;
821 if(!fhMi1) fhMi1 = new TH3D*[fNCentrBin*fNPID] ;
822 if(!fhMi2) fhMi2 = new TH3D*[fNCentrBin*fNPID] ;
823 if(!fhMi3) fhMi3 = new TH3D*[fNCentrBin*fNPID] ;
824 if(!fhReInvPt1) fhReInvPt1 = new TH3D*[fNCentrBin*fNPID] ;
825 if(!fhReInvPt2) fhReInvPt2 = new TH3D*[fNCentrBin*fNPID] ;
826 if(!fhReInvPt3) fhReInvPt3 = new TH3D*[fNCentrBin*fNPID] ;
827 if(!fhMiInvPt1) fhMiInvPt1 = new TH3D*[fNCentrBin*fNPID] ;
828 if(!fhMiInvPt2) fhMiInvPt2 = new TH3D*[fNCentrBin*fNPID] ;
829 if(!fhMiInvPt3) fhMiInvPt3 = new TH3D*[fNCentrBin*fNPID] ;
830 if(!fhReMod) fhReMod = new TH3D*[fNModules] ;
831 if(!fhReDiffMod)fhReDiffMod = new TH3D*[fNModules+1] ;
834 for(Int_t ic=0; ic<fNCentrBin; ic++){
835 for(Int_t ipid=0; ipid<fNPID; ipid++){
836 fhRe1[ic*fNPID+ipid] = (TH3D*) outputList->At(index++);
837 fhRe2[ic*fNPID+ipid] = (TH3D*) outputList->At(index++);
838 fhRe3[ic*fNPID+ipid] = (TH3D*) outputList->At(index++);
840 fhReInvPt1[ic*fNPID+ipid] = (TH3D*) outputList->At(index++);
841 fhReInvPt2[ic*fNPID+ipid] = (TH3D*) outputList->At(index++);
842 fhReInvPt3[ic*fNPID+ipid] = (TH3D*) outputList->At(index++);
844 fhMi1[ic*fNPID+ipid] = (TH3D*) outputList->At(index++);
845 fhMi2[ic*fNPID+ipid] = (TH3D*) outputList->At(index++);
846 fhMi3[ic*fNPID+ipid] = (TH3D*) outputList->At(index++);
848 fhMiInvPt1[ic*fNPID+ipid] = (TH3D*) outputList->At(index++);
849 fhMiInvPt2[ic*fNPID+ipid] = (TH3D*) outputList->At(index++);
850 fhMiInvPt3[ic*fNPID+ipid] = (TH3D*) outputList->At(index++);
856 if(!fhRePtNCellAsymCuts) fhRePtNCellAsymCuts = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
857 if(!fhRePIDBits) fhRePIDBits = new TH2D*[fNPIDBits];
859 for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
860 fhRePIDBits[ipid] = (TH2D*) outputList->At(index++);
863 for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
864 for(Int_t icell=0; icell<fNCellNCuts; icell++){
865 for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
866 fhRePtNCellAsymCuts[((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym] = (TH2D*) outputList->At(index++);
870 fhRePtMult = (TH3D*) outputList->At(index++);
871 fhRePtMultAsy07 = (TH3D*) outputList->At(index++);
872 }// multi cut analysis
874 fhEvents = (TH3D *) outputList->At(index++);
876 //Histograms filled only if MC data is requested
877 if(IsDataMC() || (GetReader()->GetDataType() == AliCaloTrackReader::kMC) ){
878 fhPrimPt = (TH1D*) outputList->At(index++);
879 fhPrimAccPt = (TH1D*) outputList->At(index++);
880 fhPrimY = (TH1D*) outputList->At(index++);
881 fhPrimAccY = (TH1D*) outputList->At(index++);
882 fhPrimPhi = (TH1D*) outputList->At(index++);
883 fhPrimAccPhi = (TH1D*) outputList->At(index++);
886 for(Int_t imod=0; imod < fNModules; imod++)
887 fhReMod[imod] = (TH3D*) outputList->At(index++);
893 //____________________________________________________________________________________________________________________________________________________
894 void AliAnaPi0::Terminate(TList* outputList)
896 //Do some calculations and plots from the final histograms.
898 printf(" *** %s Terminate:\n", GetName()) ;
900 //Recover histograms from output histograms list, needed for distributed analysis.
901 ReadHistograms(outputList);
904 printf("AliAnaPi0::Terminate() - Error: Remote output histograms not imported in AliAnaPi0 object");
908 printf("AliAnaPi0::Terminate() Mgg Real : %5.3f , RMS : %5.3f \n", fhRe1[0]->GetMean(), fhRe1[0]->GetRMS() ) ;
910 const Int_t buffersize = 255;
912 char nameIM[buffersize];
913 snprintf(nameIM, buffersize,"AliAnaPi0_%s_cPt",fCalorimeter.Data());
914 TCanvas * cIM = new TCanvas(nameIM, "", 400, 10, 600, 700) ;
919 TH1D * hIMAllPt = (TH1D*) fhRe1[0]->ProjectionZ(Form("IMPtAll_%s",fCalorimeter.Data()));
920 hIMAllPt->SetLineColor(2);
921 hIMAllPt->SetTitle("No cut on p_{T, #gamma#gamma} ");
925 TH1D * hIMPt5 = (TH1D*) fhRe1[0]->ProjectionZ(Form("IMPt0-5_%s",fCalorimeter.Data()),0, fhRe1[0]->GetXaxis()->FindBin(5.),0, -1, "");
926 // hRe1Pt5->GetXaxis()->SetRangeUser(0,5);
927 // TH1D * hIMPt5 = (TH1D*) hRe1Pt5->Project3D(Form("IMPt5_%s_pz",fCalorimeter.Data()));
928 hIMPt5->SetLineColor(2);
929 hIMPt5->SetTitle("0 < p_{T, #gamma#gamma} < 5 GeV/c");
933 TH1D * hIMPt10 = (TH1D*) fhRe1[0]->ProjectionZ(Form("IMPt5-10_%s",fCalorimeter.Data()), fhRe1[0]->GetXaxis()->FindBin(5.),fhRe1[0]->GetXaxis()->FindBin(10.),0, -1,"");
934 // hRe1Pt10->GetXaxis()->SetRangeUser(5,10);
935 // TH1D * hIMPt10 = (TH1D*) hRe1Pt10->Project3D(Form("IMPt10_%s_pz",fCalorimeter.Data()));
936 hIMPt10->SetLineColor(2);
937 hIMPt10->SetTitle("5 < p_{T, #gamma#gamma} < 10 GeV/c");
941 TH1D * hIMPt20 = (TH1D*) fhRe1[0]->ProjectionZ(Form("IMPt10-20_%s",fCalorimeter.Data()), fhRe1[0]->GetXaxis()->FindBin(10.),fhRe1[0]->GetXaxis()->FindBin(20.),0, -1,"");
942 // TH3F * hRe1Pt20 = (TH3F*)fhRe1[0]->Clone(Form("IMPt20_%s",fCalorimeter.Data()));
943 // hRe1Pt20->GetXaxis()->SetRangeUser(10,20);
944 // TH1D * hIMPt20 = (TH1D*) hRe1Pt20->Project3D(Form("IMPt20_%s_pz",fCalorimeter.Data()));
945 hIMPt20->SetLineColor(2);
946 hIMPt20->SetTitle("10 < p_{T, #gamma#gamma} < 20 GeV/c");
949 char nameIMF[buffersize];
950 snprintf(nameIMF,buffersize,"AliAnaPi0_%s_Mgg.eps",fCalorimeter.Data());
953 char namePt[buffersize];
954 snprintf(namePt,buffersize,"AliAnaPi0_%s_cPt",fCalorimeter.Data());
955 TCanvas * cPt = new TCanvas(namePt, "", 400, 10, 600, 700) ;
960 TH1D * hPt = (TH1D*) fhRe1[0]->Project3D("x");
961 hPt->SetLineColor(2);
962 hPt->SetTitle("No cut on M_{#gamma#gamma} ");
966 TH1D * hPtIM1 = (TH1D*)fhRe1[0]->ProjectionX(Form("Pt1_%s",fCalorimeter.Data()), fhRe1[0]->GetZaxis()->FindBin(0.05),fhRe1[0]->GetZaxis()->FindBin(0.21),0, -1,"");
967 // TH3F * hRe1IM1 = (TH3F*)fhRe1[0]->Clone(Form("Pt1_%s",fCalorimeter.Data()));
968 // hRe1IM1->GetZaxis()->SetRangeUser(0.05,0.21);
969 // TH1D * hPtIM1 = (TH1D*) hRe1IM1->Project3D("x");
970 hPtIM1->SetLineColor(2);
971 hPtIM1->SetTitle("0.05 < M_{#gamma#gamma} < 0.21 GeV/c^{2}");
975 TH1D * hPtIM2 = (TH1D*)fhRe1[0]->ProjectionX(Form("Pt2_%s",fCalorimeter.Data()), fhRe1[0]->GetZaxis()->FindBin(0.09),fhRe1[0]->GetZaxis()->FindBin(0.17),0, -1,"");
976 // TH3F * hRe1IM2 = (TH3F*)fhRe1[0]->Clone(Form("Pt2_%s",fCalorimeter.Data()));
977 // hRe1IM2->GetZaxis()->SetRangeUser(0.09,0.17);
978 // TH1D * hPtIM2 = (TH1D*) hRe1IM2->Project3D("x");
979 hPtIM2->SetLineColor(2);
980 hPtIM2->SetTitle("0.09 < M_{#gamma#gamma} < 0.17 GeV/c^{2}");
984 TH1D * hPtIM3 = (TH1D*)fhRe1[0]->ProjectionX(Form("Pt3_%s",fCalorimeter.Data()), fhRe1[0]->GetZaxis()->FindBin(0.11),fhRe1[0]->GetZaxis()->FindBin(0.15),0, -1,"");
985 // TH3F * hRe1IM3 = (TH3F*)fhRe1[0]->Clone(Form("Pt3_%s",fCalorimeter.Data()));
986 // hRe1IM3->GetZaxis()->SetRangeUser(0.11,0.15);
987 // TH1D * hPtIM3 = (TH1D*) hRe1IM1->Project3D("x");
988 hPtIM3->SetLineColor(2);
989 hPtIM3->SetTitle("0.11 < M_{#gamma#gamma} < 0.15 GeV/c^{2}");
992 char namePtF[buffersize];
993 snprintf(namePtF,buffersize,"AliAnaPi0_%s_Pt.eps",fCalorimeter.Data());
996 char line[buffersize] ;
997 snprintf(line,buffersize,".!tar -zcf %s_%s.tar.gz *.eps", GetName(),fCalorimeter.Data()) ;
998 gROOT->ProcessLine(line);
999 snprintf(line, buffersize,".!rm -fR AliAnaPi0_%s*.eps",fCalorimeter.Data());
1000 gROOT->ProcessLine(line);
1002 printf(" AliAnaPi0::Terminate() - !! All the eps files are in %s_%s.tar.gz !!!\n", GetName(), fCalorimeter.Data());
1005 //____________________________________________________________________________________________________________________________________________________
1006 Int_t AliAnaPi0::GetEventIndex(AliAODPWG4Particle * part, Double_t * vert)
1008 // retieves the event index and checks the vertex
1009 // in the mixed buffer returns -2 if vertex NOK
1010 // for normal events returns 0 if vertex OK and -1 if vertex NOK
1012 Int_t evtIndex = -1 ;
1013 if(GetReader()->GetDataType()!=AliCaloTrackReader::kMC){
1015 if (GetMixedEvent()){
1017 evtIndex = GetMixedEvent()->EventIndexForCaloCluster(part->GetCaloLabel(0)) ;
1018 GetVertex(vert,evtIndex);
1020 if(TMath::Abs(vert[2])> GetZvertexCut())
1021 evtIndex = -2 ; //Event can not be used (vertex, centrality,... cuts not fulfilled)
1022 } else {// Single event
1026 if(TMath::Abs(vert[2])> GetZvertexCut())
1027 evtIndex = -1 ; //Event can not be used (vertex, centrality,... cuts not fulfilled)