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