]> 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
Implement different event selection based on goodness of the reconstructed vertex...
[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"
50f39b97 29#include "TH2D.h"
1c5acb87 30//#include "Riostream.h"
6639984f 31#include "TCanvas.h"
32#include "TPad.h"
33#include "TROOT.h"
477d6cee 34#include "TClonesArray.h"
0c1383b5 35#include "TObjString.h"
1c5acb87 36
37//---- AliRoot system ----
38#include "AliAnaPi0.h"
39#include "AliCaloTrackReader.h"
40#include "AliCaloPID.h"
6639984f 41#include "AliStack.h"
ff45398a 42#include "AliFiducialCut.h"
477d6cee 43#include "TParticle.h"
477d6cee 44#include "AliVEvent.h"
6921fa00 45#include "AliESDCaloCluster.h"
46#include "AliESDEvent.h"
47#include "AliAODEvent.h"
50f39b97 48#include "AliNeutralMesonSelection.h"
c8fe2783 49#include "AliMixedEvent.h"
50
591cc579 51
1c5acb87 52ClassImp(AliAnaPi0)
53
54//________________________________________________________________________________________________________________________________________________
55AliAnaPi0::AliAnaPi0() : AliAnaPartCorrBaseClass(),
5025c139 56fDoOwnMix(kFALSE),fNCentrBin(0),//fNZvertBin(0),fNrpBin(0),
af7b3903 57fNmaxMixEv(0), fCalorimeter(""),
58fNModules(12), fUseAngleCut(kFALSE), fEventsList(0x0), fMultiCutAna(kFALSE),
398c93cc 59fNPtCuts(0),fNAsymCuts(0), fNCellNCuts(0),fNPIDBits(0), fMakeInvPtPlots(kFALSE), fSameSM(kFALSE),
d7c10d78 60fhReMod(0x0),fhReDiffMod(0x0),
5ae09196 61fhRe1(0x0), fhMi1(0x0), fhRe2(0x0), fhMi2(0x0), fhRe3(0x0), fhMi3(0x0),
62fhReInvPt1(0x0), fhMiInvPt1(0x0), fhReInvPt2(0x0), fhMiInvPt2(0x0), fhReInvPt3(0x0), fhMiInvPt3(0x0),
af7b3903 63fhRePtNCellAsymCuts(0x0), fhRePIDBits(0x0),fhRePtMult(0x0), fhRePtAsym(0x0), fhRePtAsymPi0(0x0),fhRePtAsymEta(0x0),
5ae09196 64fhEvents(0x0), fhRealOpeningAngle(0x0),fhRealCosOpeningAngle(0x0),
50f39b97 65fhPrimPt(0x0), fhPrimAccPt(0x0), fhPrimY(0x0), fhPrimAccY(0x0), fhPrimPhi(0x0), fhPrimAccPhi(0x0),
66fhPrimOpeningAngle(0x0),fhPrimCosOpeningAngle(0x0)
1c5acb87 67{
68//Default Ctor
69 InitParameters();
70
71}
1c5acb87 72
73//________________________________________________________________________________________________________________________________________________
74AliAnaPi0::~AliAnaPi0() {
477d6cee 75 // Remove event containers
7e7694bb 76
77 if(fDoOwnMix && fEventsList){
477d6cee 78 for(Int_t ic=0; ic<fNCentrBin; ic++){
5025c139 79 for(Int_t iz=0; iz<GetNZvertBin(); iz++){
80 for(Int_t irp=0; irp<GetNRPBin(); irp++){
81 fEventsList[ic*GetNZvertBin()*GetNRPBin()+iz*GetNRPBin()+irp]->Delete() ;
82 delete fEventsList[ic*GetNZvertBin()*GetNRPBin()+iz*GetNRPBin()+irp] ;
7e7694bb 83 }
477d6cee 84 }
85 }
86 delete[] fEventsList;
87 fEventsList=0 ;
88 }
591cc579 89
1c5acb87 90}
91
92//________________________________________________________________________________________________________________________________________________
93void AliAnaPi0::InitParameters()
94{
95//Init parameters when first called the analysis
96//Set default parameters
a3aebfff 97 SetInputAODName("PWG4Particle");
98
99 AddToHistogramsName("AnaPi0_");
6921fa00 100 fNModules = 12; // set maximum to maximum number of EMCAL modules
477d6cee 101 fNCentrBin = 1;
5025c139 102// fNZvertBin = 1;
103// fNrpBin = 1;
477d6cee 104 fNmaxMixEv = 10;
5025c139 105
477d6cee 106 fCalorimeter = "PHOS";
50f39b97 107 fUseAngleCut = kFALSE;
5ae09196 108
109 fMultiCutAna = kFALSE;
110
111 fNPtCuts = 3;
af7b3903 112 fPtCuts[0] = 0.; fPtCuts[1] = 0.3; fPtCuts[2] = 0.5;
113 for(Int_t i = fNPtCuts; i < 10; i++)fPtCuts[i] = 0.;
5ae09196 114
9c59b5fe 115 fNAsymCuts = 4;
116 fAsymCuts[0] = 1.; fAsymCuts[1] = 0.8; fAsymCuts[2] = 0.6; fAsymCuts[3] = 0.1;
af7b3903 117 for(Int_t i = fNAsymCuts; i < 10; i++)fAsymCuts[i] = 0.;
118
5ae09196 119 fNCellNCuts = 3;
af7b3903 120 fCellNCuts[0] = 0; fCellNCuts[1] = 1; fCellNCuts[2] = 2;
021c573a 121 for(Int_t i = fNCellNCuts; i < 10; i++)fCellNCuts[i] = 0;
af7b3903 122
123 fNPIDBits = 2;
124 fPIDBits[0] = 0; fPIDBits[1] = 2; // fPIDBits[2] = 4; fPIDBits[3] = 6;// check, no cut, dispersion, neutral, dispersion&&neutral
021c573a 125 for(Int_t i = fNPIDBits; i < 10; i++)fPIDBits[i] = 0;
af7b3903 126
1c5acb87 127}
1c5acb87 128
0c1383b5 129
130//________________________________________________________________________________________________________________________________________________
131TObjString * AliAnaPi0::GetAnalysisCuts()
132{
af7b3903 133 //Save parameters used for analysis
134 TString parList ; //this will be list of parameters used for this analysis.
135 const Int_t buffersize = 255;
136 char onePar[buffersize] ;
137 snprintf(onePar,buffersize,"--- AliAnaPi0 ---\n") ;
138 parList+=onePar ;
139 snprintf(onePar,buffersize,"Number of bins in Centrality: %d \n",fNCentrBin) ;
140 parList+=onePar ;
141 snprintf(onePar,buffersize,"Number of bins in Z vert. pos: %d \n",GetNZvertBin()) ;
142 parList+=onePar ;
143 snprintf(onePar,buffersize,"Number of bins in Reac. Plain: %d \n",GetNRPBin()) ;
144 parList+=onePar ;
145 snprintf(onePar,buffersize,"Depth of event buffer: %d \n",fNmaxMixEv) ;
146 parList+=onePar ;
398c93cc 147 snprintf(onePar,buffersize,"Pair in same Module: %d; Make InvPt plots %d \n",fSameSM, fMakeInvPtPlots) ;
af7b3903 148 parList+=onePar ;
149 snprintf(onePar,buffersize," Asymmetry cuts: n = %d, asymmetry < ",fNAsymCuts) ;
150 for(Int_t i = 0; i < fNAsymCuts; i++) snprintf(onePar,buffersize,"%s %2.2f;",onePar,fAsymCuts[i]);
151 parList+=onePar ;
152 snprintf(onePar,buffersize," PID selection bits: n = %d, PID bit =\n",fNPIDBits) ;
153 for(Int_t i = 0; i < fNPIDBits; i++) snprintf(onePar,buffersize,"%s %d;",onePar,fPIDBits[i]);
154 parList+=onePar ;
155 snprintf(onePar,buffersize,"Cuts: \n") ;
156 parList+=onePar ;
157 snprintf(onePar,buffersize,"Z vertex position: -%f < z < %f \n",GetZvertexCut(),GetZvertexCut()) ;
158 parList+=onePar ;
159 snprintf(onePar,buffersize,"Calorimeter: %s \n",fCalorimeter.Data()) ;
160 parList+=onePar ;
161 snprintf(onePar,buffersize,"Number of modules: %d \n",fNModules) ;
162 parList+=onePar ;
db2bf6fd 163 if(fMultiCutAna){
164 snprintf(onePar, buffersize," pT cuts: n = %d, pt > ",fNPtCuts) ;
165 for(Int_t i = 0; i < fNPtCuts; i++) snprintf(onePar,buffersize,"%s %2.2f;",onePar,fPtCuts[i]);
166 parList+=onePar ;
db2bf6fd 167 snprintf(onePar,buffersize, " N cell in cluster cuts: n = %d, nCell > ",fNCellNCuts) ;
168 for(Int_t i = 0; i < fNCellNCuts; i++) snprintf(onePar,buffersize,"%s %d;",onePar,fCellNCuts[i]);
169 parList+=onePar ;
db2bf6fd 170 }
171
af7b3903 172 return new TObjString(parList) ;
0c1383b5 173}
174
2e557d1c 175//________________________________________________________________________________________________________________________________________________
176TList * AliAnaPi0::GetCreateOutputObjects()
177{
477d6cee 178 // Create histograms to be saved in output file and
179 // store them in fOutputContainer
180
181 //create event containers
5025c139 182 fEventsList = new TList*[fNCentrBin*GetNZvertBin()*GetNRPBin()] ;
1c5acb87 183
477d6cee 184 for(Int_t ic=0; ic<fNCentrBin; ic++){
5025c139 185 for(Int_t iz=0; iz<GetNZvertBin(); iz++){
186 for(Int_t irp=0; irp<GetNRPBin(); irp++){
187 fEventsList[ic*GetNZvertBin()*GetNRPBin()+iz*GetNRPBin()+irp] = new TList() ;
af7b3903 188 fEventsList[ic*GetNZvertBin()*GetNRPBin()+iz*GetNRPBin()+irp]->SetOwner(kFALSE);
477d6cee 189 }
190 }
191 }
7e7694bb 192
477d6cee 193 TList * outputContainer = new TList() ;
194 outputContainer->SetName(GetName());
6921fa00 195
af7b3903 196 fhReMod = new TH2D*[fNModules] ;
197 fhReDiffMod = new TH2D*[fNModules+1] ;
198
199 fhRe1 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
200 fhRe2 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
201 fhRe3 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
202 fhMi1 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
203 fhMi2 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
204 fhMi3 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
398c93cc 205 if(fMakeInvPtPlots ) {
206 fhReInvPt1 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
207 fhReInvPt2 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
208 fhReInvPt3 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
209 fhMiInvPt1 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
210 fhMiInvPt2 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
211 fhMiInvPt3 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
212 }
5ae09196 213 const Int_t buffersize = 255;
214 char key[buffersize] ;
215 char title[buffersize] ;
477d6cee 216
5a2dbc3c 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();
226
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();
af7b3903 233 Int_t ntrmbins = GetHistoTrackMultiplicityBins();
234 Int_t ntrmmax = GetHistoTrackMultiplicityMax();
235 Int_t ntrmmin = GetHistoTrackMultiplicityMin();
236
477d6cee 237 for(Int_t ic=0; ic<fNCentrBin; ic++){
af7b3903 238 for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
239 for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
240 Int_t index = ((ic*fNPIDBits)+ipid)*fNAsymCuts + iasym;
241 //printf("cen %d, pid %d, asy %d, Index %d\n",ic,ipid,iasym,index);
7e7694bb 242 //Distance to bad module 1
af7b3903 243 snprintf(key, buffersize,"hRe_cen%d_pidbit%d_asy%d_dist1",ic,ipid,iasym) ;
244 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 1",
245 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
246 fhRe1[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
247 fhRe1[index]->SetXTitle("p_{T} (GeV/c)");
248 fhRe1[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
249 //printf("name: %s\n ",fhRe1[index]->GetName());
250 outputContainer->Add(fhRe1[index]) ;
7e7694bb 251
252 //Distance to bad module 2
af7b3903 253 snprintf(key, buffersize,"hRe_cen%d_pidbit%d_asy%d_dist2",ic,ipid,iasym) ;
254 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 2",
255 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
256 fhRe2[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
257 fhRe2[index]->SetXTitle("p_{T} (GeV/c)");
258 fhRe2[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
259 outputContainer->Add(fhRe2[index]) ;
7e7694bb 260
261 //Distance to bad module 3
af7b3903 262 snprintf(key, buffersize,"hRe_cen%d_pidbit%d_asy%d_dist3",ic,ipid,iasym) ;
263 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 3",
264 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
265 fhRe3[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
266 fhRe3[index]->SetXTitle("p_{T} (GeV/c)");
267 fhRe3[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
268 outputContainer->Add(fhRe3[index]) ;
5ae09196 269
398c93cc 270 if(fMakeInvPtPlots ) {
271 //Inverse pT
272 //Distance to bad module 1
273 snprintf(key, buffersize,"hReInvPt_cen%d_pidbit%d_asy%d_dist1",ic,ipid,iasym) ;
274 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 1",
275 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
276 fhReInvPt1[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
277 fhReInvPt1[index]->SetXTitle("p_{T} (GeV/c)");
278 fhReInvPt1[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
279 outputContainer->Add(fhReInvPt1[index]) ;
280
281 //Distance to bad module 2
282 snprintf(key, buffersize,"hReInvPt_cen%d_pidbit%d_asy%d_dist2",ic,ipid,iasym) ;
283 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 2",
284 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
285 fhReInvPt2[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
286 fhReInvPt2[index]->SetXTitle("p_{T} (GeV/c)");
287 fhReInvPt2[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
288 outputContainer->Add(fhReInvPt2[index]) ;
289
290 //Distance to bad module 3
291 snprintf(key, buffersize,"hReInvPt_cen%d_pidbit%d_asy%d_dist3",ic,ipid,iasym) ;
292 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 3",
293 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
294 fhReInvPt3[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
295 fhReInvPt3[index]->SetXTitle("p_{T} (GeV/c)");
296 fhReInvPt3[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
297 outputContainer->Add(fhReInvPt3[index]) ;
298 }
5ae09196 299
af7b3903 300 if(fDoOwnMix){
301 //Distance to bad module 1
302 snprintf(key, buffersize,"hMi_cen%d_pidbit%d_asy%d_dist1",ic,ipid,iasym) ;
303 snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 1",
304 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
305 fhMi1[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
306 fhMi1[index]->SetXTitle("p_{T} (GeV/c)");
307 fhMi1[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
308 outputContainer->Add(fhMi1[index]) ;
309
310 //Distance to bad module 2
311 snprintf(key, buffersize,"hMi_cen%d_pidbit%d_asy%d_dist2",ic,ipid,iasym) ;
312 snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 2",
313 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
314 fhMi2[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
315 fhMi2[index]->SetXTitle("p_{T} (GeV/c)");
316 fhMi2[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
317 outputContainer->Add(fhMi2[index]) ;
318
319 //Distance to bad module 3
320 snprintf(key, buffersize,"hMi_cen%d_pidbit%d_asy%d_dist3",ic,ipid,iasym) ;
321 snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 3",
322 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
323 fhMi3[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
324 fhMi3[index]->SetXTitle("p_{T} (GeV/c)");
325 fhMi3[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
326 outputContainer->Add(fhMi3[index]) ;
398c93cc 327 if(fMakeInvPtPlots ) {
328 //Inverse pT
329 //Distance to bad module 1
330 snprintf(key, buffersize,"hMiInvPt_cen%d_pidbit%d_asy%d_dist1",ic,ipid,iasym) ;
331 snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 1",
332 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
333 fhMiInvPt1[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
334 fhMiInvPt1[index]->SetXTitle("p_{T} (GeV/c)");
335 fhMiInvPt1[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
336 outputContainer->Add(fhMiInvPt1[index]) ;
337
338 //Distance to bad module 2
339 snprintf(key, buffersize,"hMiInvPt_cen%d_pidbit%d_asy%d_dist2",ic,ipid,iasym) ;
340 snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 2",
341 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
342 fhMiInvPt2[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
343 fhMiInvPt2[index]->SetXTitle("p_{T} (GeV/c)");
344 fhMiInvPt2[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
345 outputContainer->Add(fhMiInvPt2[index]) ;
346
347 //Distance to bad module 3
348 snprintf(key, buffersize,"hMiInvPt_cen%d_pidbit%d_asy%d_dist3",ic,ipid,iasym) ;
349 snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f,dist bad 3",
350 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
351 fhMiInvPt3[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
352 fhMiInvPt3[index]->SetXTitle("p_{T} (GeV/c)");
353 fhMiInvPt3[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
354 outputContainer->Add(fhMiInvPt3[index]) ;
355 }
af7b3903 356 }
7e7694bb 357 }
1c5acb87 358 }
477d6cee 359 }
360
9c59b5fe 361 fhRePtAsym = new TH2D("hRePtAsym","Asymmetry vs pt, for pairs",nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax) ;
af7b3903 362 fhRePtAsym->SetXTitle("p_{T} (GeV/c)");
363 fhRePtAsym->SetYTitle("Asymmetry");
364 outputContainer->Add(fhRePtAsym);
365
9c59b5fe 366 fhRePtAsymPi0 = new TH2D("hRePtAsymPi0","Asymmetry vs pt, for pairs close to #pi^{0} mass",nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax) ;
af7b3903 367 fhRePtAsymPi0->SetXTitle("p_{T} (GeV/c)");
368 fhRePtAsymPi0->SetYTitle("Asymmetry");
369 outputContainer->Add(fhRePtAsymPi0);
370
9c59b5fe 371 fhRePtAsymEta = new TH2D("hRePtAsymEta","Asymmetry vs pt, for pairs close to #eta mass",nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax) ;
af7b3903 372 fhRePtAsymEta->SetXTitle("p_{T} (GeV/c)");
373 fhRePtAsymEta->SetYTitle("Asymmetry");
374 outputContainer->Add(fhRePtAsymEta);
375
5ae09196 376 if(fMultiCutAna){
377
378 fhRePIDBits = new TH2D*[fNPIDBits];
379 for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
380 snprintf(key, buffersize,"hRe_pidbit%d",ipid) ;
381 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for PIDBit=%d",fPIDBits[ipid]) ;
382 fhRePIDBits[ipid] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
af7b3903 383 fhRePIDBits[ipid]->SetXTitle("p_{T} (GeV/c)");
384 fhRePIDBits[ipid]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
5ae09196 385 outputContainer->Add(fhRePIDBits[ipid]) ;
386 }// pid bit loop
387
388 fhRePtNCellAsymCuts = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
389 for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
390 for(Int_t icell=0; icell<fNCellNCuts; icell++){
391 for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
392 snprintf(key, buffersize,"hRe_pt%d_cell%d_asym%d",ipt,icell,iasym) ;
af7b3903 393 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for pt >%2.2f, ncell>%d and asym >%1.2f ",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]) ;
5ae09196 394 Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
395 //printf("ipt %d, icell %d, iassym %d, index %d\n",ipt, icell, iasym, index);
396 fhRePtNCellAsymCuts[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
af7b3903 397 fhRePtNCellAsymCuts[index]->SetXTitle("p_{T} (GeV/c)");
398 fhRePtNCellAsymCuts[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
5ae09196 399 outputContainer->Add(fhRePtNCellAsymCuts[index]) ;
400 }
401 }
402 }
821c8090 403
af7b3903 404 fhRePtMult = new TH3D*[fNAsymCuts] ;
405 for(Int_t iasym = 0; iasym<fNAsymCuts; iasym++){
406 fhRePtMult[iasym] = new TH3D(Form("hRePtMult_asym%d",iasym),Form("(p_{T},C,M)_{#gamma#gamma}, A<%1.2f",fAsymCuts[iasym]),
407 nptbins,ptmin,ptmax,ntrmbins,ntrmmin,ntrmmax,nmassbins,massmin,massmax);
408 fhRePtMult[iasym]->SetXTitle("p_{T} (GeV/c)");
409 fhRePtMult[iasym]->SetYTitle("Track multiplicity");
410 fhRePtMult[iasym]->SetZTitle("m_{#gamma,#gamma} (GeV/c^{2})");
411 outputContainer->Add(fhRePtMult[iasym]) ;
412 }
413
5ae09196 414 }// multi cuts analysis
415
477d6cee 416 fhEvents=new TH3D("hEvents","Number of events",fNCentrBin,0.,1.*fNCentrBin,
5025c139 417 GetNZvertBin(),0.,1.*GetNZvertBin(),GetNRPBin(),0.,1.*GetNRPBin()) ;
477d6cee 418 outputContainer->Add(fhEvents) ;
50f39b97 419
420 fhRealOpeningAngle = new TH2D
421 ("hRealOpeningAngle","Angle between all #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,200,0,0.5);
422 fhRealOpeningAngle->SetYTitle("#theta(rad)");
423 fhRealOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
424 outputContainer->Add(fhRealOpeningAngle) ;
7e7694bb 425
50f39b97 426 fhRealCosOpeningAngle = new TH2D
427 ("hRealCosOpeningAngle","Cosinus of angle between all #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,200,-1,1);
428 fhRealCosOpeningAngle->SetYTitle("cos (#theta) ");
429 fhRealCosOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
430 outputContainer->Add(fhRealCosOpeningAngle) ;
431
477d6cee 432 //Histograms filled only if MC data is requested
0ae57829 433 if(IsDataMC()){
5ae09196 434
7e7694bb 435 fhPrimPt = new TH1D("hPrimPt","Primary pi0 pt",nptbins,ptmin,ptmax) ;
5a2dbc3c 436 fhPrimAccPt = new TH1D("hPrimAccPt","Primary pi0 pt with both photons in acceptance",nptbins,ptmin,ptmax) ;
477d6cee 437 outputContainer->Add(fhPrimPt) ;
438 outputContainer->Add(fhPrimAccPt) ;
439
5a2dbc3c 440 fhPrimY = new TH1D("hPrimaryRapidity","Rapidity of primary pi0",netabins,etamin,etamax) ;
477d6cee 441 outputContainer->Add(fhPrimY) ;
442
5a2dbc3c 443 fhPrimAccY = new TH1D("hPrimAccRapidity","Rapidity of primary pi0",netabins,etamin,etamax) ;
477d6cee 444 outputContainer->Add(fhPrimAccY) ;
445
7e7694bb 446 fhPrimPhi = new TH1D("hPrimaryPhi","Azimithal of primary pi0",nphibins,phimin*TMath::RadToDeg(),phimax*TMath::RadToDeg()) ;
477d6cee 447 outputContainer->Add(fhPrimPhi) ;
448
5a2dbc3c 449 fhPrimAccPhi = new TH1D("hPrimAccPhi","Azimithal of primary pi0 with accepted daughters",nphibins,phimin*TMath::RadToDeg(),phimax*TMath::RadToDeg()) ;
477d6cee 450 outputContainer->Add(fhPrimAccPhi) ;
50f39b97 451
452
453 fhPrimOpeningAngle = new TH2D
7e7694bb 454 ("hPrimOpeningAngle","Angle between all primary #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,100,0,0.5);
50f39b97 455 fhPrimOpeningAngle->SetYTitle("#theta(rad)");
456 fhPrimOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
457 outputContainer->Add(fhPrimOpeningAngle) ;
458
459 fhPrimCosOpeningAngle = new TH2D
7e7694bb 460 ("hPrimCosOpeningAngle","Cosinus of angle between all primary #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,100,-1,1);
50f39b97 461 fhPrimCosOpeningAngle->SetYTitle("cos (#theta) ");
462 fhPrimCosOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
463 outputContainer->Add(fhPrimCosOpeningAngle) ;
464
477d6cee 465 }
50f39b97 466
821c8090 467 TString * pairname = new TString[fNModules];
468 if(fCalorimeter=="EMCAL"){
469 pairname[0]="A side (0-2)";
470 pairname[1]="C side (1-3)";
471 pairname[2]="Sector 0 (0-1)";
472 pairname[3]="Sector 1 (2-3)";
473 for(Int_t i = 4 ; i < fNModules ; i++) pairname[i]="";}
474 if(fCalorimeter=="PHOS") {
475 pairname[0]="(0-1)";
476 pairname[1]="(0-2)";
af7b3903 477 pairname[2]="(1-2)";
821c8090 478 for(Int_t i = 3 ; i < fNModules ; i++) pairname[i]="";}
479
6921fa00 480 for(Int_t imod=0; imod<fNModules; imod++){
50f39b97 481 //Module dependent invariant mass
5ae09196 482 snprintf(key, buffersize,"hReMod_%d",imod) ;
483 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for Module %d",imod) ;
af7b3903 484 fhReMod[imod] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
485 fhReMod[imod]->SetXTitle("p_{T} (GeV/c)");
486 fhReMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
50f39b97 487 outputContainer->Add(fhReMod[imod]) ;
821c8090 488
489 snprintf(key, buffersize,"hReDiffMod_%d",imod) ;
490 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for Different Modules: %s",(pairname[imod]).Data()) ;
af7b3903 491 fhReDiffMod[imod] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
492 fhReDiffMod[imod]->SetXTitle("p_{T} (GeV/c)");
493 fhReDiffMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
821c8090 494 outputContainer->Add(fhReDiffMod[imod]) ;
6921fa00 495 }
50f39b97 496
821c8090 497 delete [] pairname;
498
499 snprintf(key, buffersize,"hReDiffMod_%d",fNModules) ;
54769bc0 500 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for all Modules Combination") ;
af7b3903 501 fhReDiffMod[fNModules] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
821c8090 502 outputContainer->Add(fhReDiffMod[fNModules]) ;
503
504
eee5fcf1 505// for(Int_t i = 0; i < outputContainer->GetEntries() ; i++){
506//
507// printf("Histogram %d, name: %s\n ",i, outputContainer->At(i)->GetName());
508//
509// }
510
477d6cee 511 return outputContainer;
1c5acb87 512}
513
514//_________________________________________________________________________________________________________________________________________________
515void AliAnaPi0::Print(const Option_t * /*opt*/) const
516{
477d6cee 517 //Print some relevant parameters set for the analysis
a3aebfff 518 printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
477d6cee 519 AliAnaPartCorrBaseClass::Print(" ");
a3aebfff 520
477d6cee 521 printf("Number of bins in Centrality: %d \n",fNCentrBin) ;
5025c139 522 printf("Number of bins in Z vert. pos: %d \n",GetNZvertBin()) ;
523 printf("Number of bins in Reac. Plain: %d \n",GetNRPBin()) ;
477d6cee 524 printf("Depth of event buffer: %d \n",fNmaxMixEv) ;
af7b3903 525 printf("Pair in same Module: %d \n",fSameSM) ;
477d6cee 526 printf("Cuts: \n") ;
5025c139 527 printf("Z vertex position: -%2.3f < z < %2.3f \n",GetZvertexCut(),GetZvertexCut()) ;
50f39b97 528 printf("Number of modules: %d \n",fNModules) ;
529 printf("Select pairs with their angle: %d \n",fUseAngleCut) ;
af7b3903 530 printf("Asymmetry cuts: n = %d, \n",fNAsymCuts) ;
531 printf("\tasymmetry < ");
532 for(Int_t i = 0; i < fNAsymCuts; i++) printf("%2.2f ",fAsymCuts[i]);
533 printf("\n");
534
535 printf("PID selection bits: n = %d, \n",fNPIDBits) ;
536 printf("\tPID bit = ");
537 for(Int_t i = 0; i < fNPIDBits; i++) printf("%d ",fPIDBits[i]);
538 printf("\n");
539
db2bf6fd 540 if(fMultiCutAna){
541 printf("pT cuts: n = %d, \n",fNPtCuts) ;
542 printf("\tpT > ");
543 for(Int_t i = 0; i < fNPtCuts; i++) printf("%2.2f ",fPtCuts[i]);
544 printf("GeV/c\n");
545
546 printf("N cell in cluster cuts: n = %d, \n",fNCellNCuts) ;
547 printf("\tnCell > ");
548 for(Int_t i = 0; i < fNCellNCuts; i++) printf("%d ",fCellNCuts[i]);
549 printf("\n");
550
db2bf6fd 551 }
477d6cee 552 printf("------------------------------------------------------\n") ;
1c5acb87 553}
554
5ae09196 555//_____________________________________________________________
556void AliAnaPi0::FillAcceptanceHistograms(){
557 //Fill acceptance histograms if MC data is available
c8fe2783 558
5ae09196 559 if(IsDataMC() && GetReader()->ReadStack()){
560 AliStack * stack = GetMCStack();
561 if(stack && (IsDataMC() || (GetReader()->GetDataType() == AliCaloTrackReader::kMC)) ){
562 for(Int_t i=0 ; i<stack->GetNprimary(); i++){
563 TParticle * prim = stack->Particle(i) ;
564 if(prim->GetPdgCode() == 111){
565 Double_t pi0Pt = prim->Pt() ;
566 //printf("pi0, pt %2.2f\n",pi0Pt);
567 if(prim->Energy() == TMath::Abs(prim->Pz())) continue ; //Protection against floating point exception
568 Double_t pi0Y = 0.5*TMath::Log((prim->Energy()-prim->Pz())/(prim->Energy()+prim->Pz())) ;
569 Double_t phi = TMath::RadToDeg()*prim->Phi() ;
570 if(TMath::Abs(pi0Y) < 0.5){
571 fhPrimPt->Fill(pi0Pt) ;
572 }
573 fhPrimY ->Fill(pi0Y) ;
574 fhPrimPhi->Fill(phi) ;
575
576 //Check if both photons hit Calorimeter
577 Int_t iphot1=prim->GetFirstDaughter() ;
578 Int_t iphot2=prim->GetLastDaughter() ;
579 if(iphot1>-1 && iphot1<stack->GetNtrack() && iphot2>-1 && iphot2<stack->GetNtrack()){
580 TParticle * phot1 = stack->Particle(iphot1) ;
581 TParticle * phot2 = stack->Particle(iphot2) ;
582 if(phot1 && phot2 && phot1->GetPdgCode()==22 && phot2->GetPdgCode()==22){
583 //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",
584 // phot1->Pt(), phot1->Phi()*180./3.1415, phot1->Eta(), phot2->Pt(), phot2->Phi()*180./3.1415, phot2->Eta());
585
586 TLorentzVector lv1, lv2;
587 phot1->Momentum(lv1);
588 phot2->Momentum(lv2);
589
590 Bool_t inacceptance = kFALSE;
591 if(fCalorimeter == "PHOS"){
592 if(GetPHOSGeometry() && GetCaloUtils()->IsPHOSGeoMatrixSet()){
593 Int_t mod ;
594 Double_t x,z ;
595 if(GetPHOSGeometry()->ImpactOnEmc(phot1,mod,z,x) && GetPHOSGeometry()->ImpactOnEmc(phot2,mod,z,x))
596 inacceptance = kTRUE;
597 if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
598 }
599 else{
600
601 if(GetFiducialCut()->IsInFiducialCut(lv1,fCalorimeter) && GetFiducialCut()->IsInFiducialCut(lv2,fCalorimeter))
602 inacceptance = kTRUE ;
603 if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
604 }
605
606 }
607 else if(fCalorimeter == "EMCAL" && GetCaloUtils()->IsEMCALGeoMatrixSet()){
608 if(GetEMCALGeometry()){
609 if(GetEMCALGeometry()->Impact(phot1) && GetEMCALGeometry()->Impact(phot2))
610 inacceptance = kTRUE;
611 if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
612 }
613 else{
614 if(GetFiducialCut()->IsInFiducialCut(lv1,fCalorimeter) && GetFiducialCut()->IsInFiducialCut(lv2,fCalorimeter))
615 inacceptance = kTRUE ;
616 if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
617 }
618 }
619
620 if(inacceptance){
621
622 fhPrimAccPt->Fill(pi0Pt) ;
623 fhPrimAccPhi->Fill(phi) ;
624 fhPrimAccY->Fill(pi0Y) ;
625 Double_t angle = lv1.Angle(lv2.Vect());
626 fhPrimOpeningAngle ->Fill(pi0Pt,angle);
627 fhPrimCosOpeningAngle->Fill(pi0Pt,TMath::Cos(angle));
628
629 }//Accepted
630 }// 2 photons
631 }//Check daughters exist
632 }// Primary pi0
633 }//loop on primaries
634 }//stack exists and data is MC
635 }//read stack
636 else if(GetReader()->ReadAODMCParticles()){
637 if(GetDebug() >= 0) printf("AliAnaPi0::FillAcceptanceHistograms() - Acceptance calculation with MCParticles not implemented yet\n");
638 }
639}
640
1c5acb87 641//____________________________________________________________________________________________________________________________________________________
6639984f 642void AliAnaPi0::MakeAnalysisFillHistograms()
1c5acb87 643{
477d6cee 644 //Process one event and extract photons from AOD branch
645 // filled with AliAnaPhoton and fill histos with invariant mass
646
5ae09196 647 //In case of MC data, fill acceptance histograms
648 FillAcceptanceHistograms();
649
477d6cee 650 //Apply some cuts on event: vertex position and centrality range
651 Int_t iRun=(GetReader()->GetInputEvent())->GetRunNumber() ;
652 if(IsBadRun(iRun)) return ;
653
477d6cee 654 Int_t nPhot = GetInputAODBranch()->GetEntriesFast() ;
c8fe2783 655 if(GetDebug() > 1)
656 printf("AliAnaPi0::MakeAnalysisFillHistograms() - Photon entries %d\n", nPhot);
657 if(nPhot < 2 )
658 return ;
6921fa00 659 Int_t module1 = -1;
660 Int_t module2 = -1;
7e7694bb 661 Double_t vert[] = {0.0, 0.0, 0.0} ; //vertex
662 Int_t evtIndex1 = 0 ;
663 Int_t currentEvtIndex = -1 ;
664 Int_t curCentrBin = 0 ;
665 Int_t curRPBin = 0 ;
666 Int_t curZvertBin = 0 ;
667
477d6cee 668 for(Int_t i1=0; i1<nPhot-1; i1++){
669 AliAODPWG4Particle * p1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i1)) ;
7e7694bb 670 // get the event index in the mixed buffer where the photon comes from
671 // in case of mixing with analysis frame, not own mixing
c8fe2783 672 evtIndex1 = GetEventIndex(p1, vert) ;
5025c139 673 //printf("charge = %d\n", track->Charge());
c8fe2783 674 if ( evtIndex1 == -1 )
675 return ;
676 if ( evtIndex1 == -2 )
677 continue ;
2244659d 678 if(TMath::Abs(vert[2]) > GetZvertexCut()) continue ; //vertex cut
c8fe2783 679 if (evtIndex1 != currentEvtIndex) {
3ff76907 680 curCentrBin = GetEventCentrality();
e6cb049a 681 //In case too few centrality bins defined with respect the selected centrality class, put in the last bin the higher centrality classes
682 if(curCentrBin >= fNCentrBin) curCentrBin=fNCentrBin-1;
c8fe2783 683 curRPBin = 0 ;
5025c139 684 curZvertBin = (Int_t)(0.5*GetNZvertBin()*(vert[2]+GetZvertexCut())/GetZvertexCut()) ;
c8fe2783 685 fhEvents->Fill(curCentrBin+0.5,curZvertBin+0.5,curRPBin+0.5) ;
686 currentEvtIndex = evtIndex1 ;
ca468d44 687 if(GetDebug() > 1)
3ff76907 688 printf("AliAnaPi0::MakeAnalysisFillHistograms() - Centrality %d, Vertex Bin %d, RP bin %d\n",curCentrBin,curRPBin,curZvertBin);
c8fe2783 689 }
7e7694bb 690
f8006433 691 //printf("AliAnaPi0::MakeAnalysisFillHistograms(): Photon 1 Evt %d Vertex : %f,%f,%f\n",evtIndex1, GetVertex(evtIndex1)[0] ,GetVertex(evtIndex1)[1],GetVertex(evtIndex1)[2]);
af7b3903 692
477d6cee 693 TLorentzVector photon1(p1->Px(),p1->Py(),p1->Pz(),p1->E());
59b6bd99 694 //Get Module number
695 module1 = GetModuleNumber(p1);
477d6cee 696 for(Int_t i2=i1+1; i2<nPhot; i2++){
697 AliAODPWG4Particle * p2 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i2)) ;
c8fe2783 698 Int_t evtIndex2 = GetEventIndex(p2, vert) ;
699 if ( evtIndex2 == -1 )
700 return ;
701 if ( evtIndex2 == -2 )
702 continue ;
703 if (GetMixedEvent() && (evtIndex1 == evtIndex2))
7e7694bb 704 continue ;
f8006433 705 //printf("AliAnaPi0::MakeAnalysisFillHistograms(): Photon 2 Evt %d Vertex : %f,%f,%f\n",evtIndex2, GetVertex(evtIndex2)[0] ,GetVertex(evtIndex2)[1],GetVertex(evtIndex2)[2]);
477d6cee 706 TLorentzVector photon2(p2->Px(),p2->Py(),p2->Pz(),p2->E());
59b6bd99 707 //Get module number
708 module2 = GetModuleNumber(p2);
477d6cee 709 Double_t m = (photon1 + photon2).M() ;
710 Double_t pt = (photon1 + photon2).Pt();
711 Double_t a = TMath::Abs(p1->E()-p2->E())/(p1->E()+p2->E()) ;
712 if(GetDebug() > 2)
c8fe2783 713 printf("AliAnaPi0::MakeAnalysisFillHistograms() - Current Event: pT: photon1 %2.2f, photon2 %2.2f; Pair: pT %2.2f, mass %2.3f, a %f2.3\n",
714 p1->Pt(), p2->Pt(), pt,m,a);
50f39b97 715 //Check if opening angle is too large or too small compared to what is expected
716 Double_t angle = photon1.Angle(photon2.Vect());
717 //if(fUseAngleCut && !GetNeutralMesonSelection()->IsAngleInWindow((photon1+photon2).E(),angle)) continue;
718 //printf("angle %f\n",angle);
c8fe2783 719 if(fUseAngleCut && angle < 0.1)
720 continue;
af7b3903 721
722 //Fill module dependent histograms, put a cut on assymmetry on the first available cut in the array
723 if(a < fAsymCuts[0]){
724 if(module1==module2 && module1 >=0 && module1<fNModules)
725 fhReMod[module1]->Fill(pt,m) ;
726 else
727 fhReDiffMod[fNModules]->Fill(pt,m) ;
728
729 if(fCalorimeter=="EMCAL"){
730 if((module1==0 && module2==2) || (module1==2 && module2==0)) fhReDiffMod[0]->Fill(pt,m) ;
731 if((module1==1 && module2==3) || (module1==3 && module2==1)) fhReDiffMod[1]->Fill(pt,m) ;
732 if((module1==0 && module2==1) || (module1==1 && module2==0)) fhReDiffMod[2]->Fill(pt,m) ;
733 if((module1==2 && module2==3) || (module1==3 && module2==2)) fhReDiffMod[3]->Fill(pt,m) ;
734 }
735 else {
736 if((module1==0 && module2==1) || (module1==1 && module2==0)) fhReDiffMod[0]->Fill(pt,m) ;
737 if((module1==0 && module2==2) || (module1==2 && module2==0)) fhReDiffMod[1]->Fill(pt,m) ;
738 if((module1==1 && module2==2) || (module1==2 && module2==1)) fhReDiffMod[2]->Fill(pt,m) ;
739 }
821c8090 740 }
7e7694bb 741
af7b3903 742 //In case we want only pairs in same (super) module, check their origin.
743 Bool_t ok = kTRUE;
744 if(fSameSM && module1!=module2) ok=kFALSE;
745 if(ok){
746 //Fill histograms for different bad channel distance, centrality, assymmetry cut and pid bit
5ae09196 747 for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
af7b3903 748 if((p1->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton)) && (p2->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton))){
749 for(Int_t iasym=0; iasym < fNAsymCuts; iasym++){
750 if(a < fAsymCuts[iasym]){
751 Int_t index = ((curCentrBin*fNPIDBits)+ipid)*fNAsymCuts + iasym;
752 //printf("cen %d, pid %d, asy %d, Index %d\n",curCentrBin,ipid,iasym,index);
753 fhRe1 [index]->Fill(pt,m);
398c93cc 754 if(fMakeInvPtPlots)fhReInvPt1[index]->Fill(pt,m,1./pt) ;
af7b3903 755 if(p1->DistToBad()>0 && p2->DistToBad()>0){
756 fhRe2 [index]->Fill(pt,m) ;
398c93cc 757 if(fMakeInvPtPlots)fhReInvPt2[index]->Fill(pt,m,1./pt) ;
af7b3903 758 if(p1->DistToBad()>1 && p2->DistToBad()>1){
759 fhRe3 [index]->Fill(pt,m) ;
398c93cc 760 if(fMakeInvPtPlots)fhReInvPt3[index]->Fill(pt,m,1./pt) ;
af7b3903 761 }//assymetry cut
762 }// asymmetry cut loop
763 }// bad 3
764 }// bad2
765 }// bad 1
766 }// pid bit loop
5ae09196 767
af7b3903 768 //Fill histograms with opening angle
769 fhRealOpeningAngle ->Fill(pt,angle);
770 fhRealCosOpeningAngle->Fill(pt,TMath::Cos(angle));
771
772 //Fill histograms with pair assymmetry
773 fhRePtAsym->Fill(pt,a);
774 if(m > 0.10 && m < 0.16) fhRePtAsymPi0->Fill(pt,a);
775 if(m > 0.45 && m < 0.65) fhRePtAsymEta->Fill(pt,a);
776
777 //Multi cuts analysis
778 if(fMultiCutAna){
779 //Histograms for different PID bits selection
780 for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
52f4577d 781
af7b3903 782 if(p1->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton) &&
783 p2->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton)) fhRePIDBits[ipid]->Fill(pt,m) ;
52f4577d 784
af7b3903 785 //printf("ipt %d, ipid%d, name %s\n",ipt, ipid, fhRePtPIDCuts[ipt*fNPIDBitsBits+ipid]->GetName());
786 } // pid bit cut loop
787
788 //Several pt,ncell and asymmetry cuts
789 //Get the number of cells
790 Int_t ncell1 = 0;
791 Int_t ncell2 = 0;
792 AliVEvent * event = GetReader()->GetInputEvent();
793 if(event){
794 for(Int_t iclus = 0; iclus < event->GetNumberOfCaloClusters(); iclus++){
795 AliVCluster *cluster = event->GetCaloCluster(iclus);
5ae09196 796
af7b3903 797 Bool_t is = kFALSE;
798 if (fCalorimeter == "EMCAL" && GetReader()->IsEMCALCluster(cluster)) is = kTRUE;
799 else if(fCalorimeter == "PHOS" && GetReader()->IsPHOSCluster (cluster)) is = kTRUE;
5ae09196 800
af7b3903 801 if(is){
802 if (p1->GetCaloLabel(0) == cluster->GetID()) ncell1 = cluster->GetNCells();
803 else if (p2->GetCaloLabel(0) == cluster->GetID()) ncell2 = cluster->GetNCells();
804 } // PHOS or EMCAL cluster as requested in analysis
805
806 if(ncell2 > 0 && ncell1 > 0) break; // No need to continue the iteration
807
808 }
809 //printf("e 1: %2.2f, e 2: %2.2f, ncells: n1 %d, n2 %d\n", p1->E(), p2->E(),ncell1,ncell2);
810 }
811 for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
812 for(Int_t icell=0; icell<fNCellNCuts; icell++){
813 for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
814 Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
815 if(p1->Pt() > fPtCuts[ipt] && p2->Pt() > fPtCuts[ipt] &&
816 a < fAsymCuts[iasym] &&
817 ncell1 >= fCellNCuts[icell] && ncell2 >= fCellNCuts[icell]) fhRePtNCellAsymCuts[index]->Fill(pt,m) ;
818
819 //printf("ipt %d, icell%d, iasym %d, name %s\n",ipt, icell, iasym, fhRePtNCellAsymCuts[((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym]->GetName());
820 }// pid bit cut loop
821 }// icell loop
822 }// pt cut loop
823 for(Int_t iasym = 0; iasym < fNAsymCuts; iasym++){
824 if(a < fAsymCuts[iasym])fhRePtMult[iasym]->Fill(pt,GetTrackMultiplicity(),m) ;
825 }
826
827 }// multiple cuts analysis
828 }// ok if same sm
7e7694bb 829 }// second same event particle
830 }// first cluster
831
832 if(fDoOwnMix){
833 //Fill mixed
5025c139 834 TList * evMixList=fEventsList[curCentrBin*GetNZvertBin()*GetNRPBin()+curZvertBin*GetNRPBin()+curRPBin] ;
7e7694bb 835 Int_t nMixed = evMixList->GetSize() ;
836 for(Int_t ii=0; ii<nMixed; ii++){
837 TClonesArray* ev2= (TClonesArray*) (evMixList->At(ii));
838 Int_t nPhot2=ev2->GetEntriesFast() ;
839 Double_t m = -999;
840 if(GetDebug() > 1) printf("AliAnaPi0::MakeAnalysisFillHistograms() - Mixed event %d photon entries %d\n", ii, nPhot);
841
842 for(Int_t i1=0; i1<nPhot; i1++){
843 AliAODPWG4Particle * p1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i1)) ;
844 TLorentzVector photon1(p1->Px(),p1->Py(),p1->Pz(),p1->E());
af7b3903 845 module1 = GetModuleNumber(p1);
7e7694bb 846 for(Int_t i2=0; i2<nPhot2; i2++){
847 AliAODPWG4Particle * p2 = (AliAODPWG4Particle*) (ev2->At(i2)) ;
848
849 TLorentzVector photon2(p2->Px(),p2->Py(),p2->Pz(),p2->E());
850 m = (photon1+photon2).M() ;
851 Double_t pt = (photon1 + photon2).Pt();
852 Double_t a = TMath::Abs(p1->E()-p2->E())/(p1->E()+p2->E()) ;
853
854 //Check if opening angle is too large or too small compared to what is expected
855 Double_t angle = photon1.Angle(photon2.Vect());
856 //if(fUseAngleCut && !GetNeutralMesonSelection()->IsAngleInWindow((photon1+photon2).E(),angle)) continue;
857 if(fUseAngleCut && angle < 0.1) continue;
858
859 if(GetDebug() > 2)
860 printf("AliAnaPi0::MakeAnalysisFillHistograms() - Mixed Event: pT: photon1 %2.2f, photon2 %2.2f; Pair: pT %2.2f, mass %2.3f, a %f2.3\n",
af7b3903 861 p1->Pt(), p2->Pt(), pt,m,a);
862 //In case we want only pairs in same (super) module, check their origin.
863 module2 = GetModuleNumber(p2);
864 Bool_t ok = kTRUE;
865 if(fSameSM && module1!=module2) ok=kFALSE;
866 if(ok){
867 for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
868 if((p1->IsPIDOK(ipid,AliCaloPID::kPhoton)) && (p2->IsPIDOK(ipid,AliCaloPID::kPhoton))){
869 for(Int_t iasym=0; iasym < fNAsymCuts; iasym++){
870 if(a < fAsymCuts[iasym]){
871 Int_t index = ((curCentrBin*fNPIDBits)+ipid)*fNAsymCuts + iasym;
872 fhMi1 [index]->Fill(pt,m) ;
398c93cc 873 if(fMakeInvPtPlots)fhMiInvPt1[index]->Fill(pt,m,1./pt) ;
af7b3903 874 if(p1->DistToBad()>0 && p2->DistToBad()>0){
875 fhMi2 [index]->Fill(pt,m) ;
398c93cc 876 if(fMakeInvPtPlots)fhMiInvPt2[index]->Fill(pt,m,1./pt) ;
af7b3903 877 if(p1->DistToBad()>1 && p2->DistToBad()>1){
878 fhMi3 [index]->Fill(pt,m) ;
398c93cc 879 if(fMakeInvPtPlots)fhMiInvPt3[index]->Fill(pt,m,1./pt) ;
af7b3903 880 }
881 }
882 }//Asymmetry cut
883 }// Asymmetry loop
884 }//PID cut
885 }//loop for histograms
886 }//ok
7e7694bb 887 }// second cluster loop
888 }//first cluster loop
889 }//loop on mixed events
890
891 TClonesArray *currentEvent = new TClonesArray(*GetInputAODBranch());
af7b3903 892 //Add current event to buffer and Remove redundant events
7e7694bb 893 if(currentEvent->GetEntriesFast()>0){
894 evMixList->AddFirst(currentEvent) ;
895 currentEvent=0 ; //Now list of particles belongs to buffer and it will be deleted with buffer
896 if(evMixList->GetSize()>=fNmaxMixEv)
897 {
898 TClonesArray * tmp = (TClonesArray*) (evMixList->Last()) ;
899 evMixList->RemoveLast() ;
900 delete tmp ;
901 }
902 }
903 else{ //empty event
904 delete currentEvent ;
905 currentEvent=0 ;
477d6cee 906 }
7e7694bb 907 }// DoOwnMix
c8fe2783 908
1c5acb87 909}
910
a5cc4f03 911//________________________________________________________________________
912void AliAnaPi0::ReadHistograms(TList* outputList)
913{
50f39b97 914 // Needed when Terminate is executed in distributed environment
915 // Refill analysis histograms of this class with corresponding histograms in output list.
916
917 // Histograms of this analsys are kept in the same list as other analysis, recover the position of
918 // the first one and then add the next.
919 Int_t index = outputList->IndexOf(outputList->FindObject(GetAddedHistogramsStringToName()+"hRe_cen0_pid0_dist1"));
920
af7b3903 921 if(!fhRe1) fhRe1 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
922 if(!fhRe2) fhRe2 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
923 if(!fhRe3) fhRe3 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
924 if(!fhMi1) fhMi1 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
925 if(!fhMi2) fhMi2 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
926 if(!fhMi3) fhMi3 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
398c93cc 927 if(fMakeInvPtPlots){
928 if(!fhReInvPt1) fhReInvPt1 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
929 if(!fhReInvPt2) fhReInvPt2 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
930 if(!fhReInvPt3) fhReInvPt3 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
931 if(!fhMiInvPt1) fhMiInvPt1 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
932 if(!fhMiInvPt2) fhMiInvPt2 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
933 if(!fhMiInvPt3) fhMiInvPt3 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
934 }
af7b3903 935 if(!fhReMod) fhReMod = new TH2D*[fNModules] ;
936 if(!fhReDiffMod)fhReDiffMod = new TH2D*[fNModules+1] ;
821c8090 937
50f39b97 938 for(Int_t ic=0; ic<fNCentrBin; ic++){
af7b3903 939 for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
940 for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
941 Int_t ihisto = ((ic*fNPIDBits)+ipid)*fNAsymCuts + iasym;
942
943 fhRe1[ihisto] = (TH2D*) outputList->At(index++);
944 fhRe2[ihisto] = (TH2D*) outputList->At(index++);
945 fhRe3[ihisto] = (TH2D*) outputList->At(index++);
398c93cc 946 if(fMakeInvPtPlots){
947 fhReInvPt1[ihisto] = (TH2D*) outputList->At(index++);
948 fhReInvPt2[ihisto] = (TH2D*) outputList->At(index++);
949 fhReInvPt3[ihisto] = (TH2D*) outputList->At(index++);
950 }
5ae09196 951
af7b3903 952 if(fDoOwnMix){
953 fhMi1[ihisto] = (TH2D*) outputList->At(index++);
954 fhMi2[ihisto] = (TH2D*) outputList->At(index++);
955 fhMi3[ihisto] = (TH2D*) outputList->At(index++);
398c93cc 956 if(fMakeInvPtPlots){
957 fhMiInvPt1[ihisto] = (TH2D*) outputList->At(index++);
958 fhMiInvPt2[ihisto] = (TH2D*) outputList->At(index++);
959 fhMiInvPt3[ihisto] = (TH2D*) outputList->At(index++);
960 }
af7b3903 961 }//Own mix
962 }//asymmetry loop
963 }// pid loop
964 }// centrality loop
965
966 fhRePtAsym = (TH2D*)outputList->At(index++);
967 fhRePtAsymPi0 = (TH2D*)outputList->At(index++);
968 fhRePtAsymEta = (TH2D*)outputList->At(index++);
eee5fcf1 969
5ae09196 970 if(fMultiCutAna){
971
eee5fcf1 972 if(!fhRePtNCellAsymCuts) fhRePtNCellAsymCuts = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
973 if(!fhRePIDBits) fhRePIDBits = new TH2D*[fNPIDBits];
974
5ae09196 975 for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
976 fhRePIDBits[ipid] = (TH2D*) outputList->At(index++);
977 }// ipid loop
978
979 for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
980 for(Int_t icell=0; icell<fNCellNCuts; icell++){
981 for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
982 fhRePtNCellAsymCuts[((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym] = (TH2D*) outputList->At(index++);
983 }// iasym
984 }// icell loop
985 }// ipt loop
af7b3903 986
987 if(!fhRePtMult) fhRePtMult = new TH3D*[fNAsymCuts] ;
988 for(Int_t iasym = 0; iasym < fNAsymCuts; iasym++)
989 fhRePtMult[iasym] = (TH3D*) outputList->At(index++);
5ae09196 990 }// multi cut analysis
50f39b97 991
992 fhEvents = (TH3D *) outputList->At(index++);
993
af7b3903 994 fhRealOpeningAngle = (TH2D*) outputList->At(index++);
995 fhRealCosOpeningAngle = (TH2D*) outputList->At(index++);
996
50f39b97 997 //Histograms filled only if MC data is requested
998 if(IsDataMC() || (GetReader()->GetDataType() == AliCaloTrackReader::kMC) ){
999 fhPrimPt = (TH1D*) outputList->At(index++);
1000 fhPrimAccPt = (TH1D*) outputList->At(index++);
1001 fhPrimY = (TH1D*) outputList->At(index++);
1002 fhPrimAccY = (TH1D*) outputList->At(index++);
1003 fhPrimPhi = (TH1D*) outputList->At(index++);
1004 fhPrimAccPhi = (TH1D*) outputList->At(index++);
1005 }
1006
1007 for(Int_t imod=0; imod < fNModules; imod++)
af7b3903 1008 fhReMod[imod] = (TH2D*) outputList->At(index++);
50f39b97 1009
eee5fcf1 1010
a5cc4f03 1011}
1012
1013
6639984f 1014//____________________________________________________________________________________________________________________________________________________
a5cc4f03 1015void AliAnaPi0::Terminate(TList* outputList)
6639984f 1016{
1017 //Do some calculations and plots from the final histograms.
477d6cee 1018
fbeaf916 1019 printf(" *** %s Terminate:\n", GetName()) ;
50f39b97 1020
a5cc4f03 1021 //Recover histograms from output histograms list, needed for distributed analysis.
1022 ReadHistograms(outputList);
50f39b97 1023
2e557d1c 1024 if (!fhRe1) {
50f39b97 1025 printf("AliAnaPi0::Terminate() - Error: Remote output histograms not imported in AliAnaPi0 object");
1026 return;
2e557d1c 1027 }
50f39b97 1028
a3aebfff 1029 printf("AliAnaPi0::Terminate() Mgg Real : %5.3f , RMS : %5.3f \n", fhRe1[0]->GetMean(), fhRe1[0]->GetRMS() ) ;
5ae09196 1030
1031 const Int_t buffersize = 255;
1032
1033 char nameIM[buffersize];
1034 snprintf(nameIM, buffersize,"AliAnaPi0_%s_cPt",fCalorimeter.Data());
71dd883b 1035 TCanvas * cIM = new TCanvas(nameIM, "", 400, 10, 600, 700) ;
6639984f 1036 cIM->Divide(2, 2);
50f39b97 1037
6639984f 1038 cIM->cd(1) ;
1039 //gPad->SetLogy();
af7b3903 1040 TH1D * hIMAllPt = (TH1D*) fhRe1[0]->ProjectionY(Form("IMPtAll_%s",fCalorimeter.Data()));
6639984f 1041 hIMAllPt->SetLineColor(2);
1042 hIMAllPt->SetTitle("No cut on p_{T, #gamma#gamma} ");
1043 hIMAllPt->Draw();
1044
1045 cIM->cd(2) ;
af7b3903 1046 TH1D * hIMPt5 = (TH1D*) fhRe1[0]->ProjectionY(Form("IMPt0-5_%s",fCalorimeter.Data()),0, fhRe1[0]->GetXaxis()->FindBin(5.));
2244659d 1047// hRe1Pt5->GetXaxis()->SetRangeUser(0,5);
1048// TH1D * hIMPt5 = (TH1D*) hRe1Pt5->Project3D(Form("IMPt5_%s_pz",fCalorimeter.Data()));
6639984f 1049 hIMPt5->SetLineColor(2);
1050 hIMPt5->SetTitle("0 < p_{T, #gamma#gamma} < 5 GeV/c");
1051 hIMPt5->Draw();
1052
1053 cIM->cd(3) ;
af7b3903 1054 TH1D * hIMPt10 = (TH1D*) fhRe1[0]->ProjectionY(Form("IMPt5-10_%s",fCalorimeter.Data()), fhRe1[0]->GetXaxis()->FindBin(5.),fhRe1[0]->GetXaxis()->FindBin(10.));
2244659d 1055// hRe1Pt10->GetXaxis()->SetRangeUser(5,10);
1056// TH1D * hIMPt10 = (TH1D*) hRe1Pt10->Project3D(Form("IMPt10_%s_pz",fCalorimeter.Data()));
6639984f 1057 hIMPt10->SetLineColor(2);
1058 hIMPt10->SetTitle("5 < p_{T, #gamma#gamma} < 10 GeV/c");
1059 hIMPt10->Draw();
1060
1061 cIM->cd(4) ;
af7b3903 1062 TH1D * hIMPt20 = (TH1D*) fhRe1[0]->ProjectionY(Form("IMPt10-20_%s",fCalorimeter.Data()), fhRe1[0]->GetXaxis()->FindBin(10.),fhRe1[0]->GetXaxis()->FindBin(20.));
2244659d 1063 // TH3F * hRe1Pt20 = (TH3F*)fhRe1[0]->Clone(Form("IMPt20_%s",fCalorimeter.Data()));
1064// hRe1Pt20->GetXaxis()->SetRangeUser(10,20);
1065// TH1D * hIMPt20 = (TH1D*) hRe1Pt20->Project3D(Form("IMPt20_%s_pz",fCalorimeter.Data()));
6639984f 1066 hIMPt20->SetLineColor(2);
1067 hIMPt20->SetTitle("10 < p_{T, #gamma#gamma} < 20 GeV/c");
1068 hIMPt20->Draw();
1069
5ae09196 1070 char nameIMF[buffersize];
1071 snprintf(nameIMF,buffersize,"AliAnaPi0_%s_Mgg.eps",fCalorimeter.Data());
71dd883b 1072 cIM->Print(nameIMF);
6639984f 1073
5ae09196 1074 char namePt[buffersize];
1075 snprintf(namePt,buffersize,"AliAnaPi0_%s_cPt",fCalorimeter.Data());
71dd883b 1076 TCanvas * cPt = new TCanvas(namePt, "", 400, 10, 600, 700) ;
6639984f 1077 cPt->Divide(2, 2);
1078
1079 cPt->cd(1) ;
1080 //gPad->SetLogy();
af7b3903 1081 TH1D * hPt = (TH1D*) fhRe1[0]->ProjectionX(Form("Pt0_%s",fCalorimeter.Data()),-1,-1);
6639984f 1082 hPt->SetLineColor(2);
1083 hPt->SetTitle("No cut on M_{#gamma#gamma} ");
1084 hPt->Draw();
1085
1086 cPt->cd(2) ;
af7b3903 1087 TH1D * hPtIM1 = (TH1D*)fhRe1[0]->ProjectionX(Form("Pt1_%s",fCalorimeter.Data()), fhRe1[0]->GetZaxis()->FindBin(0.05),fhRe1[0]->GetZaxis()->FindBin(0.21));
2244659d 1088// TH3F * hRe1IM1 = (TH3F*)fhRe1[0]->Clone(Form("Pt1_%s",fCalorimeter.Data()));
1089// hRe1IM1->GetZaxis()->SetRangeUser(0.05,0.21);
1090// TH1D * hPtIM1 = (TH1D*) hRe1IM1->Project3D("x");
6639984f 1091 hPtIM1->SetLineColor(2);
1092 hPtIM1->SetTitle("0.05 < M_{#gamma#gamma} < 0.21 GeV/c^{2}");
1093 hPtIM1->Draw();
1094
1095 cPt->cd(3) ;
af7b3903 1096 TH1D * hPtIM2 = (TH1D*)fhRe1[0]->ProjectionX(Form("Pt2_%s",fCalorimeter.Data()), fhRe1[0]->GetZaxis()->FindBin(0.09),fhRe1[0]->GetZaxis()->FindBin(0.17));
2244659d 1097// TH3F * hRe1IM2 = (TH3F*)fhRe1[0]->Clone(Form("Pt2_%s",fCalorimeter.Data()));
1098// hRe1IM2->GetZaxis()->SetRangeUser(0.09,0.17);
1099// TH1D * hPtIM2 = (TH1D*) hRe1IM2->Project3D("x");
6639984f 1100 hPtIM2->SetLineColor(2);
1101 hPtIM2->SetTitle("0.09 < M_{#gamma#gamma} < 0.17 GeV/c^{2}");
1102 hPtIM2->Draw();
1103
1104 cPt->cd(4) ;
af7b3903 1105 TH1D * hPtIM3 = (TH1D*)fhRe1[0]->ProjectionX(Form("Pt3_%s",fCalorimeter.Data()), fhRe1[0]->GetZaxis()->FindBin(0.11),fhRe1[0]->GetZaxis()->FindBin(0.15));
2244659d 1106// TH3F * hRe1IM3 = (TH3F*)fhRe1[0]->Clone(Form("Pt3_%s",fCalorimeter.Data()));
1107// hRe1IM3->GetZaxis()->SetRangeUser(0.11,0.15);
1108// TH1D * hPtIM3 = (TH1D*) hRe1IM1->Project3D("x");
6639984f 1109 hPtIM3->SetLineColor(2);
1110 hPtIM3->SetTitle("0.11 < M_{#gamma#gamma} < 0.15 GeV/c^{2}");
1111 hPtIM3->Draw();
1112
164a1d84 1113 char namePtF[buffersize];
5ae09196 1114 snprintf(namePtF,buffersize,"AliAnaPi0_%s_Pt.eps",fCalorimeter.Data());
71dd883b 1115 cPt->Print(namePtF);
1c5acb87 1116
5ae09196 1117 char line[buffersize] ;
1118 snprintf(line,buffersize,".!tar -zcf %s_%s.tar.gz *.eps", GetName(),fCalorimeter.Data()) ;
6639984f 1119 gROOT->ProcessLine(line);
5ae09196 1120 snprintf(line, buffersize,".!rm -fR AliAnaPi0_%s*.eps",fCalorimeter.Data());
6639984f 1121 gROOT->ProcessLine(line);
1122
71dd883b 1123 printf(" AliAnaPi0::Terminate() - !! All the eps files are in %s_%s.tar.gz !!!\n", GetName(), fCalorimeter.Data());
1c5acb87 1124
6639984f 1125}
c8fe2783 1126 //____________________________________________________________________________________________________________________________________________________
1127Int_t AliAnaPi0::GetEventIndex(AliAODPWG4Particle * part, Double_t * vert)
1128{
f8006433 1129 // retieves the event index and checks the vertex
1130 // in the mixed buffer returns -2 if vertex NOK
1131 // for normal events returns 0 if vertex OK and -1 if vertex NOK
1132
1133 Int_t evtIndex = -1 ;
1134 if(GetReader()->GetDataType()!=AliCaloTrackReader::kMC){
1135
1136 if (GetMixedEvent()){
1137
1138 evtIndex = GetMixedEvent()->EventIndexForCaloCluster(part->GetCaloLabel(0)) ;
1139 GetVertex(vert,evtIndex);
1140
5025c139 1141 if(TMath::Abs(vert[2])> GetZvertexCut())
f8006433 1142 evtIndex = -2 ; //Event can not be used (vertex, centrality,... cuts not fulfilled)
1143 } else {// Single event
1144
1145 GetVertex(vert);
1146
5025c139 1147 if(TMath::Abs(vert[2])> GetZvertexCut())
f8006433 1148 evtIndex = -1 ; //Event can not be used (vertex, centrality,... cuts not fulfilled)
1149 else
1150 evtIndex = 0 ;
c8fe2783 1151 }
0ae57829 1152 }//No MC reader
f8006433 1153 else {
1154 evtIndex = 0;
1155 vert[0] = 0. ;
1156 vert[1] = 0. ;
1157 vert[2] = 0. ;
1158 }
0ae57829 1159
f8006433 1160 return evtIndex ;
c8fe2783 1161}
f8006433 1162