]> 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
AliAnaPi0: Added option to fill mass histograms with different cuts
[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(),
7e7694bb 56fDoOwnMix(kFALSE),fNCentrBin(0),fNZvertBin(0),fNrpBin(0),
1c5acb87 57fNPID(0),fNmaxMixEv(0), fZvtxCut(0.),fCalorimeter(""),
5ae09196 58fNModules(12), fUseAngleCut(kFALSE), fEventsList(0x0), fMultiCutAna(kFALSE),
59fNPtCuts(0),fPtCuts(0x0),fNAsymCuts(0),fAsymCuts(0x0),
60fNCellNCuts(0),fCellNCuts(0x0),fNPIDBits(0),fPIDBits(0x0),
61fHistoInvPtBins(0), fHistoInvPtMax(0), fHistoInvPtMin(0), fhReMod(0x0),
62fhRe1(0x0), fhMi1(0x0), fhRe2(0x0), fhMi2(0x0), fhRe3(0x0), fhMi3(0x0),
63fhReInvPt1(0x0), fhMiInvPt1(0x0), fhReInvPt2(0x0), fhMiInvPt2(0x0), fhReInvPt3(0x0), fhMiInvPt3(0x0),
64fhRePtNCellAsymCuts(0x0), fhRePIDBits(0x0),
65fhEvents(0x0), fhRealOpeningAngle(0x0),fhRealCosOpeningAngle(0x0),
50f39b97 66fhPrimPt(0x0), fhPrimAccPt(0x0), fhPrimY(0x0), fhPrimAccY(0x0), fhPrimPhi(0x0), fhPrimAccPhi(0x0),
67fhPrimOpeningAngle(0x0),fhPrimCosOpeningAngle(0x0)
1c5acb87 68{
69//Default Ctor
70 InitParameters();
71
72}
1c5acb87 73
74//________________________________________________________________________________________________________________________________________________
75AliAnaPi0::~AliAnaPi0() {
477d6cee 76 // Remove event containers
7e7694bb 77
78 if(fDoOwnMix && fEventsList){
477d6cee 79 for(Int_t ic=0; ic<fNCentrBin; ic++){
80 for(Int_t iz=0; iz<fNZvertBin; iz++){
7e7694bb 81 for(Int_t irp=0; irp<fNrpBin; irp++){
82 fEventsList[ic*fNZvertBin*fNrpBin+iz*fNrpBin+irp]->Delete() ;
83 delete fEventsList[ic*fNZvertBin*fNrpBin+iz*fNrpBin+irp] ;
84 }
477d6cee 85 }
86 }
87 delete[] fEventsList;
88 fEventsList=0 ;
89 }
591cc579 90
1c5acb87 91}
92
93//________________________________________________________________________________________________________________________________________________
94void AliAnaPi0::InitParameters()
95{
96//Init parameters when first called the analysis
97//Set default parameters
a3aebfff 98 SetInputAODName("PWG4Particle");
99
100 AddToHistogramsName("AnaPi0_");
6921fa00 101 fNModules = 12; // set maximum to maximum number of EMCAL modules
477d6cee 102 fNCentrBin = 1;
103 fNZvertBin = 1;
104 fNrpBin = 1;
105 fNPID = 9;
106 fNmaxMixEv = 10;
107 fZvtxCut = 40;
108 fCalorimeter = "PHOS";
50f39b97 109 fUseAngleCut = kFALSE;
5ae09196 110
111 fMultiCutAna = kFALSE;
112
113 fNPtCuts = 3;
114 fPtCuts = new Float_t[fNPtCuts];
115 fPtCuts[0] = 0.; fPtCuts[1] = 0.2; fPtCuts[2] = 0.3;
116
117 fNAsymCuts = 3;
118 fAsymCuts = new Float_t[fNAsymCuts];
119 fAsymCuts[0] = 0.7; fAsymCuts[2] = 0.8; fAsymCuts[2] = 1.;
120
121 fNCellNCuts = 3;
122 fCellNCuts = new Int_t[fNCellNCuts];
123 fCellNCuts[0] = 1; fCellNCuts[1] = 2; fCellNCuts[2] = 3;
124
125 fNPIDBits = 3;
126 fPIDBits = new Int_t[fNPIDBits];
127 fPIDBits[0] = 2; fPIDBits[1] = 4; fPIDBits[2] = 6; // check dispersion, neutral, dispersion&&neutral
128
129 fHistoInvPtMax = 10.;
130 fHistoInvPtMin = 0.;
131 fHistoInvPtBins = 200;
132
1c5acb87 133}
1c5acb87 134
0c1383b5 135
136//________________________________________________________________________________________________________________________________________________
137TObjString * AliAnaPi0::GetAnalysisCuts()
138{
139 //Save parameters used for analysis
140 TString parList ; //this will be list of parameters used for this analysis.
5ae09196 141 const Int_t buffersize = 255;
142 char onePar[buffersize] ;
143 snprintf(onePar,buffersize,"--- AliAnaPi0 ---\n") ;
0c1383b5 144 parList+=onePar ;
5ae09196 145 snprintf(onePar,buffersize,"Number of bins in Centrality: %d \n",fNCentrBin) ;
0c1383b5 146 parList+=onePar ;
5ae09196 147 snprintf(onePar,buffersize,"Number of bins in Z vert. pos: %d \n",fNZvertBin) ;
0c1383b5 148 parList+=onePar ;
5ae09196 149 snprintf(onePar,buffersize,"Number of bins in Reac. Plain: %d \n",fNrpBin) ;
0c1383b5 150 parList+=onePar ;
5ae09196 151 snprintf(onePar,buffersize,"Depth of event buffer: %d \n",fNmaxMixEv) ;
0c1383b5 152 parList+=onePar ;
5ae09196 153 snprintf(onePar,buffersize,"Number of different PID used: %d \n",fNPID) ;
0c1383b5 154 parList+=onePar ;
5ae09196 155 snprintf(onePar,buffersize,"Cuts: \n") ;
0c1383b5 156 parList+=onePar ;
5ae09196 157 snprintf(onePar,buffersize,"Z vertex position: -%f < z < %f \n",fZvtxCut,fZvtxCut) ;
0c1383b5 158 parList+=onePar ;
5ae09196 159 snprintf(onePar,buffersize,"Calorimeter: %s \n",fCalorimeter.Data()) ;
0c1383b5 160 parList+=onePar ;
5ae09196 161 snprintf(onePar,buffersize,"Number of modules: %d \n",fNModules) ;
0c1383b5 162 parList+=onePar ;
163
164 return new TObjString(parList) ;
165}
166
2e557d1c 167//________________________________________________________________________________________________________________________________________________
168TList * AliAnaPi0::GetCreateOutputObjects()
169{
477d6cee 170 // Create histograms to be saved in output file and
171 // store them in fOutputContainer
172
173 //create event containers
174 fEventsList = new TList*[fNCentrBin*fNZvertBin*fNrpBin] ;
1c5acb87 175
477d6cee 176 for(Int_t ic=0; ic<fNCentrBin; ic++){
177 for(Int_t iz=0; iz<fNZvertBin; iz++){
178 for(Int_t irp=0; irp<fNrpBin; irp++){
7e7694bb 179 fEventsList[ic*fNZvertBin*fNrpBin+iz*fNrpBin+irp] = new TList() ;
477d6cee 180 }
181 }
182 }
7e7694bb 183
477d6cee 184 TList * outputContainer = new TList() ;
185 outputContainer->SetName(GetName());
6921fa00 186
187 fhReMod = new TH3D*[fNModules] ;
50305ae9 188 fhRe1 = new TH3D*[fNCentrBin*fNPID] ;
189 fhRe2 = new TH3D*[fNCentrBin*fNPID] ;
190 fhRe3 = new TH3D*[fNCentrBin*fNPID] ;
191 fhMi1 = new TH3D*[fNCentrBin*fNPID] ;
192 fhMi2 = new TH3D*[fNCentrBin*fNPID] ;
193 fhMi3 = new TH3D*[fNCentrBin*fNPID] ;
5ae09196 194
195 fhReInvPt1 = new TH3D*[fNCentrBin*fNPID] ;
196 fhReInvPt2 = new TH3D*[fNCentrBin*fNPID] ;
197 fhReInvPt3 = new TH3D*[fNCentrBin*fNPID] ;
198 fhMiInvPt1 = new TH3D*[fNCentrBin*fNPID] ;
199 fhMiInvPt2 = new TH3D*[fNCentrBin*fNPID] ;
200 fhMiInvPt3 = new TH3D*[fNCentrBin*fNPID] ;
201
202 const Int_t buffersize = 255;
203 char key[buffersize] ;
204 char title[buffersize] ;
477d6cee 205
5a2dbc3c 206 Int_t nptbins = GetHistoPtBins();
5ae09196 207 Int_t niptbins = GetHistoInvPtBins();
5a2dbc3c 208 Int_t nphibins = GetHistoPhiBins();
209 Int_t netabins = GetHistoEtaBins();
210 Float_t ptmax = GetHistoPtMax();
5ae09196 211 Float_t iptmax = GetHistoInvPtMax();
5a2dbc3c 212 Float_t phimax = GetHistoPhiMax();
213 Float_t etamax = GetHistoEtaMax();
214 Float_t ptmin = GetHistoPtMin();
5ae09196 215 Float_t iptmin = GetHistoInvPtMin();
5a2dbc3c 216 Float_t phimin = GetHistoPhiMin();
217 Float_t etamin = GetHistoEtaMin();
218
219 Int_t nmassbins = GetHistoMassBins();
220 Int_t nasymbins = GetHistoAsymmetryBins();
221 Float_t massmax = GetHistoMassMax();
222 Float_t asymmax = GetHistoAsymmetryMax();
223 Float_t massmin = GetHistoMassMin();
224 Float_t asymmin = GetHistoAsymmetryMin();
225
477d6cee 226 for(Int_t ic=0; ic<fNCentrBin; ic++){
227 for(Int_t ipid=0; ipid<fNPID; ipid++){
7e7694bb 228
477d6cee 229 //Distance to bad module 1
5ae09196 230 snprintf(key, buffersize,"hRe_cen%d_pid%d_dist1",ic,ipid) ;
231 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
7e7694bb 232 fhRe1[ic*fNPID+ipid] = new TH3D(key,title,nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax,nmassbins,massmin,massmax) ;
477d6cee 233 outputContainer->Add(fhRe1[ic*fNPID+ipid]) ;
7e7694bb 234
477d6cee 235 //Distance to bad module 2
5ae09196 236 snprintf(key, buffersize,"hRe_cen%d_pid%d_dist2",ic,ipid) ;
237 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
7e7694bb 238 fhRe2[ic*fNPID+ipid] = new TH3D(key,title,nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax,nmassbins,massmin,massmax) ;
477d6cee 239 outputContainer->Add(fhRe2[ic*fNPID+ipid]) ;
240
477d6cee 241 //Distance to bad module 3
5ae09196 242 snprintf(key, buffersize,"hRe_cen%d_pid%d_dist3",ic,ipid) ;
243 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
7e7694bb 244 fhRe3[ic*fNPID+ipid] = new TH3D(key,title,nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax,nmassbins,massmin,massmax) ;
477d6cee 245 outputContainer->Add(fhRe3[ic*fNPID+ipid]) ;
246
5ae09196 247 //Inverse pT
248 //Distance to bad module 1
249 snprintf(key, buffersize,"hReInvPt_cen%d_pid%d_dist1",ic,ipid) ;
250 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
251 fhReInvPt1[ic*fNPID+ipid] = new TH3D(key,title,niptbins,iptmin,iptmax,nasymbins,asymmin,asymmax,nmassbins,massmin,massmax) ;
252 outputContainer->Add(fhReInvPt1[ic*fNPID+ipid]) ;
253
254 //Distance to bad module 2
255 snprintf(key, buffersize,"hReInvPt_cen%d_pid%d_dist2",ic,ipid) ;
256 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
257 fhReInvPt2[ic*fNPID+ipid] = new TH3D(key,title,niptbins,iptmin,iptmax,nasymbins,asymmin,asymmax,nmassbins,massmin,massmax) ;
258 outputContainer->Add(fhReInvPt2[ic*fNPID+ipid]) ;
259
260 //Distance to bad module 3
261 snprintf(key, buffersize,"hReInvPt_cen%d_pid%d_dist3",ic,ipid) ;
262 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
263 fhReInvPt3[ic*fNPID+ipid] = new TH3D(key,title,niptbins,iptmin,iptmax,nasymbins,asymmin,asymmax,nmassbins,massmin,massmax) ;
264 outputContainer->Add(fhReInvPt3[ic*fNPID+ipid]) ;
265
266
7e7694bb 267 if(fDoOwnMix){
268 //Distance to bad module 1
5ae09196 269 snprintf(key, buffersize,"hMi_cen%d_pid%d_dist1",ic,ipid) ;
270 snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
7e7694bb 271 fhMi1[ic*fNPID+ipid] = new TH3D(key,title,nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax,nmassbins,massmin,massmax) ;
272 outputContainer->Add(fhMi1[ic*fNPID+ipid]) ;
273
274 //Distance to bad module 2
5ae09196 275 snprintf(key, buffersize,"hMi_cen%d_pid%d_dist2",ic,ipid) ;
276 snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
7e7694bb 277 fhMi2[ic*fNPID+ipid] = new TH3D(key,title,nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax,nmassbins,massmin,massmax) ;
278 outputContainer->Add(fhMi2[ic*fNPID+ipid]) ;
279
280 //Distance to bad module 3
5ae09196 281 snprintf(key, buffersize,"hMi_cen%d_pid%d_dist3",ic,ipid) ;
282 snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
7e7694bb 283 fhMi3[ic*fNPID+ipid] = new TH3D(key,title,nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax,nmassbins,massmin,massmax) ;
284 outputContainer->Add(fhMi3[ic*fNPID+ipid]) ;
5ae09196 285
286 //Inverse pT
287 //Distance to bad module 1
288 snprintf(key, buffersize,"hMiInvPt_cen%d_pid%d_dist1",ic,ipid) ;
289 snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
290 fhMiInvPt1[ic*fNPID+ipid] = new TH3D(key,title,niptbins,iptmin,iptmax,nasymbins,asymmin,asymmax,nmassbins,massmin,massmax) ;
291 outputContainer->Add(fhMiInvPt1[ic*fNPID+ipid]) ;
292
293 //Distance to bad module 2
294 snprintf(key, buffersize,"hMiInvPt_cen%d_pid%d_dist2",ic,ipid) ;
295 snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
296 fhMiInvPt2[ic*fNPID+ipid] = new TH3D(key,title,niptbins,iptmin,iptmax,nasymbins,asymmin,asymmax,nmassbins,massmin,massmax) ;
297 outputContainer->Add(fhMiInvPt2[ic*fNPID+ipid]) ;
298
299 //Distance to bad module 3
300 snprintf(key, buffersize,"hMiInvPt_cen%d_pid%d_dist3",ic,ipid) ;
301 snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
302 fhMiInvPt3[ic*fNPID+ipid] = new TH3D(key,title,niptbins,iptmin,iptmax,nasymbins,asymmin,asymmax,nmassbins,massmin,massmax) ;
303 outputContainer->Add(fhMiInvPt3[ic*fNPID+ipid]) ;
304
305
7e7694bb 306 }
1c5acb87 307 }
477d6cee 308 }
309
5ae09196 310 if(fMultiCutAna){
311
312 fhRePIDBits = new TH2D*[fNPIDBits];
313 for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
314 snprintf(key, buffersize,"hRe_pidbit%d",ipid) ;
315 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for PIDBit=%d",fPIDBits[ipid]) ;
316 fhRePIDBits[ipid] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
317 outputContainer->Add(fhRePIDBits[ipid]) ;
318 }// pid bit loop
319
320 fhRePtNCellAsymCuts = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
321 for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
322 for(Int_t icell=0; icell<fNCellNCuts; icell++){
323 for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
324 snprintf(key, buffersize,"hRe_pt%d_cell%d_asym%d",ipt,icell,iasym) ;
325 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for pt >%2.2f, ncell>%d and asym >%2.2f ",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]) ;
326 Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
327 //printf("ipt %d, icell %d, iassym %d, index %d\n",ipt, icell, iasym, index);
328 fhRePtNCellAsymCuts[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
329 outputContainer->Add(fhRePtNCellAsymCuts[index]) ;
330 }
331 }
332 }
333 }// multi cuts analysis
334
477d6cee 335 fhEvents=new TH3D("hEvents","Number of events",fNCentrBin,0.,1.*fNCentrBin,
7e7694bb 336 fNZvertBin,0.,1.*fNZvertBin,fNrpBin,0.,1.*fNrpBin) ;
477d6cee 337 outputContainer->Add(fhEvents) ;
50f39b97 338
339 fhRealOpeningAngle = new TH2D
340 ("hRealOpeningAngle","Angle between all #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,200,0,0.5);
341 fhRealOpeningAngle->SetYTitle("#theta(rad)");
342 fhRealOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
343 outputContainer->Add(fhRealOpeningAngle) ;
7e7694bb 344
50f39b97 345 fhRealCosOpeningAngle = new TH2D
346 ("hRealCosOpeningAngle","Cosinus of angle between all #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,200,-1,1);
347 fhRealCosOpeningAngle->SetYTitle("cos (#theta) ");
348 fhRealCosOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
349 outputContainer->Add(fhRealCosOpeningAngle) ;
350
351
477d6cee 352 //Histograms filled only if MC data is requested
0ae57829 353 if(IsDataMC()){
5ae09196 354
7e7694bb 355 fhPrimPt = new TH1D("hPrimPt","Primary pi0 pt",nptbins,ptmin,ptmax) ;
5a2dbc3c 356 fhPrimAccPt = new TH1D("hPrimAccPt","Primary pi0 pt with both photons in acceptance",nptbins,ptmin,ptmax) ;
477d6cee 357 outputContainer->Add(fhPrimPt) ;
358 outputContainer->Add(fhPrimAccPt) ;
359
5a2dbc3c 360 fhPrimY = new TH1D("hPrimaryRapidity","Rapidity of primary pi0",netabins,etamin,etamax) ;
477d6cee 361 outputContainer->Add(fhPrimY) ;
362
5a2dbc3c 363 fhPrimAccY = new TH1D("hPrimAccRapidity","Rapidity of primary pi0",netabins,etamin,etamax) ;
477d6cee 364 outputContainer->Add(fhPrimAccY) ;
365
7e7694bb 366 fhPrimPhi = new TH1D("hPrimaryPhi","Azimithal of primary pi0",nphibins,phimin*TMath::RadToDeg(),phimax*TMath::RadToDeg()) ;
477d6cee 367 outputContainer->Add(fhPrimPhi) ;
368
5a2dbc3c 369 fhPrimAccPhi = new TH1D("hPrimAccPhi","Azimithal of primary pi0 with accepted daughters",nphibins,phimin*TMath::RadToDeg(),phimax*TMath::RadToDeg()) ;
477d6cee 370 outputContainer->Add(fhPrimAccPhi) ;
50f39b97 371
372
373 fhPrimOpeningAngle = new TH2D
7e7694bb 374 ("hPrimOpeningAngle","Angle between all primary #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,100,0,0.5);
50f39b97 375 fhPrimOpeningAngle->SetYTitle("#theta(rad)");
376 fhPrimOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
377 outputContainer->Add(fhPrimOpeningAngle) ;
378
379 fhPrimCosOpeningAngle = new TH2D
7e7694bb 380 ("hPrimCosOpeningAngle","Cosinus of angle between all primary #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,100,-1,1);
50f39b97 381 fhPrimCosOpeningAngle->SetYTitle("cos (#theta) ");
382 fhPrimCosOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
383 outputContainer->Add(fhPrimCosOpeningAngle) ;
384
477d6cee 385 }
50f39b97 386
6921fa00 387 for(Int_t imod=0; imod<fNModules; imod++){
50f39b97 388 //Module dependent invariant mass
5ae09196 389 snprintf(key, buffersize,"hReMod_%d",imod) ;
390 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for Module %d",imod) ;
50f39b97 391 fhReMod[imod] = new TH3D(key,title,nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax,nmassbins,massmin,massmax) ;
392 outputContainer->Add(fhReMod[imod]) ;
6921fa00 393 }
50f39b97 394
477d6cee 395 return outputContainer;
1c5acb87 396}
397
398//_________________________________________________________________________________________________________________________________________________
399void AliAnaPi0::Print(const Option_t * /*opt*/) const
400{
477d6cee 401 //Print some relevant parameters set for the analysis
a3aebfff 402 printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
477d6cee 403 AliAnaPartCorrBaseClass::Print(" ");
a3aebfff 404
477d6cee 405 printf("Number of bins in Centrality: %d \n",fNCentrBin) ;
406 printf("Number of bins in Z vert. pos: %d \n",fNZvertBin) ;
407 printf("Number of bins in Reac. Plain: %d \n",fNrpBin) ;
408 printf("Depth of event buffer: %d \n",fNmaxMixEv) ;
409 printf("Number of different PID used: %d \n",fNPID) ;
410 printf("Cuts: \n") ;
411 printf("Z vertex position: -%2.3f < z < %2.3f \n",fZvtxCut,fZvtxCut) ;
50f39b97 412 printf("Number of modules: %d \n",fNModules) ;
413 printf("Select pairs with their angle: %d \n",fUseAngleCut) ;
477d6cee 414 printf("------------------------------------------------------\n") ;
1c5acb87 415}
416
5ae09196 417//_____________________________________________________________
418void AliAnaPi0::FillAcceptanceHistograms(){
419 //Fill acceptance histograms if MC data is available
c8fe2783 420
5ae09196 421 if(IsDataMC() && GetReader()->ReadStack()){
422 AliStack * stack = GetMCStack();
423 if(stack && (IsDataMC() || (GetReader()->GetDataType() == AliCaloTrackReader::kMC)) ){
424 for(Int_t i=0 ; i<stack->GetNprimary(); i++){
425 TParticle * prim = stack->Particle(i) ;
426 if(prim->GetPdgCode() == 111){
427 Double_t pi0Pt = prim->Pt() ;
428 //printf("pi0, pt %2.2f\n",pi0Pt);
429 if(prim->Energy() == TMath::Abs(prim->Pz())) continue ; //Protection against floating point exception
430 Double_t pi0Y = 0.5*TMath::Log((prim->Energy()-prim->Pz())/(prim->Energy()+prim->Pz())) ;
431 Double_t phi = TMath::RadToDeg()*prim->Phi() ;
432 if(TMath::Abs(pi0Y) < 0.5){
433 fhPrimPt->Fill(pi0Pt) ;
434 }
435 fhPrimY ->Fill(pi0Y) ;
436 fhPrimPhi->Fill(phi) ;
437
438 //Check if both photons hit Calorimeter
439 Int_t iphot1=prim->GetFirstDaughter() ;
440 Int_t iphot2=prim->GetLastDaughter() ;
441 if(iphot1>-1 && iphot1<stack->GetNtrack() && iphot2>-1 && iphot2<stack->GetNtrack()){
442 TParticle * phot1 = stack->Particle(iphot1) ;
443 TParticle * phot2 = stack->Particle(iphot2) ;
444 if(phot1 && phot2 && phot1->GetPdgCode()==22 && phot2->GetPdgCode()==22){
445 //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",
446 // phot1->Pt(), phot1->Phi()*180./3.1415, phot1->Eta(), phot2->Pt(), phot2->Phi()*180./3.1415, phot2->Eta());
447
448 TLorentzVector lv1, lv2;
449 phot1->Momentum(lv1);
450 phot2->Momentum(lv2);
451
452 Bool_t inacceptance = kFALSE;
453 if(fCalorimeter == "PHOS"){
454 if(GetPHOSGeometry() && GetCaloUtils()->IsPHOSGeoMatrixSet()){
455 Int_t mod ;
456 Double_t x,z ;
457 if(GetPHOSGeometry()->ImpactOnEmc(phot1,mod,z,x) && GetPHOSGeometry()->ImpactOnEmc(phot2,mod,z,x))
458 inacceptance = kTRUE;
459 if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
460 }
461 else{
462
463 if(GetFiducialCut()->IsInFiducialCut(lv1,fCalorimeter) && GetFiducialCut()->IsInFiducialCut(lv2,fCalorimeter))
464 inacceptance = kTRUE ;
465 if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
466 }
467
468 }
469 else if(fCalorimeter == "EMCAL" && GetCaloUtils()->IsEMCALGeoMatrixSet()){
470 if(GetEMCALGeometry()){
471 if(GetEMCALGeometry()->Impact(phot1) && GetEMCALGeometry()->Impact(phot2))
472 inacceptance = kTRUE;
473 if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
474 }
475 else{
476 if(GetFiducialCut()->IsInFiducialCut(lv1,fCalorimeter) && GetFiducialCut()->IsInFiducialCut(lv2,fCalorimeter))
477 inacceptance = kTRUE ;
478 if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
479 }
480 }
481
482 if(inacceptance){
483
484 fhPrimAccPt->Fill(pi0Pt) ;
485 fhPrimAccPhi->Fill(phi) ;
486 fhPrimAccY->Fill(pi0Y) ;
487 Double_t angle = lv1.Angle(lv2.Vect());
488 fhPrimOpeningAngle ->Fill(pi0Pt,angle);
489 fhPrimCosOpeningAngle->Fill(pi0Pt,TMath::Cos(angle));
490
491 }//Accepted
492 }// 2 photons
493 }//Check daughters exist
494 }// Primary pi0
495 }//loop on primaries
496 }//stack exists and data is MC
497 }//read stack
498 else if(GetReader()->ReadAODMCParticles()){
499 if(GetDebug() >= 0) printf("AliAnaPi0::FillAcceptanceHistograms() - Acceptance calculation with MCParticles not implemented yet\n");
500 }
501}
502
1c5acb87 503//____________________________________________________________________________________________________________________________________________________
6639984f 504void AliAnaPi0::MakeAnalysisFillHistograms()
1c5acb87 505{
477d6cee 506 //Process one event and extract photons from AOD branch
507 // filled with AliAnaPhoton and fill histos with invariant mass
508
5ae09196 509 //In case of MC data, fill acceptance histograms
510 FillAcceptanceHistograms();
511
477d6cee 512 //Apply some cuts on event: vertex position and centrality range
513 Int_t iRun=(GetReader()->GetInputEvent())->GetRunNumber() ;
514 if(IsBadRun(iRun)) return ;
515
477d6cee 516 Int_t nPhot = GetInputAODBranch()->GetEntriesFast() ;
c8fe2783 517 if(GetDebug() > 1)
518 printf("AliAnaPi0::MakeAnalysisFillHistograms() - Photon entries %d\n", nPhot);
519 if(nPhot < 2 )
520 return ;
6921fa00 521 Int_t module1 = -1;
522 Int_t module2 = -1;
7e7694bb 523 Double_t vert[] = {0.0, 0.0, 0.0} ; //vertex
524 Int_t evtIndex1 = 0 ;
525 Int_t currentEvtIndex = -1 ;
526 Int_t curCentrBin = 0 ;
527 Int_t curRPBin = 0 ;
528 Int_t curZvertBin = 0 ;
529
477d6cee 530 for(Int_t i1=0; i1<nPhot-1; i1++){
531 AliAODPWG4Particle * p1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i1)) ;
7e7694bb 532 // get the event index in the mixed buffer where the photon comes from
533 // in case of mixing with analysis frame, not own mixing
c8fe2783 534 evtIndex1 = GetEventIndex(p1, vert) ;
7e7694bb 535 if(vert[2]<-fZvtxCut || vert[2]> fZvtxCut) return ; //Event can not be used (vertex, centrality,... cuts not fulfilled)
c8fe2783 536 if ( evtIndex1 == -1 )
537 return ;
538 if ( evtIndex1 == -2 )
539 continue ;
540 if (evtIndex1 != currentEvtIndex) {
7e7694bb 541 //Get Reaction Plan position and calculate RP bin
542 //does not exist in ESD yet????
c8fe2783 543 curCentrBin = 0 ;
544 curRPBin = 0 ;
545 curZvertBin = (Int_t)(0.5*fNZvertBin*(vert[2]+fZvtxCut)/fZvtxCut) ;
546 fhEvents->Fill(curCentrBin+0.5,curZvertBin+0.5,curRPBin+0.5) ;
547 currentEvtIndex = evtIndex1 ;
548 }
7e7694bb 549
477d6cee 550 TLorentzVector photon1(p1->Px(),p1->Py(),p1->Pz(),p1->E());
59b6bd99 551 //Get Module number
552 module1 = GetModuleNumber(p1);
477d6cee 553 for(Int_t i2=i1+1; i2<nPhot; i2++){
554 AliAODPWG4Particle * p2 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i2)) ;
c8fe2783 555 Int_t evtIndex2 = GetEventIndex(p2, vert) ;
556 if ( evtIndex2 == -1 )
557 return ;
558 if ( evtIndex2 == -2 )
559 continue ;
560 if (GetMixedEvent() && (evtIndex1 == evtIndex2))
7e7694bb 561 continue ;
477d6cee 562 TLorentzVector photon2(p2->Px(),p2->Py(),p2->Pz(),p2->E());
59b6bd99 563 //Get module number
564 module2 = GetModuleNumber(p2);
477d6cee 565 Double_t m = (photon1 + photon2).M() ;
566 Double_t pt = (photon1 + photon2).Pt();
567 Double_t a = TMath::Abs(p1->E()-p2->E())/(p1->E()+p2->E()) ;
568 if(GetDebug() > 2)
c8fe2783 569 printf("AliAnaPi0::MakeAnalysisFillHistograms() - Current Event: pT: photon1 %2.2f, photon2 %2.2f; Pair: pT %2.2f, mass %2.3f, a %f2.3\n",
570 p1->Pt(), p2->Pt(), pt,m,a);
50f39b97 571 //Check if opening angle is too large or too small compared to what is expected
572 Double_t angle = photon1.Angle(photon2.Vect());
573 //if(fUseAngleCut && !GetNeutralMesonSelection()->IsAngleInWindow((photon1+photon2).E(),angle)) continue;
574 //printf("angle %f\n",angle);
c8fe2783 575 if(fUseAngleCut && angle < 0.1)
576 continue;
50f39b97 577 fhRealOpeningAngle ->Fill(pt,angle);
578 fhRealCosOpeningAngle->Fill(pt,TMath::Cos(angle));
59b6bd99 579 //Fill module dependent histograms
580 //if(module1==module2) printf("mod1 %d\n",module1);
581 if(module1==module2 && module1 >=0 && module1<fNModules)
c8fe2783 582 fhReMod[module1]->Fill(pt,a,m) ;
7e7694bb 583
c8fe2783 584 for(Int_t ipid=0; ipid<fNPID; ipid++){
585 if((p1->IsPIDOK(ipid,AliCaloPID::kPhoton)) && (p2->IsPIDOK(ipid,AliCaloPID::kPhoton))){
5ae09196 586 fhRe1 [curCentrBin*fNPID+ipid]->Fill(pt, a,m) ;
587 fhReInvPt1[curCentrBin*fNPID+ipid]->Fill(1./pt,a,m) ;
c8fe2783 588 if(p1->DistToBad()>0 && p2->DistToBad()>0){
5ae09196 589 fhRe2 [curCentrBin*fNPID+ipid]->Fill(pt, a,m) ;
590 fhReInvPt2[curCentrBin*fNPID+ipid]->Fill(1./pt,a,m) ;
c8fe2783 591 if(p1->DistToBad()>1 && p2->DistToBad()>1){
5ae09196 592 fhRe3 [curCentrBin*fNPID+ipid]->Fill(pt, a,m) ;
593 fhReInvPt3[curCentrBin*fNPID+ipid]->Fill(1./pt,a,m) ;
594 }// bad 3
595 }// bad2
596 }// bad 1
597 }// pid loop
598
599 //Multi cuts analysis
600 if(fMultiCutAna){
601 //Histograms for different PID bits selection
602 for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
603
604 if(p1->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton) &&
605 p2->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton)) fhRePIDBits[ipid]->Fill(pt,m) ;
606
607 //printf("ipt %d, ipid%d, name %s\n",ipt, ipid, fhRePtPIDCuts[ipt*fNPIDBits+ipid]->GetName());
608 } // pid bit cut loop
609
610 //Several pt,ncell and asymetry cuts
611 //Get the number of cells
612 Int_t ncell1 = 0;
613 Int_t ncell2 = 0;
614 if(GetReader()->GetInputEvent()){
615 AliVCluster *cluster1 = (GetReader()->GetInputEvent())->GetCaloCluster(p1->GetCaloLabel(0));
616 ncell1 = cluster1->GetNCells();
617 AliVCluster *cluster2 = (GetReader()->GetInputEvent())->GetCaloCluster(p2->GetCaloLabel(0));
618 ncell2 = cluster2->GetNCells();
619 //printf("e 1: %2.2f - %2.2f, e 2: %2.2f - %2.2f, ncells: %d, %d\n",cluster1->E(),p1->E(),cluster2->E(), p2->E(),ncell1,ncell2);
c8fe2783 620 }
5ae09196 621 for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
622 for(Int_t icell=0; icell<fNCellNCuts; icell++){
623 for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
624
625 if(p1->Pt() > fPtCuts[ipt] && p2->Pt() > fPtCuts[ipt] &&
626 a < fAsymCuts[iasym] &&
627 ncell1 >= fCellNCuts[icell] && ncell2 >= fCellNCuts[icell]) fhRePtNCellAsymCuts[((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym]->Fill(pt,m) ;
628
629 //printf("ipt %d, icell%d, iasym %d, name %s\n",ipt, icell, iasym, fhRePtNCellAsymCuts[((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym]->GetName());
630 }// pid bit cut loop
631 }// icell loop
632 }// pt cut loop
633
634 }// multiple cuts analysis
7e7694bb 635 }// second same event particle
636 }// first cluster
637
638 if(fDoOwnMix){
639 //Fill mixed
640 TList * evMixList=fEventsList[curCentrBin*fNZvertBin*fNrpBin+curZvertBin*fNrpBin+curRPBin] ;
641 Int_t nMixed = evMixList->GetSize() ;
642 for(Int_t ii=0; ii<nMixed; ii++){
643 TClonesArray* ev2= (TClonesArray*) (evMixList->At(ii));
644 Int_t nPhot2=ev2->GetEntriesFast() ;
645 Double_t m = -999;
646 if(GetDebug() > 1) printf("AliAnaPi0::MakeAnalysisFillHistograms() - Mixed event %d photon entries %d\n", ii, nPhot);
647
648 for(Int_t i1=0; i1<nPhot; i1++){
649 AliAODPWG4Particle * p1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i1)) ;
650 TLorentzVector photon1(p1->Px(),p1->Py(),p1->Pz(),p1->E());
651 for(Int_t i2=0; i2<nPhot2; i2++){
652 AliAODPWG4Particle * p2 = (AliAODPWG4Particle*) (ev2->At(i2)) ;
653
654 TLorentzVector photon2(p2->Px(),p2->Py(),p2->Pz(),p2->E());
655 m = (photon1+photon2).M() ;
656 Double_t pt = (photon1 + photon2).Pt();
657 Double_t a = TMath::Abs(p1->E()-p2->E())/(p1->E()+p2->E()) ;
658
659 //Check if opening angle is too large or too small compared to what is expected
660 Double_t angle = photon1.Angle(photon2.Vect());
661 //if(fUseAngleCut && !GetNeutralMesonSelection()->IsAngleInWindow((photon1+photon2).E(),angle)) continue;
662 if(fUseAngleCut && angle < 0.1) continue;
663
664 if(GetDebug() > 2)
665 printf("AliAnaPi0::MakeAnalysisFillHistograms() - Mixed Event: pT: photon1 %2.2f, photon2 %2.2f; Pair: pT %2.2f, mass %2.3f, a %f2.3\n",
666 p1->Pt(), p2->Pt(), pt,m,a);
667 for(Int_t ipid=0; ipid<fNPID; ipid++){
668 if((p1->IsPIDOK(ipid,AliCaloPID::kPhoton)) && (p2->IsPIDOK(ipid,AliCaloPID::kPhoton))){
5ae09196 669 fhMi1 [curCentrBin*fNPID+ipid]->Fill(pt, a,m) ;
670 fhMiInvPt1[curCentrBin*fNPID+ipid]->Fill(1./pt,a,m) ;
7e7694bb 671 if(p1->DistToBad()>0 && p2->DistToBad()>0){
5ae09196 672 fhMi2 [curCentrBin*fNPID+ipid]->Fill(pt, a,m) ;
673 fhMiInvPt2[curCentrBin*fNPID+ipid]->Fill(1./pt,a,m) ;
7e7694bb 674 if(p1->DistToBad()>1 && p2->DistToBad()>1){
5ae09196 675 fhMi3 [curCentrBin*fNPID+ipid]->Fill(pt, a,m) ;
676 fhMiInvPt3[curCentrBin*fNPID+ipid]->Fill(1./pt,a,m) ;
7e7694bb 677 }
678
679 }
680 }
681 }//loop for histograms
682 }// second cluster loop
683 }//first cluster loop
684 }//loop on mixed events
685
686 TClonesArray *currentEvent = new TClonesArray(*GetInputAODBranch());
687 //Ad d current event to buffer and Remove redundant events
688 if(currentEvent->GetEntriesFast()>0){
689 evMixList->AddFirst(currentEvent) ;
690 currentEvent=0 ; //Now list of particles belongs to buffer and it will be deleted with buffer
691 if(evMixList->GetSize()>=fNmaxMixEv)
692 {
693 TClonesArray * tmp = (TClonesArray*) (evMixList->Last()) ;
694 evMixList->RemoveLast() ;
695 delete tmp ;
696 }
697 }
698 else{ //empty event
699 delete currentEvent ;
700 currentEvent=0 ;
477d6cee 701 }
7e7694bb 702 }// DoOwnMix
c8fe2783 703
1c5acb87 704}
705
a5cc4f03 706//________________________________________________________________________
707void AliAnaPi0::ReadHistograms(TList* outputList)
708{
50f39b97 709 // Needed when Terminate is executed in distributed environment
710 // Refill analysis histograms of this class with corresponding histograms in output list.
711
712 // Histograms of this analsys are kept in the same list as other analysis, recover the position of
713 // the first one and then add the next.
714 Int_t index = outputList->IndexOf(outputList->FindObject(GetAddedHistogramsStringToName()+"hRe_cen0_pid0_dist1"));
715
716 if(!fhRe1) fhRe1 = new TH3D*[fNCentrBin*fNPID] ;
717 if(!fhRe2) fhRe2 = new TH3D*[fNCentrBin*fNPID] ;
718 if(!fhRe3) fhRe3 = new TH3D*[fNCentrBin*fNPID] ;
719 if(!fhMi1) fhMi1 = new TH3D*[fNCentrBin*fNPID] ;
720 if(!fhMi2) fhMi2 = new TH3D*[fNCentrBin*fNPID] ;
721 if(!fhMi3) fhMi3 = new TH3D*[fNCentrBin*fNPID] ;
5ae09196 722 if(!fhReInvPt1) fhRe1 = new TH3D*[fNCentrBin*fNPID] ;
723 if(!fhReInvPt2) fhRe2 = new TH3D*[fNCentrBin*fNPID] ;
724 if(!fhReInvPt3) fhRe3 = new TH3D*[fNCentrBin*fNPID] ;
725 if(!fhMiInvPt1) fhMi1 = new TH3D*[fNCentrBin*fNPID] ;
726 if(!fhMiInvPt2) fhMi2 = new TH3D*[fNCentrBin*fNPID] ;
727 if(!fhMiInvPt3) fhMi3 = new TH3D*[fNCentrBin*fNPID] ;
50f39b97 728 if(!fhReMod) fhReMod = new TH3D*[fNModules] ;
729
730 for(Int_t ic=0; ic<fNCentrBin; ic++){
731 for(Int_t ipid=0; ipid<fNPID; ipid++){
732 fhRe1[ic*fNPID+ipid] = (TH3D*) outputList->At(index++);
50f39b97 733 fhRe2[ic*fNPID+ipid] = (TH3D*) outputList->At(index++);
50f39b97 734 fhRe3[ic*fNPID+ipid] = (TH3D*) outputList->At(index++);
5ae09196 735
736 fhReInvPt1[ic*fNPID+ipid] = (TH3D*) outputList->At(index++);
737 fhReInvPt2[ic*fNPID+ipid] = (TH3D*) outputList->At(index++);
738 fhReInvPt3[ic*fNPID+ipid] = (TH3D*) outputList->At(index++);
739
740 fhMi1[ic*fNPID+ipid] = (TH3D*) outputList->At(index++);
741 fhMi2[ic*fNPID+ipid] = (TH3D*) outputList->At(index++);
50f39b97 742 fhMi3[ic*fNPID+ipid] = (TH3D*) outputList->At(index++);
5ae09196 743
744 fhMiInvPt1[ic*fNPID+ipid] = (TH3D*) outputList->At(index++);
745 fhMiInvPt2[ic*fNPID+ipid] = (TH3D*) outputList->At(index++);
746 fhMiInvPt3[ic*fNPID+ipid] = (TH3D*) outputList->At(index++);
50f39b97 747 }
748 }
5ae09196 749 if(fMultiCutAna){
750
751 for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
752 fhRePIDBits[ipid] = (TH2D*) outputList->At(index++);
753 }// ipid loop
754
755 for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
756 for(Int_t icell=0; icell<fNCellNCuts; icell++){
757 for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
758 fhRePtNCellAsymCuts[((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym] = (TH2D*) outputList->At(index++);
759 }// iasym
760 }// icell loop
761 }// ipt loop
762 }// multi cut analysis
50f39b97 763
764 fhEvents = (TH3D *) outputList->At(index++);
765
766 //Histograms filled only if MC data is requested
767 if(IsDataMC() || (GetReader()->GetDataType() == AliCaloTrackReader::kMC) ){
768 fhPrimPt = (TH1D*) outputList->At(index++);
769 fhPrimAccPt = (TH1D*) outputList->At(index++);
770 fhPrimY = (TH1D*) outputList->At(index++);
771 fhPrimAccY = (TH1D*) outputList->At(index++);
772 fhPrimPhi = (TH1D*) outputList->At(index++);
773 fhPrimAccPhi = (TH1D*) outputList->At(index++);
774 }
775
776 for(Int_t imod=0; imod < fNModules; imod++)
777 fhReMod[imod] = (TH3D*) outputList->At(index++);
778
a5cc4f03 779}
780
781
6639984f 782//____________________________________________________________________________________________________________________________________________________
a5cc4f03 783void AliAnaPi0::Terminate(TList* outputList)
6639984f 784{
785 //Do some calculations and plots from the final histograms.
477d6cee 786
fbeaf916 787 printf(" *** %s Terminate:\n", GetName()) ;
50f39b97 788
a5cc4f03 789 //Recover histograms from output histograms list, needed for distributed analysis.
790 ReadHistograms(outputList);
50f39b97 791
2e557d1c 792 if (!fhRe1) {
50f39b97 793 printf("AliAnaPi0::Terminate() - Error: Remote output histograms not imported in AliAnaPi0 object");
794 return;
2e557d1c 795 }
50f39b97 796
a3aebfff 797 printf("AliAnaPi0::Terminate() Mgg Real : %5.3f , RMS : %5.3f \n", fhRe1[0]->GetMean(), fhRe1[0]->GetRMS() ) ;
5ae09196 798
799 const Int_t buffersize = 255;
800
801 char nameIM[buffersize];
802 snprintf(nameIM, buffersize,"AliAnaPi0_%s_cPt",fCalorimeter.Data());
71dd883b 803 TCanvas * cIM = new TCanvas(nameIM, "", 400, 10, 600, 700) ;
6639984f 804 cIM->Divide(2, 2);
50f39b97 805
6639984f 806 cIM->cd(1) ;
807 //gPad->SetLogy();
583c48f1 808 TH1D * hIMAllPt = (TH1D*) fhRe1[0]->ProjectionZ(Form("IMPtAll_%s",fCalorimeter.Data()));
6639984f 809 hIMAllPt->SetLineColor(2);
810 hIMAllPt->SetTitle("No cut on p_{T, #gamma#gamma} ");
811 hIMAllPt->Draw();
812
813 cIM->cd(2) ;
583c48f1 814 TH3F * hRe1Pt5 = (TH3F*)fhRe1[0]->Clone(Form("IMPt5_%s",fCalorimeter.Data()));
6639984f 815 hRe1Pt5->GetXaxis()->SetRangeUser(0,5);
583c48f1 816 TH1D * hIMPt5 = (TH1D*) hRe1Pt5->Project3D(Form("IMPt5_%s_pz",fCalorimeter.Data()));
6639984f 817 hIMPt5->SetLineColor(2);
818 hIMPt5->SetTitle("0 < p_{T, #gamma#gamma} < 5 GeV/c");
819 hIMPt5->Draw();
820
821 cIM->cd(3) ;
583c48f1 822 TH3F * hRe1Pt10 = (TH3F*)fhRe1[0]->Clone(Form("IMPt10_%s",fCalorimeter.Data()));
6639984f 823 hRe1Pt10->GetXaxis()->SetRangeUser(5,10);
583c48f1 824 TH1D * hIMPt10 = (TH1D*) hRe1Pt10->Project3D(Form("IMPt10_%s_pz",fCalorimeter.Data()));
6639984f 825 hIMPt10->SetLineColor(2);
826 hIMPt10->SetTitle("5 < p_{T, #gamma#gamma} < 10 GeV/c");
827 hIMPt10->Draw();
828
829 cIM->cd(4) ;
583c48f1 830 TH3F * hRe1Pt20 = (TH3F*)fhRe1[0]->Clone(Form("IMPt20_%s",fCalorimeter.Data()));
6639984f 831 hRe1Pt20->GetXaxis()->SetRangeUser(10,20);
583c48f1 832 TH1D * hIMPt20 = (TH1D*) hRe1Pt20->Project3D(Form("IMPt20_%s_pz",fCalorimeter.Data()));
6639984f 833 hIMPt20->SetLineColor(2);
834 hIMPt20->SetTitle("10 < p_{T, #gamma#gamma} < 20 GeV/c");
835 hIMPt20->Draw();
836
5ae09196 837 char nameIMF[buffersize];
838 snprintf(nameIMF,buffersize,"AliAnaPi0_%s_Mgg.eps",fCalorimeter.Data());
71dd883b 839 cIM->Print(nameIMF);
6639984f 840
5ae09196 841 char namePt[buffersize];
842 snprintf(namePt,buffersize,"AliAnaPi0_%s_cPt",fCalorimeter.Data());
71dd883b 843 TCanvas * cPt = new TCanvas(namePt, "", 400, 10, 600, 700) ;
6639984f 844 cPt->Divide(2, 2);
845
846 cPt->cd(1) ;
847 //gPad->SetLogy();
848 TH1D * hPt = (TH1D*) fhRe1[0]->Project3D("x");
849 hPt->SetLineColor(2);
850 hPt->SetTitle("No cut on M_{#gamma#gamma} ");
851 hPt->Draw();
852
853 cPt->cd(2) ;
583c48f1 854 TH3F * hRe1IM1 = (TH3F*)fhRe1[0]->Clone(Form("Pt1_%s",fCalorimeter.Data()));
6639984f 855 hRe1IM1->GetZaxis()->SetRangeUser(0.05,0.21);
856 TH1D * hPtIM1 = (TH1D*) hRe1IM1->Project3D("x");
857 hPtIM1->SetLineColor(2);
858 hPtIM1->SetTitle("0.05 < M_{#gamma#gamma} < 0.21 GeV/c^{2}");
859 hPtIM1->Draw();
860
861 cPt->cd(3) ;
583c48f1 862 TH3F * hRe1IM2 = (TH3F*)fhRe1[0]->Clone(Form("Pt2_%s",fCalorimeter.Data()));
6639984f 863 hRe1IM2->GetZaxis()->SetRangeUser(0.09,0.17);
864 TH1D * hPtIM2 = (TH1D*) hRe1IM2->Project3D("x");
865 hPtIM2->SetLineColor(2);
866 hPtIM2->SetTitle("0.09 < M_{#gamma#gamma} < 0.17 GeV/c^{2}");
867 hPtIM2->Draw();
868
869 cPt->cd(4) ;
583c48f1 870 TH3F * hRe1IM3 = (TH3F*)fhRe1[0]->Clone(Form("Pt3_%s",fCalorimeter.Data()));
6639984f 871 hRe1IM3->GetZaxis()->SetRangeUser(0.11,0.15);
872 TH1D * hPtIM3 = (TH1D*) hRe1IM1->Project3D("x");
873 hPtIM3->SetLineColor(2);
874 hPtIM3->SetTitle("0.11 < M_{#gamma#gamma} < 0.15 GeV/c^{2}");
875 hPtIM3->Draw();
876
71dd883b 877 char namePtF[128];
5ae09196 878 snprintf(namePtF,buffersize,"AliAnaPi0_%s_Pt.eps",fCalorimeter.Data());
71dd883b 879 cPt->Print(namePtF);
1c5acb87 880
5ae09196 881 char line[buffersize] ;
882 snprintf(line,buffersize,".!tar -zcf %s_%s.tar.gz *.eps", GetName(),fCalorimeter.Data()) ;
6639984f 883 gROOT->ProcessLine(line);
5ae09196 884 snprintf(line, buffersize,".!rm -fR AliAnaPi0_%s*.eps",fCalorimeter.Data());
6639984f 885 gROOT->ProcessLine(line);
886
71dd883b 887 printf(" AliAnaPi0::Terminate() - !! All the eps files are in %s_%s.tar.gz !!!\n", GetName(), fCalorimeter.Data());
1c5acb87 888
6639984f 889}
c8fe2783 890 //____________________________________________________________________________________________________________________________________________________
891Int_t AliAnaPi0::GetEventIndex(AliAODPWG4Particle * part, Double_t * vert)
892{
893 // retieves the event index and checks the vertex
894 // in the mixed buffer returns -2 if vertex NOK
895 // for normal events returns 0 if vertex OK and -1 if vertex NOK
896
897 Int_t rv = -1 ;
898 if (GetMixedEvent()){
899 TObjArray * pl = 0x0;
900 if (part->GetDetector().Contains("PHOS")) {
901 pl = GetAODPHOS();
902 } else if (part->GetDetector().Contains("EMCAL")) {
903 pl = GetAODEMCAL();
904 } else {
905 AliFatal(Form("%s is an unknown calorimeter", part->GetDetector().Data())) ;
906 }
907 rv = GetMixedEvent()->EventIndexForCaloCluster(part->GetCaloLabel(0)) ;
908 GetMixedEvent()->GetVertexOfEvent(rv)->GetXYZ(vert);
909 if(vert[2] < -fZvtxCut || vert[2] > fZvtxCut)
910 rv = -2 ; //Event can not be used (vertex, centrality,... cuts not fulfilled)
0ae57829 911 } else if(GetReader()->GetDataType()!=AliCaloTrackReader::kMC){
c8fe2783 912 Double_t * tempo = GetReader()->GetVertex() ;
913 vert[0] = tempo[0] ;
914 vert[1] = tempo[1] ;
915 vert[2] = tempo[2] ;
916 if(vert[2] < -fZvtxCut || vert[2] > fZvtxCut)
917 rv = -1 ; //Event can not be used (vertex, centrality,... cuts not fulfilled)
918 else
919 rv = 0 ;
0ae57829 920 }//No MC reader
921 else rv = 0;
922
c8fe2783 923 return rv ;
924}