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