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