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