1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
17 //_________________________________________________________________________
18 // Class to collect two-photon invariant mass distributions for
19 // extracting raw pi0 yield.
20 // Input is produced by AliAnaPhoton (or any other analysis producing output AliAODPWG4Particles),
21 // it will do nothing if executed alone
23 //-- Author: Dmitri Peressounko (RRC "KI")
24 //-- Adapted to PartCorr frame by Lamia Benhabib (SUBATECH)
25 //-- and Gustavo Conesa (INFN-Frascati)
26 //_________________________________________________________________________
29 // --- ROOT system ---
32 //#include "Riostream.h"
36 #include "TClonesArray.h"
37 #include "TObjString.h"
38 #include "TDatabasePDG.h"
40 //---- AliRoot system ----
41 #include "AliAnaPi0.h"
42 #include "AliCaloTrackReader.h"
43 #include "AliCaloPID.h"
45 #include "AliFiducialCut.h"
46 #include "TParticle.h"
47 #include "AliVEvent.h"
48 #include "AliESDCaloCluster.h"
49 #include "AliESDEvent.h"
50 #include "AliAODEvent.h"
51 #include "AliNeutralMesonSelection.h"
52 #include "AliMixedEvent.h"
53 #include "AliAODMCParticle.h"
57 //________________________________________________________________________________________________________________________________________________
58 AliAnaPi0::AliAnaPi0() : AliAnaPartCorrBaseClass(),
59 fDoOwnMix(kFALSE),fNCentrBin(0),//fNZvertBin(0),fNrpBin(0),
60 fNmaxMixEv(0), fCalorimeter(""),
61 fNModules(12), fUseAngleCut(kFALSE), fUseAngleEDepCut(kFALSE),fAngleCut(0), fAngleMaxCut(7.),fEventsList(0x0), fMultiCutAna(kFALSE), fMultiCutAnaSim(kFALSE),
62 fNPtCuts(0),fNAsymCuts(0), fNCellNCuts(0),fNPIDBits(0), fMakeInvPtPlots(kFALSE), fSameSM(kFALSE),
63 fUseTrackMultBins(kFALSE),fUsePhotonMultBins(kFALSE),fUseAverCellEBins(kFALSE), fUseAverClusterEBins(kFALSE),
64 fUseAverClusterEDenBins(0), //fUseAverClusterPairRBins(0), fUseAverClusterPairRWeightBins(0), fUseEMaxBins(0),
65 fFillBadDistHisto(kFALSE),
66 fhAverTotECluster(0), fhAverTotECell(0), fhAverTotECellvsCluster(0),
67 fhEDensityCluster(0), fhEDensityCell(0), fhEDensityCellvsCluster(0),
68 //fhClusterPairDist(0), fhClusterPairDistWeight(0), fhAverClusterPairDist(0), fhAverClusterPairDistWeight(0),
69 //fhAverClusterPairDistvsAverE(0), fhAverClusterPairDistWeightvsAverE(0),fhAverClusterPairDistvsN(0), fhAverClusterPairDistWeightvsN(0),
70 //fhMaxEvsClustMult(0), fhMaxEvsClustEDen(0),
71 fhReMod(0x0), fhReDiffMod(0x0), fhMiMod(0x0), fhMiDiffMod(0x0),
72 fhReConv(0x0), fhMiConv(0x0), fhReConv2(0x0), fhMiConv2(0x0),
73 fhRe1(0x0), fhMi1(0x0), fhRe2(0x0), fhMi2(0x0), fhRe3(0x0), fhMi3(0x0),
74 fhReInvPt1(0x0), fhMiInvPt1(0x0), fhReInvPt2(0x0), fhMiInvPt2(0x0), fhReInvPt3(0x0), fhMiInvPt3(0x0),
75 fhRePtNCellAsymCuts(0x0), fhRePtNCellAsymCutsSM0(0x0), fhRePtNCellAsymCutsSM1(0x0), fhRePtNCellAsymCutsSM2(0x0), fhRePtNCellAsymCutsSM3(0x0), fhMiPtNCellAsymCuts(0x0),
76 fhRePIDBits(0x0),fhRePtMult(0x0), fhRePtAsym(0x0), fhRePtAsymPi0(0x0),fhRePtAsymEta(0x0),
77 fhEvents(0x0), fhRealOpeningAngle(0x0),fhRealCosOpeningAngle(0x0), fhMixedOpeningAngle(0x0),fhMixedCosOpeningAngle(0x0),
78 fhPrimPi0Pt(0x0), fhPrimPi0AccPt(0x0), fhPrimPi0Y(0x0), fhPrimPi0AccY(0x0), fhPrimPi0Phi(0x0), fhPrimPi0AccPhi(0x0),
79 fhPrimPi0OpeningAngle(0x0),fhPrimPi0CosOpeningAngle(0x0),
80 fhPrimEtaPt(0x0), fhPrimEtaAccPt(0x0), fhPrimEtaY(0x0), fhPrimEtaAccY(0x0), fhPrimEtaPhi(0x0), fhPrimEtaAccPhi(0x0),
81 fhMCOrgMass(),fhMCOrgAsym(), fhMCOrgDeltaEta(),fhMCOrgDeltaPhi(),
82 fhMCPi0MassPtRec(), fhMCPi0MassPtTrue(), fhMCPi0PtTruePtRec(), fhMCEtaMassPtRec(), fhMCEtaMassPtTrue(), fhMCEtaPtTruePtRec()
89 //________________________________________________________________________________________________________________________________________________
90 AliAnaPi0::~AliAnaPi0() {
91 // Remove event containers
93 if(fDoOwnMix && fEventsList){
94 for(Int_t ic=0; ic<fNCentrBin; ic++){
95 for(Int_t iz=0; iz<GetNZvertBin(); iz++){
96 for(Int_t irp=0; irp<GetNRPBin(); irp++){
97 fEventsList[ic*GetNZvertBin()*GetNRPBin()+iz*GetNRPBin()+irp]->Delete() ;
98 delete fEventsList[ic*GetNZvertBin()*GetNRPBin()+iz*GetNRPBin()+irp] ;
102 delete[] fEventsList;
108 //________________________________________________________________________________________________________________________________________________
109 void AliAnaPi0::InitParameters()
111 //Init parameters when first called the analysis
112 //Set default parameters
113 SetInputAODName("PWG4Particle");
115 AddToHistogramsName("AnaPi0_");
116 fNModules = 12; // set maximum to maximum number of EMCAL modules
122 fCalorimeter = "PHOS";
123 fUseAngleCut = kFALSE;
124 fUseAngleEDepCut = kFALSE;
126 fAngleMaxCut = TMath::Pi();
128 fMultiCutAna = kFALSE;
131 fPtCuts[0] = 0.; fPtCuts[1] = 0.3; fPtCuts[2] = 0.5;
132 for(Int_t i = fNPtCuts; i < 10; i++)fPtCuts[i] = 0.;
135 fAsymCuts[0] = 1.; fAsymCuts[1] = 0.8; fAsymCuts[2] = 0.6; fAsymCuts[3] = 0.1;
136 for(Int_t i = fNAsymCuts; i < 10; i++)fAsymCuts[i] = 0.;
139 fCellNCuts[0] = 0; fCellNCuts[1] = 1; fCellNCuts[2] = 2;
140 for(Int_t i = fNCellNCuts; i < 10; i++)fCellNCuts[i] = 0;
143 fPIDBits[0] = 0; fPIDBits[1] = 2; // fPIDBits[2] = 4; fPIDBits[3] = 6;// check, no cut, dispersion, neutral, dispersion&&neutral
144 for(Int_t i = fNPIDBits; i < 10; i++)fPIDBits[i] = 0;
149 //________________________________________________________________________________________________________________________________________________
150 TObjString * AliAnaPi0::GetAnalysisCuts()
152 //Save parameters used for analysis
153 TString parList ; //this will be list of parameters used for this analysis.
154 const Int_t buffersize = 255;
155 char onePar[buffersize] ;
156 snprintf(onePar,buffersize,"--- AliAnaPi0 ---\n") ;
158 snprintf(onePar,buffersize,"Number of bins in Centrality: %d \n",fNCentrBin) ;
160 snprintf(onePar,buffersize,"Number of bins in Z vert. pos: %d \n",GetNZvertBin()) ;
162 snprintf(onePar,buffersize,"Number of bins in Reac. Plain: %d \n",GetNRPBin()) ;
164 snprintf(onePar,buffersize,"Depth of event buffer: %d \n",fNmaxMixEv) ;
166 snprintf(onePar,buffersize,"Pair in same Module: %d ; TrackMult as centrality: %d; PhotonMult as centrality: %d; cluster E as centrality: %d; cell as centrality: %d; Fill InvPt histos %d\n",
167 fSameSM, fUseTrackMultBins, fUsePhotonMultBins, fUseAverClusterEBins, fUseAverCellEBins, fMakeInvPtPlots) ;
169 snprintf(onePar,buffersize,"Select pairs with their angle: %d, edep %d, min angle %2.3f, max angle %2.3f,\n",fUseAngleCut, fUseAngleEDepCut,fAngleCut,fAngleMaxCut) ;
171 snprintf(onePar,buffersize," Asymmetry cuts: n = %d, asymmetry < ",fNAsymCuts) ;
172 for(Int_t i = 0; i < fNAsymCuts; i++) snprintf(onePar,buffersize,"%s %2.2f;",onePar,fAsymCuts[i]);
174 snprintf(onePar,buffersize," PID selection bits: n = %d, PID bit =\n",fNPIDBits) ;
175 for(Int_t i = 0; i < fNPIDBits; i++) snprintf(onePar,buffersize,"%s %d;",onePar,fPIDBits[i]);
177 snprintf(onePar,buffersize,"Cuts: \n") ;
179 snprintf(onePar,buffersize,"Z vertex position: -%f < z < %f \n",GetZvertexCut(),GetZvertexCut()) ;
181 snprintf(onePar,buffersize,"Calorimeter: %s \n",fCalorimeter.Data()) ;
183 snprintf(onePar,buffersize,"Number of modules: %d \n",fNModules) ;
186 snprintf(onePar, buffersize," pT cuts: n = %d, pt > ",fNPtCuts) ;
187 for(Int_t i = 0; i < fNPtCuts; i++) snprintf(onePar,buffersize,"%s %2.2f;",onePar,fPtCuts[i]);
189 snprintf(onePar,buffersize, " N cell in cluster cuts: n = %d, nCell > ",fNCellNCuts) ;
190 for(Int_t i = 0; i < fNCellNCuts; i++) snprintf(onePar,buffersize,"%s %d;",onePar,fCellNCuts[i]);
194 return new TObjString(parList) ;
197 //________________________________________________________________________________________________________________________________________________
198 TList * AliAnaPi0::GetCreateOutputObjects()
200 // Create histograms to be saved in output file and
201 // store them in fOutputContainer
203 //create event containers
204 fEventsList = new TList*[fNCentrBin*GetNZvertBin()*GetNRPBin()] ;
206 for(Int_t ic=0; ic<fNCentrBin; ic++){
207 for(Int_t iz=0; iz<GetNZvertBin(); iz++){
208 for(Int_t irp=0; irp<GetNRPBin(); irp++){
209 fEventsList[ic*GetNZvertBin()*GetNRPBin()+iz*GetNRPBin()+irp] = new TList() ;
210 fEventsList[ic*GetNZvertBin()*GetNRPBin()+iz*GetNRPBin()+irp]->SetOwner(kFALSE);
215 TList * outputContainer = new TList() ;
216 outputContainer->SetName(GetName());
218 fhReMod = new TH2D*[fNModules] ;
219 fhReDiffMod = new TH2D*[fNModules+3] ;
221 fhMiMod = new TH2D*[fNModules] ;
222 fhMiDiffMod = new TH2D*[fNModules+3] ;
224 fhRe1 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
225 fhMi1 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
226 if(fFillBadDistHisto){
227 fhRe2 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
228 fhRe3 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
229 fhMi2 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
230 fhMi3 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
232 if(fMakeInvPtPlots) {
233 fhReInvPt1 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
234 fhMiInvPt1 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
235 if(fFillBadDistHisto){
236 fhReInvPt2 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
237 fhReInvPt3 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
238 fhMiInvPt2 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
239 fhMiInvPt3 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
243 const Int_t buffersize = 255;
244 char key[buffersize] ;
245 char title[buffersize] ;
247 Int_t nptbins = GetHistoPtBins();
248 Int_t nphibins = GetHistoPhiBins();
249 Int_t netabins = GetHistoEtaBins();
250 Float_t ptmax = GetHistoPtMax();
251 Float_t phimax = GetHistoPhiMax();
252 Float_t etamax = GetHistoEtaMax();
253 Float_t ptmin = GetHistoPtMin();
254 Float_t phimin = GetHistoPhiMin();
255 Float_t etamin = GetHistoEtaMin();
257 Int_t nmassbins = GetHistoMassBins();
258 Int_t nasymbins = GetHistoAsymmetryBins();
259 Float_t massmax = GetHistoMassMax();
260 Float_t asymmax = GetHistoAsymmetryMax();
261 Float_t massmin = GetHistoMassMin();
262 Float_t asymmin = GetHistoAsymmetryMin();
263 Int_t ntrmbins = GetHistoTrackMultiplicityBins();
264 Int_t ntrmmax = GetHistoTrackMultiplicityMax();
265 Int_t ntrmmin = GetHistoTrackMultiplicityMin();
267 fhAverTotECluster = new TH1F("hAverTotECluster","hAverTotECluster",200,0,50) ;
268 fhAverTotECluster->SetXTitle("E_{cluster, aver. SM} (GeV)");
269 outputContainer->Add(fhAverTotECluster) ;
271 fhAverTotECell = new TH1F("hAverTotECell","hAverTotECell",200,0,50) ;
272 fhAverTotECell->SetXTitle("E_{cell, aver. SM} (GeV)");
273 outputContainer->Add(fhAverTotECell) ;
275 fhAverTotECellvsCluster = new TH2F("hAverTotECellvsCluster","hAverTotECellvsCluster",200,0,50,200,0,50) ;
276 fhAverTotECellvsCluster->SetYTitle("E_{cell, aver. SM} (GeV)");
277 fhAverTotECellvsCluster->SetXTitle("E_{cluster, aver. SM} (GeV)");
278 outputContainer->Add(fhAverTotECellvsCluster) ;
280 fhEDensityCluster = new TH1F("hEDensityCluster","hEDensityCluster",200,0,50) ;
281 fhEDensityCluster->SetXTitle("#Sigma E_{cluster} / N_{cluster} (GeV)");
282 outputContainer->Add(fhEDensityCluster) ;
284 fhEDensityCell = new TH1F("hEDensityCell","hEDensityCell",200,0,50) ;
285 fhEDensityCell->SetXTitle("#Sigma E_{cell} / N_{cell} (GeV)");
286 outputContainer->Add(fhEDensityCell) ;
288 fhEDensityCellvsCluster = new TH2F("hEDensityCellvsCluster","hEDensityCellvsCluster",200,0,50,200,0,50) ;
289 fhEDensityCellvsCluster->SetYTitle("#Sigma E_{cell} / N_{cell} (GeV)");
290 fhEDensityCellvsCluster->SetXTitle("#Sigma E_{cluster} / N_{cluster} (GeV)");
291 outputContainer->Add(fhEDensityCellvsCluster) ;
293 // fhClusterPairDist = new TH1F("hClusterPairDist","Distance between clusters",250,0,750) ;
294 // fhClusterPairDist->SetXTitle("#sqrt{(x_{1}-x_{2})^2+(z_{1}-z_{2})^2} (cm)");
295 // outputContainer->Add(fhClusterPairDist) ;
297 // fhClusterPairDistWeight = new TH1F("hClusterPairDistWeighted","Distance between clusters, weighted by pair energy",200,0,400) ;
298 // fhClusterPairDistWeight->SetXTitle("#sqrt{(x_{1}E_{1}-x_{2}E_{2})^{2}+(z_{1}E_{1}-z_{2}E_{2})^{2}}/ (E_{1}+E_{2}) (cm)");
299 // outputContainer->Add(fhClusterPairDistWeight) ;
301 // fhAverClusterPairDist = new TH1F("hAverClusterPairDist","Average distance between clusters",250,0,750) ;
302 // fhAverClusterPairDist->SetXTitle("#Sigma (#sqrt{(x_{1}-x_{2})^{2}+(z_{1}-z_{2})^{2}}) / N_{pairs} (cm)");
303 // outputContainer->Add(fhAverClusterPairDist) ;
305 // fhAverClusterPairDistWeight = new TH1F("hAverClusterPairDistWeighted","Average distance between clusters, weighted by pair energy",200,0,400) ;
306 // fhAverClusterPairDistWeight->SetXTitle("#Sigma (#sqrt{(x_{1}E_{1}-x_{2}E_{2})^{2}+(z_{1}E_{1}-z_{2}E_{2})^{2}}/ (E_{1}+E_{2})) / N_{pairs} (cm)");
307 // outputContainer->Add(fhAverClusterPairDistWeight) ;
309 // fhAverClusterPairDistvsAverE = new TH2F("hAverClusterPairDistvsAverE","Average distance between clusters",250,0,750,200,0,50) ;
310 // fhAverClusterPairDistvsAverE->SetXTitle("#Sigma (#sqrt{(x_{1}-x_{2})^{2}+(z_{1}-z_{2})^{2}}) / N_{pairs} (cm)");
311 // fhAverClusterPairDistvsAverE->SetYTitle("#Sigma E_{cluster} / N_{cluster} (GeV)");
312 // outputContainer->Add(fhAverClusterPairDistvsAverE) ;
314 // fhAverClusterPairDistWeightvsAverE = new TH2F("hAverClusterPairDistWeightedvsAverE","Average distance between clusters, weighted by pair energy",200,0,400,200,0,50) ;
315 // fhAverClusterPairDistWeightvsAverE->SetXTitle("#Sigma (#sqrt{(x_{1}E_{1}-x_{2}E_{2})^2+(z_{1}E_{1}-z_{2}E_{2})^2}/ (E_{1}+E_{2})) / N_{pairs} (cm/GeV)");
316 // fhAverClusterPairDistWeightvsAverE->SetYTitle("#Sigma E_{cluster} / N_{cluster} (GeV)");
317 // outputContainer->Add(fhAverClusterPairDistWeightvsAverE) ;
319 // fhAverClusterPairDistvsN = new TH2F("hAverClusterPairDistvsN","Average distance between clusters",250,0,750,200,0,50) ;
320 // fhAverClusterPairDistvsN->SetXTitle("#Sigma (#sqrt{(x_{1}-x_{2})^{2}+(z_{1}-z_{2})^{2}}) / N_{pairs} (cm)");
321 // fhAverClusterPairDistvsN->SetYTitle("N_{cluster}");
322 // outputContainer->Add(fhAverClusterPairDistvsN) ;
324 // fhAverClusterPairDistWeightvsN = new TH2F("hAverClusterPairDistWeightedvsN","Average distance between clusters, weighted by pair energy",200,0,400,200,0,50) ;
325 // fhAverClusterPairDistWeightvsN->SetXTitle("#Sigma (#sqrt{(x_{1}E_{1}-x_{2}E_{2})^{2}+(z_{1}E_{1}-z_{2}E_{2})^{2}}/ (E_{1}+E_{2})) / N_{pairs} (cm)");
326 // fhAverClusterPairDistWeightvsN->SetYTitle("N_{cluster}");
327 // outputContainer->Add(fhAverClusterPairDistWeightvsN) ;
329 // fhMaxEvsClustMult = new TH2F("hMaxEvsClustMult","",nptbins,ptmin,ptmax,50,0,50) ;
330 // fhMaxEvsClustMult->SetXTitle("E_{max}");
331 // fhMaxEvsClustMult->SetYTitle("N_{cluster}");
332 // outputContainer->Add(fhMaxEvsClustMult) ;
334 // fhMaxEvsClustEDen = new TH2F("hMaxEvsClustEDen","",nptbins,ptmin,ptmax,200,0,50) ;
335 // fhMaxEvsClustEDen->SetXTitle("E_{max}");
336 // fhMaxEvsClustEDen->SetYTitle("#Sigma E_{cluster} / N_{cluster} (GeV)");
337 // outputContainer->Add(fhMaxEvsClustEDen) ;
339 fhReConv = new TH2D("hReConv","Real Pair with one recombined conversion ",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
340 fhReConv->SetXTitle("p_{T} (GeV/c)");
341 fhReConv->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
342 outputContainer->Add(fhReConv) ;
344 fhReConv2 = new TH2D("hReConv2","Real Pair with 2 recombined conversion ",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
345 fhReConv2->SetXTitle("p_{T} (GeV/c)");
346 fhReConv2->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
347 outputContainer->Add(fhReConv2) ;
350 fhMiConv = new TH2D("hMiConv","Mixed Pair with one recombined conversion ",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
351 fhMiConv->SetXTitle("p_{T} (GeV/c)");
352 fhMiConv->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
353 outputContainer->Add(fhMiConv) ;
355 fhMiConv2 = new TH2D("hMiConv2","Mixed Pair with 2 recombined conversion ",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
356 fhMiConv2->SetXTitle("p_{T} (GeV/c)");
357 fhMiConv2->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
358 outputContainer->Add(fhMiConv2) ;
361 for(Int_t ic=0; ic<fNCentrBin; ic++){
362 for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
363 for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
364 Int_t index = ((ic*fNPIDBits)+ipid)*fNAsymCuts + iasym;
365 //printf("cen %d, pid %d, asy %d, Index %d\n",ic,ipid,iasym,index);
366 //Distance to bad module 1
367 snprintf(key, buffersize,"hRe_cen%d_pidbit%d_asy%d_dist1",ic,ipid,iasym) ;
368 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 1",
369 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
370 fhRe1[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
371 fhRe1[index]->SetXTitle("p_{T} (GeV/c)");
372 fhRe1[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
373 //printf("name: %s\n ",fhRe1[index]->GetName());
374 outputContainer->Add(fhRe1[index]) ;
376 if(fFillBadDistHisto){
377 //Distance to bad module 2
378 snprintf(key, buffersize,"hRe_cen%d_pidbit%d_asy%d_dist2",ic,ipid,iasym) ;
379 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 2",
380 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
381 fhRe2[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
382 fhRe2[index]->SetXTitle("p_{T} (GeV/c)");
383 fhRe2[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
384 outputContainer->Add(fhRe2[index]) ;
386 //Distance to bad module 3
387 snprintf(key, buffersize,"hRe_cen%d_pidbit%d_asy%d_dist3",ic,ipid,iasym) ;
388 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 3",
389 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
390 fhRe3[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
391 fhRe3[index]->SetXTitle("p_{T} (GeV/c)");
392 fhRe3[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
393 outputContainer->Add(fhRe3[index]) ;
398 //Distance to bad module 1
399 snprintf(key, buffersize,"hReInvPt_cen%d_pidbit%d_asy%d_dist1",ic,ipid,iasym) ;
400 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 1",
401 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
402 fhReInvPt1[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
403 fhReInvPt1[index]->SetXTitle("p_{T} (GeV/c)");
404 fhReInvPt1[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
405 outputContainer->Add(fhReInvPt1[index]) ;
407 if(fFillBadDistHisto){
408 //Distance to bad module 2
409 snprintf(key, buffersize,"hReInvPt_cen%d_pidbit%d_asy%d_dist2",ic,ipid,iasym) ;
410 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 2",
411 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
412 fhReInvPt2[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
413 fhReInvPt2[index]->SetXTitle("p_{T} (GeV/c)");
414 fhReInvPt2[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
415 outputContainer->Add(fhReInvPt2[index]) ;
417 //Distance to bad module 3
418 snprintf(key, buffersize,"hReInvPt_cen%d_pidbit%d_asy%d_dist3",ic,ipid,iasym) ;
419 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 3",
420 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
421 fhReInvPt3[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
422 fhReInvPt3[index]->SetXTitle("p_{T} (GeV/c)");
423 fhReInvPt3[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
424 outputContainer->Add(fhReInvPt3[index]) ;
428 //Distance to bad module 1
429 snprintf(key, buffersize,"hMi_cen%d_pidbit%d_asy%d_dist1",ic,ipid,iasym) ;
430 snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 1",
431 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
432 fhMi1[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
433 fhMi1[index]->SetXTitle("p_{T} (GeV/c)");
434 fhMi1[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
435 outputContainer->Add(fhMi1[index]) ;
436 if(fFillBadDistHisto){
437 //Distance to bad module 2
438 snprintf(key, buffersize,"hMi_cen%d_pidbit%d_asy%d_dist2",ic,ipid,iasym) ;
439 snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 2",
440 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
441 fhMi2[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
442 fhMi2[index]->SetXTitle("p_{T} (GeV/c)");
443 fhMi2[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
444 outputContainer->Add(fhMi2[index]) ;
446 //Distance to bad module 3
447 snprintf(key, buffersize,"hMi_cen%d_pidbit%d_asy%d_dist3",ic,ipid,iasym) ;
448 snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 3",
449 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
450 fhMi3[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
451 fhMi3[index]->SetXTitle("p_{T} (GeV/c)");
452 fhMi3[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
453 outputContainer->Add(fhMi3[index]) ;
457 //Distance to bad module 1
458 snprintf(key, buffersize,"hMiInvPt_cen%d_pidbit%d_asy%d_dist1",ic,ipid,iasym) ;
459 snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 1",
460 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
461 fhMiInvPt1[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
462 fhMiInvPt1[index]->SetXTitle("p_{T} (GeV/c)");
463 fhMiInvPt1[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
464 outputContainer->Add(fhMiInvPt1[index]) ;
465 if(fFillBadDistHisto){
466 //Distance to bad module 2
467 snprintf(key, buffersize,"hMiInvPt_cen%d_pidbit%d_asy%d_dist2",ic,ipid,iasym) ;
468 snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 2",
469 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
470 fhMiInvPt2[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
471 fhMiInvPt2[index]->SetXTitle("p_{T} (GeV/c)");
472 fhMiInvPt2[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
473 outputContainer->Add(fhMiInvPt2[index]) ;
475 //Distance to bad module 3
476 snprintf(key, buffersize,"hMiInvPt_cen%d_pidbit%d_asy%d_dist3",ic,ipid,iasym) ;
477 snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f,dist bad 3",
478 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
479 fhMiInvPt3[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
480 fhMiInvPt3[index]->SetXTitle("p_{T} (GeV/c)");
481 fhMiInvPt3[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
482 outputContainer->Add(fhMiInvPt3[index]) ;
490 fhRePtAsym = new TH2D("hRePtAsym","Asymmetry vs pt, for pairs",nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax) ;
491 fhRePtAsym->SetXTitle("p_{T} (GeV/c)");
492 fhRePtAsym->SetYTitle("Asymmetry");
493 outputContainer->Add(fhRePtAsym);
495 fhRePtAsymPi0 = new TH2D("hRePtAsymPi0","Asymmetry vs pt, for pairs close to #pi^{0} mass",nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax) ;
496 fhRePtAsymPi0->SetXTitle("p_{T} (GeV/c)");
497 fhRePtAsymPi0->SetYTitle("Asymmetry");
498 outputContainer->Add(fhRePtAsymPi0);
500 fhRePtAsymEta = new TH2D("hRePtAsymEta","Asymmetry vs pt, for pairs close to #eta mass",nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax) ;
501 fhRePtAsymEta->SetXTitle("p_{T} (GeV/c)");
502 fhRePtAsymEta->SetYTitle("Asymmetry");
503 outputContainer->Add(fhRePtAsymEta);
507 fhRePIDBits = new TH2D*[fNPIDBits];
508 for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
509 snprintf(key, buffersize,"hRe_pidbit%d",ipid) ;
510 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for PIDBit=%d",fPIDBits[ipid]) ;
511 fhRePIDBits[ipid] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
512 fhRePIDBits[ipid]->SetXTitle("p_{T} (GeV/c)");
513 fhRePIDBits[ipid]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
514 outputContainer->Add(fhRePIDBits[ipid]) ;
517 fhRePtNCellAsymCuts = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
518 fhRePtNCellAsymCutsSM0 = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
519 fhRePtNCellAsymCutsSM1 = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
520 fhRePtNCellAsymCutsSM2 = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
521 fhRePtNCellAsymCutsSM3 = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
522 fhMiPtNCellAsymCuts = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
523 for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
524 for(Int_t icell=0; icell<fNCellNCuts; icell++){
525 for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
526 snprintf(key, buffersize,"hRe_pt%d_cell%d_asym%d",ipt,icell,iasym) ;
527 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for pt >%2.2f, ncell>%d and asym >%1.2f ",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]) ;
528 Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
529 //printf("ipt %d, icell %d, iassym %d, index %d\n",ipt, icell, iasym, index);
530 fhRePtNCellAsymCuts[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
531 fhRePtNCellAsymCuts[index]->SetXTitle("p_{T} (GeV/c)");
532 fhRePtNCellAsymCuts[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
533 outputContainer->Add(fhRePtNCellAsymCuts[index]) ;
535 snprintf(key, buffersize,"hRe_pt%d_cell%d_asym%d_SM0",ipt,icell,iasym) ;
536 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for pt >%2.2f, ncell>%d and asym >%1.2f, SM 0 ",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]) ;
537 fhRePtNCellAsymCutsSM0[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
538 fhRePtNCellAsymCutsSM0[index]->SetXTitle("p_{T} (GeV/c)");
539 fhRePtNCellAsymCutsSM0[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
540 outputContainer->Add(fhRePtNCellAsymCutsSM0[index]) ;
542 snprintf(key, buffersize,"hRe_pt%d_cell%d_asym%d_SM1",ipt,icell,iasym) ;
543 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for pt >%2.2f, ncell>%d and asym >%1.2f, SM 1 ",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]) ;
544 fhRePtNCellAsymCutsSM1[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
545 fhRePtNCellAsymCutsSM1[index]->SetXTitle("p_{T} (GeV/c)");
546 fhRePtNCellAsymCutsSM1[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
547 outputContainer->Add(fhRePtNCellAsymCutsSM1[index]) ;
549 snprintf(key, buffersize,"hRe_pt%d_cell%d_asym%d_SM2",ipt,icell,iasym) ;
550 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for pt >%2.2f, ncell>%d and asym >%1.2f, SM 2 ",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]) ;
551 fhRePtNCellAsymCutsSM2[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
552 fhRePtNCellAsymCutsSM2[index]->SetXTitle("p_{T} (GeV/c)");
553 fhRePtNCellAsymCutsSM2[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
554 outputContainer->Add(fhRePtNCellAsymCutsSM2[index]) ;
556 snprintf(key, buffersize,"hRe_pt%d_cell%d_asym%d_SM3",ipt,icell,iasym) ;
557 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for pt >%2.2f, ncell>%d and asym >%1.2f, SM 3 ",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]) ;
558 fhRePtNCellAsymCutsSM3[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
559 fhRePtNCellAsymCutsSM3[index]->SetXTitle("p_{T} (GeV/c)");
560 fhRePtNCellAsymCutsSM3[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
561 outputContainer->Add(fhRePtNCellAsymCutsSM3[index]) ;
563 snprintf(key, buffersize,"hMi_pt%d_cell%d_asym%d",ipt,icell,iasym) ;
564 snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for pt >%2.2f, ncell>%d and asym >%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]) ;
565 fhMiPtNCellAsymCuts[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
566 fhMiPtNCellAsymCuts[index]->SetXTitle("p_{T} (GeV/c)");
567 fhMiPtNCellAsymCuts[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
568 outputContainer->Add(fhMiPtNCellAsymCuts[index]) ;
574 fhRePtMult = new TH3D*[fNAsymCuts] ;
575 for(Int_t iasym = 0; iasym<fNAsymCuts; iasym++){
576 fhRePtMult[iasym] = new TH3D(Form("hRePtMult_asym%d",iasym),Form("(p_{T},C,M)_{#gamma#gamma}, A<%1.2f",fAsymCuts[iasym]),
577 nptbins,ptmin,ptmax,ntrmbins,ntrmmin,ntrmmax,nmassbins,massmin,massmax);
578 fhRePtMult[iasym]->SetXTitle("p_{T} (GeV/c)");
579 fhRePtMult[iasym]->SetYTitle("Track multiplicity");
580 fhRePtMult[iasym]->SetZTitle("m_{#gamma,#gamma} (GeV/c^{2})");
581 outputContainer->Add(fhRePtMult[iasym]) ;
584 }// multi cuts analysis
586 fhEvents=new TH3D("hEvents","Number of events",fNCentrBin,0.,1.*fNCentrBin,
587 GetNZvertBin(),0.,1.*GetNZvertBin(),GetNRPBin(),0.,1.*GetNRPBin()) ;
588 outputContainer->Add(fhEvents) ;
590 fhRealOpeningAngle = new TH2D
591 ("hRealOpeningAngle","Angle between all #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,300,0,TMath::Pi());
592 fhRealOpeningAngle->SetYTitle("#theta(rad)");
593 fhRealOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
594 outputContainer->Add(fhRealOpeningAngle) ;
596 fhRealCosOpeningAngle = new TH2D
597 ("hRealCosOpeningAngle","Cosinus of angle between all #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,100,0,1);
598 fhRealCosOpeningAngle->SetYTitle("cos (#theta) ");
599 fhRealCosOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
600 outputContainer->Add(fhRealCosOpeningAngle) ;
604 fhMixedOpeningAngle = new TH2D
605 ("hMixedOpeningAngle","Angle between all #gamma pair vs E_{#pi^{0}}, Mixed pairs",nptbins,ptmin,ptmax,300,0,TMath::Pi());
606 fhMixedOpeningAngle->SetYTitle("#theta(rad)");
607 fhMixedOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
608 outputContainer->Add(fhMixedOpeningAngle) ;
610 fhMixedCosOpeningAngle = new TH2D
611 ("hMixedCosOpeningAngle","Cosinus of angle between all #gamma pair vs E_{#pi^{0}}, Mixed pairs",nptbins,ptmin,ptmax,100,0,1);
612 fhMixedCosOpeningAngle->SetYTitle("cos (#theta) ");
613 fhMixedCosOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
614 outputContainer->Add(fhMixedCosOpeningAngle) ;
618 //Histograms filled only if MC data is requested
621 fhPrimPi0Pt = new TH1D("hPrimPi0Pt","Primary pi0 pt",nptbins,ptmin,ptmax) ;
622 fhPrimPi0AccPt = new TH1D("hPrimPi0AccPt","Primary pi0 pt with both photons in acceptance",nptbins,ptmin,ptmax) ;
623 fhPrimPi0Pt ->SetXTitle("p_{T} (GeV/c)");
624 fhPrimPi0AccPt->SetXTitle("p_{T} (GeV/c)");
625 outputContainer->Add(fhPrimPi0Pt) ;
626 outputContainer->Add(fhPrimPi0AccPt) ;
628 fhPrimPi0Y = new TH2D("hPrimPi0Rapidity","Rapidity of primary pi0",nptbins,ptmin,ptmax, netabins,etamin,etamax) ;
629 fhPrimPi0Y ->SetYTitle("Rapidity");
630 fhPrimPi0Y ->SetXTitle("p_{T} (GeV/c)");
631 outputContainer->Add(fhPrimPi0Y) ;
633 fhPrimPi0AccY = new TH2D("hPrimPi0AccRapidity","Rapidity of primary pi0",nptbins,ptmin,ptmax, netabins,etamin,etamax) ;
634 fhPrimPi0AccY->SetYTitle("Rapidity");
635 fhPrimPi0AccY->SetXTitle("p_{T} (GeV/c)");
636 outputContainer->Add(fhPrimPi0AccY) ;
638 fhPrimPi0Phi = new TH2D("hPrimPi0Phi","Azimuthal of primary pi0",nptbins,ptmin,ptmax, nphibins,phimin*TMath::RadToDeg(),phimax*TMath::RadToDeg()) ;
639 fhPrimPi0Phi->SetYTitle("#phi (deg)");
640 fhPrimPi0Phi->SetXTitle("p_{T} (GeV/c)");
641 outputContainer->Add(fhPrimPi0Phi) ;
643 fhPrimPi0AccPhi = new TH2D("hPrimPi0AccPhi","Azimuthal of primary pi0 with accepted daughters",nptbins,ptmin,ptmax, nphibins,phimin*TMath::RadToDeg(),phimax*TMath::RadToDeg()) ;
644 fhPrimPi0AccPhi->SetYTitle("#phi (deg)");
645 fhPrimPi0AccPhi->SetXTitle("p_{T} (GeV/c)");
646 outputContainer->Add(fhPrimPi0AccPhi) ;
649 fhPrimEtaPt = new TH1D("hPrimEtaPt","Primary eta pt",nptbins,ptmin,ptmax) ;
650 fhPrimEtaAccPt = new TH1D("hPrimEtaAccPt","Primary eta pt with both photons in acceptance",nptbins,ptmin,ptmax) ;
651 fhPrimEtaPt ->SetXTitle("p_{T} (GeV/c)");
652 fhPrimEtaAccPt->SetXTitle("p_{T} (GeV/c)");
653 outputContainer->Add(fhPrimEtaPt) ;
654 outputContainer->Add(fhPrimEtaAccPt) ;
656 fhPrimEtaY = new TH2D("hPrimEtaRapidity","Rapidity of primary eta",nptbins,ptmin,ptmax,netabins,etamin,etamax) ;
657 fhPrimEtaY->SetYTitle("Rapidity");
658 fhPrimEtaY->SetXTitle("p_{T} (GeV/c)");
659 outputContainer->Add(fhPrimEtaY) ;
661 fhPrimEtaAccY = new TH2D("hPrimEtaAccRapidity","Rapidity of primary eta",nptbins,ptmin,ptmax, netabins,etamin,etamax) ;
662 fhPrimEtaAccY->SetYTitle("Rapidity");
663 fhPrimEtaAccY->SetXTitle("p_{T} (GeV/c)");
664 outputContainer->Add(fhPrimEtaAccY) ;
666 fhPrimEtaPhi = new TH2D("hPrimEtaPhi","Azimuthal of primary eta",nptbins,ptmin,ptmax, nphibins,phimin*TMath::RadToDeg(),phimax*TMath::RadToDeg()) ;
667 fhPrimEtaPhi->SetYTitle("#phi (deg)");
668 fhPrimEtaPhi->SetXTitle("p_{T} (GeV/c)");
669 outputContainer->Add(fhPrimEtaPhi) ;
671 fhPrimEtaAccPhi = new TH2D("hPrimEtaAccPhi","Azimuthal of primary eta with accepted daughters",nptbins,ptmin,ptmax, nphibins,phimin*TMath::RadToDeg(),phimax*TMath::RadToDeg()) ;
672 fhPrimEtaAccPhi->SetYTitle("#phi (deg)");
673 fhPrimEtaAccPhi->SetXTitle("p_{T} (GeV/c)");
674 outputContainer->Add(fhPrimEtaAccPhi) ;
677 fhPrimPi0OpeningAngle = new TH2D
678 ("hPrimPi0OpeningAngle","Angle between all primary #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,100,0,0.5);
679 fhPrimPi0OpeningAngle->SetYTitle("#theta(rad)");
680 fhPrimPi0OpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
681 outputContainer->Add(fhPrimPi0OpeningAngle) ;
683 fhPrimPi0CosOpeningAngle = new TH2D
684 ("hPrimPi0CosOpeningAngle","Cosinus of angle between all primary #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,100,-1,1);
685 fhPrimPi0CosOpeningAngle->SetYTitle("cos (#theta) ");
686 fhPrimPi0CosOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
687 outputContainer->Add(fhPrimPi0CosOpeningAngle) ;
689 for(Int_t i = 0; i<13; i++){
690 fhMCOrgMass[i] = new TH2D(Form("hMCOrgMass_%d",i),Form("mass vs pt, origin %d",i),nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
691 fhMCOrgMass[i]->SetXTitle("p_{T} (GeV/c)");
692 fhMCOrgMass[i]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
693 outputContainer->Add(fhMCOrgMass[i]) ;
695 fhMCOrgAsym[i]= new TH2D(Form("hMCOrgAsym_%d",i),Form("asymmetry vs pt, origin %d",i),nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax) ;
696 fhMCOrgAsym[i]->SetXTitle("p_{T} (GeV/c)");
697 fhMCOrgAsym[i]->SetYTitle("A");
698 outputContainer->Add(fhMCOrgAsym[i]) ;
700 fhMCOrgDeltaEta[i] = new TH2D(Form("hMCOrgDeltaEta_%d",i),Form("#Delta #eta of pair vs pt, origin %d",i),nptbins,ptmin,ptmax,netabins,-1.4,1.4) ;
701 fhMCOrgDeltaEta[i]->SetXTitle("p_{T} (GeV/c)");
702 fhMCOrgDeltaEta[i]->SetYTitle("#Delta #eta");
703 outputContainer->Add(fhMCOrgDeltaEta[i]) ;
705 fhMCOrgDeltaPhi[i]= new TH2D(Form("hMCOrgDeltaPhi_%d",i),Form("#Delta #phi of pair vs p_{T}, origin %d",i),nptbins,ptmin,ptmax,nphibins,-0.7,0.7) ;
706 fhMCOrgDeltaPhi[i]->SetXTitle("p_{T} (GeV/c)");
707 fhMCOrgDeltaPhi[i]->SetYTitle("#Delta #phi (rad)");
708 outputContainer->Add(fhMCOrgDeltaPhi[i]) ;
713 fhMCPi0MassPtTrue = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
714 fhMCPi0MassPtRec = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
715 fhMCPi0PtTruePtRec = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
716 fhMCEtaMassPtRec = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
717 fhMCEtaMassPtTrue = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
718 fhMCEtaPtTruePtRec = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
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 Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
724 fhMCPi0MassPtRec[index] = new TH2D(Form("hMCPi0MassPtRec_pt%d_cell%d_asym%d",ipt,icell,iasym),
725 Form("Reconstructed Mass vs reconstructed p_T of true #pi^{0} cluster pairs for pt >%2.2f, ncell>%d and asym >%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]),
726 nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
727 fhMCPi0MassPtRec[index]->SetXTitle("p_{T, reconstructed} (GeV/c)");
728 fhMCPi0MassPtRec[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
729 outputContainer->Add(fhMCPi0MassPtRec[index]) ;
731 fhMCPi0MassPtTrue[index] = new TH2D(Form("hMCPi0MassPtTrue_pt%d_cell%d_asym%d",ipt,icell,iasym),
732 Form("Reconstructed Mass vs generated p_T of true #pi^{0} cluster pairs for pt >%2.2f, ncell>%d and asym >%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]),
733 nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
734 fhMCPi0MassPtTrue[index]->SetXTitle("p_{T, generated} (GeV/c)");
735 fhMCPi0MassPtTrue[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
736 outputContainer->Add(fhMCPi0MassPtTrue[index]) ;
738 fhMCPi0PtTruePtRec[index] = new TH2D(Form("hMCPi0PtTruePtRec_pt%d_cell%d_asym%d",ipt,icell,iasym),
739 Form("Generated vs reconstructed p_T of true #pi^{0} cluster pairs, 0.01 < rec. mass < 0.17 MeV/c^{2} for pt >%2.2f, ncell>%d and asym >%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]),
740 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax) ;
741 fhMCPi0PtTruePtRec[index]->SetXTitle("p_{T, generated} (GeV/c)");
742 fhMCPi0PtTruePtRec[index]->SetYTitle("p_{T, reconstructed} (GeV/c)");
743 outputContainer->Add(fhMCPi0PtTruePtRec[index]) ;
745 fhMCEtaMassPtRec[index] = new TH2D(Form("hMCEtaMassPtRec_pt%d_cell%d_asym%d",ipt,icell,iasym),
746 Form("Reconstructed Mass vs reconstructed p_T of true #eta cluster pairs for pt >%2.2f, ncell>%d and asym >%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]),
747 nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
748 fhMCEtaMassPtRec[index]->SetXTitle("p_{T, generated} (GeV/c)");
749 fhMCEtaMassPtRec[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
750 outputContainer->Add(fhMCEtaMassPtRec[index]) ;
752 fhMCEtaMassPtTrue[index] = new TH2D(Form("hMCEtaMassPtTrue_pt%d_cell%d_asym%d",ipt,icell,iasym),
753 Form("Reconstructed Mass vs generated p_T of true #eta cluster pairs for pt >%2.2f, ncell>%d and asym >%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]),
754 nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
755 fhMCEtaMassPtTrue[index]->SetXTitle("p_{T, generated} (GeV/c)");
756 fhMCEtaMassPtTrue[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
757 outputContainer->Add(fhMCEtaMassPtTrue[index]) ;
759 fhMCEtaPtTruePtRec[index] = new TH2D(Form("hMCEtaPtTruePtRec_pt%d_cell%d_asym%d",ipt,icell,iasym),
760 Form("Generated vs reconstructed p_T of true #eta cluster pairs, 0.01 < rec. mass < 0.17 MeV/c^{2} for pt >%2.2f, ncell>%d and asym >%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]),
761 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax) ;
762 fhMCEtaPtTruePtRec[index]->SetXTitle("p_{T, generated} (GeV/c)");
763 fhMCEtaPtTruePtRec[index]->SetYTitle("p_{T, reconstructed} (GeV/c)");
764 outputContainer->Add(fhMCEtaPtTruePtRec[index]) ;
770 fhMCPi0MassPtTrue = new TH2D*[1];
771 fhMCPi0PtTruePtRec = new TH2D*[1];
772 fhMCEtaMassPtTrue = new TH2D*[1];
773 fhMCEtaPtTruePtRec = new TH2D*[1];
775 fhMCPi0MassPtTrue[0] = new TH2D("hMCPi0MassPtTrue","Reconstructed Mass vs generated p_T of true #pi^{0} cluster pairs",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
776 fhMCPi0MassPtTrue[0]->SetXTitle("p_{T, generated} (GeV/c)");
777 fhMCPi0MassPtTrue[0]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
778 outputContainer->Add(fhMCPi0MassPtTrue[0]) ;
780 fhMCPi0PtTruePtRec[0]= new TH2D("hMCPi0PtTruePtRec","Generated vs reconstructed p_T of true #pi^{0} cluster pairs, 0.01 < rec. mass < 0.17 MeV/c^{2}",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax) ;
781 fhMCPi0PtTruePtRec[0]->SetXTitle("p_{T, generated} (GeV/c)");
782 fhMCPi0PtTruePtRec[0]->SetYTitle("p_{T, reconstructed} (GeV/c)");
783 outputContainer->Add(fhMCPi0PtTruePtRec[0]) ;
785 fhMCEtaMassPtTrue[0] = new TH2D("hMCEtaMassPtTrue","Reconstructed Mass vs generated p_T of true #eta cluster pairs",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
786 fhMCEtaMassPtTrue[0]->SetXTitle("p_{T, generated} (GeV/c)");
787 fhMCEtaMassPtTrue[0]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
788 outputContainer->Add(fhMCEtaMassPtTrue[0]) ;
790 fhMCEtaPtTruePtRec[0]= new TH2D("hMCEtaPtTruePtRec","Generated vs reconstructed p_T of true #eta cluster pairs, 0.01 < rec. mass < 0.17 MeV/c^{2}",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax) ;
791 fhMCEtaPtTruePtRec[0]->SetXTitle("p_{T, generated} (GeV/c)");
792 fhMCEtaPtTruePtRec[0]->SetYTitle("p_{T, reconstructed} (GeV/c)");
793 outputContainer->Add(fhMCEtaPtTruePtRec[0]) ;
797 TString * pairname = new TString[fNModules+3];
798 if(fCalorimeter=="EMCAL"){
799 pairname[0]="A side (0-2)";
800 pairname[1]="C side (1-3)";
801 pairname[2]="Sector 0 (0-1)";
802 pairname[3]="Sector 1 (2-3)";
803 pairname[4]="Cluster in different SM";
804 pairname[5]="SM 0 and SM3";
805 pairname[6]="SM 1 and SM2";
806 for(Int_t i = 7 ; i < fNModules ; i++) pairname[i]="";}
807 if(fCalorimeter=="PHOS") {
811 for(Int_t i = 3 ; i < fNModules ; i++) pairname[i]="";}
813 for(Int_t imod=0; imod<fNModules; imod++){
814 //Module dependent invariant mass
815 snprintf(key, buffersize,"hReMod_%d",imod) ;
816 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for Module %d",imod) ;
817 fhReMod[imod] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
818 fhReMod[imod]->SetXTitle("p_{T} (GeV/c)");
819 fhReMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
820 outputContainer->Add(fhReMod[imod]) ;
822 snprintf(key, buffersize,"hReDiffMod_%d",imod) ;
823 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for Different Modules: %s",(pairname[imod]).Data()) ;
824 fhReDiffMod[imod] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
825 fhReDiffMod[imod]->SetXTitle("p_{T} (GeV/c)");
826 fhReDiffMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
827 outputContainer->Add(fhReDiffMod[imod]) ;
830 snprintf(key, buffersize,"hMiMod_%d",imod) ;
831 snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for Module %d",imod) ;
832 fhMiMod[imod] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
833 fhMiMod[imod]->SetXTitle("p_{T} (GeV/c)");
834 fhMiMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
835 outputContainer->Add(fhMiMod[imod]) ;
837 snprintf(key, buffersize,"hMiDiffMod_%d",imod) ;
838 snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for Different Modules: %s",(pairname[imod]).Data()) ;
839 fhMiDiffMod[imod] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
840 fhMiDiffMod[imod]->SetXTitle("p_{T} (GeV/c)");
841 fhMiDiffMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
842 outputContainer->Add(fhMiDiffMod[imod]) ;
847 for (Int_t imod=4; imod<7; imod++) {
849 snprintf(key, buffersize,"hReDiffMod_%d",imod) ;
850 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for Different Modules: %s",(pairname[imod]).Data()) ;
851 fhReDiffMod[imod] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
852 fhReDiffMod[imod]->SetXTitle("p_{T} (GeV/c)");
853 fhReDiffMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
854 outputContainer->Add(fhReDiffMod[imod]) ;
857 snprintf(key, buffersize,"hMiDiffMod_%d",imod) ;
858 snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for Different Modules: %s",(pairname[imod]).Data()) ;
859 fhMiDiffMod[imod] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
860 fhMiDiffMod[imod]->SetXTitle("p_{T} (GeV/c)");
861 fhMiDiffMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
862 outputContainer->Add(fhMiDiffMod[imod]) ;
868 // for(Int_t i = 0; i < outputContainer->GetEntries() ; i++){
870 // printf("Histogram %d, name: %s\n ",i, outputContainer->At(i)->GetName());
874 return outputContainer;
877 //_________________________________________________________________________________________________________________________________________________
878 void AliAnaPi0::Print(const Option_t * /*opt*/) const
880 //Print some relevant parameters set for the analysis
881 printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
882 AliAnaPartCorrBaseClass::Print(" ");
884 printf("Number of bins in Centrality: %d \n",fNCentrBin) ;
885 printf("Number of bins in Z vert. pos: %d \n",GetNZvertBin()) ;
886 printf("Number of bins in Reac. Plain: %d \n",GetNRPBin()) ;
887 printf("Depth of event buffer: %d \n",fNmaxMixEv) ;
888 printf("Pair in same Module: %d \n",fSameSM) ;
890 printf("Z vertex position: -%2.3f < z < %2.3f \n",GetZvertexCut(),GetZvertexCut()) ;
891 printf("Number of modules: %d \n",fNModules) ;
892 printf("Select pairs with their angle: %d, edep %d, min angle %2.3f, max angle %2.3f \n",fUseAngleCut, fUseAngleEDepCut, fAngleCut, fAngleMaxCut) ;
893 printf("Asymmetry cuts: n = %d, \n",fNAsymCuts) ;
894 printf("\tasymmetry < ");
895 for(Int_t i = 0; i < fNAsymCuts; i++) printf("%2.2f ",fAsymCuts[i]);
898 printf("PID selection bits: n = %d, \n",fNPIDBits) ;
899 printf("\tPID bit = ");
900 for(Int_t i = 0; i < fNPIDBits; i++) printf("%d ",fPIDBits[i]);
904 printf("pT cuts: n = %d, \n",fNPtCuts) ;
906 for(Int_t i = 0; i < fNPtCuts; i++) printf("%2.2f ",fPtCuts[i]);
909 printf("N cell in cluster cuts: n = %d, \n",fNCellNCuts) ;
910 printf("\tnCell > ");
911 for(Int_t i = 0; i < fNCellNCuts; i++) printf("%d ",fCellNCuts[i]);
915 printf("------------------------------------------------------\n") ;
918 //_____________________________________________________________
919 void AliAnaPi0::FillAcceptanceHistograms(){
920 //Fill acceptance histograms if MC data is available
922 if(GetReader()->ReadStack()){
923 AliStack * stack = GetMCStack();
925 for(Int_t i=0 ; i<stack->GetNprimary(); i++){
926 TParticle * prim = stack->Particle(i) ;
927 Int_t pdg = prim->GetPdgCode();
928 if( pdg == 111 || pdg == 221){
929 Double_t pi0Pt = prim->Pt() ;
930 //printf("pi0, pt %2.2f\n",pi0Pt);
931 if(prim->Energy() == TMath::Abs(prim->Pz())) continue ; //Protection against floating point exception
932 Double_t pi0Y = 0.5*TMath::Log((prim->Energy()-prim->Pz())/(prim->Energy()+prim->Pz())) ;
933 Double_t phi = TMath::RadToDeg()*prim->Phi() ;
935 if(TMath::Abs(pi0Y) < 0.5){
936 fhPrimPi0Pt->Fill(pi0Pt) ;
938 fhPrimPi0Y ->Fill(pi0Pt, pi0Y) ;
939 fhPrimPi0Phi->Fill(pi0Pt, phi) ;
942 if(TMath::Abs(pi0Y) < 0.5){
943 fhPrimEtaPt->Fill(pi0Pt) ;
945 fhPrimEtaY ->Fill(pi0Pt, pi0Y) ;
946 fhPrimEtaPhi->Fill(pi0Pt, phi) ;
948 //Check if both photons hit Calorimeter
949 if(prim->GetNDaughters()!=2) return; //Only interested in 2 gamma decay
950 Int_t iphot1=prim->GetFirstDaughter() ;
951 Int_t iphot2=prim->GetLastDaughter() ;
952 if(iphot1>-1 && iphot1<stack->GetNtrack() && iphot2>-1 && iphot2<stack->GetNtrack()){
953 TParticle * phot1 = stack->Particle(iphot1) ;
954 TParticle * phot2 = stack->Particle(iphot2) ;
955 if(phot1 && phot2 && phot1->GetPdgCode()==22 && phot2->GetPdgCode()==22){
956 //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",
957 // phot1->Pt(), phot1->Phi()*180./3.1415, phot1->Eta(), phot2->Pt(), phot2->Phi()*180./3.1415, phot2->Eta());
959 TLorentzVector lv1, lv2;
960 phot1->Momentum(lv1);
961 phot2->Momentum(lv2);
963 Bool_t inacceptance = kFALSE;
964 if(fCalorimeter == "PHOS"){
965 if(GetPHOSGeometry() && GetCaloUtils()->IsPHOSGeoMatrixSet()){
968 if(GetPHOSGeometry()->ImpactOnEmc(phot1,mod,z,x) && GetPHOSGeometry()->ImpactOnEmc(phot2,mod,z,x))
969 inacceptance = kTRUE;
970 if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
974 if(GetFiducialCut()->IsInFiducialCut(lv1,fCalorimeter) && GetFiducialCut()->IsInFiducialCut(lv2,fCalorimeter))
975 inacceptance = kTRUE ;
976 if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
980 else if(fCalorimeter == "EMCAL" && GetCaloUtils()->IsEMCALGeoMatrixSet()){
981 if(GetEMCALGeometry()){
986 GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(phot1->Eta(),phot1->Phi(),absID1);
987 GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(phot2->Eta(),phot2->Phi(),absID2);
989 if( absID1 >= 0 && absID2 >= 0)
990 inacceptance = kTRUE;
992 // if(GetEMCALGeometry()->Impact(phot1) && GetEMCALGeometry()->Impact(phot2))
993 // inacceptance = kTRUE;
994 if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
997 if(GetFiducialCut()->IsInFiducialCut(lv1,fCalorimeter) && GetFiducialCut()->IsInFiducialCut(lv2,fCalorimeter))
998 inacceptance = kTRUE ;
999 if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
1005 fhPrimPi0AccPt ->Fill(pi0Pt) ;
1006 fhPrimPi0AccPhi->Fill(pi0Pt, phi) ;
1007 fhPrimPi0AccY ->Fill(pi0Pt, pi0Y) ;
1008 Double_t angle = lv1.Angle(lv2.Vect());
1009 fhPrimPi0OpeningAngle ->Fill(pi0Pt,angle);
1010 fhPrimPi0CosOpeningAngle->Fill(pi0Pt,TMath::Cos(angle));
1013 fhPrimEtaAccPt ->Fill(pi0Pt) ;
1014 fhPrimEtaAccPhi->Fill(pi0Pt, phi) ;
1015 fhPrimEtaAccY ->Fill(pi0Pt, pi0Y) ;
1019 }//Check daughters exist
1020 }// Primary pi0 or eta
1021 }//loop on primaries
1022 }//stack exists and data is MC
1024 else if(GetReader()->ReadAODMCParticles()){
1026 TClonesArray * mcparticles = GetReader()->GetAODMCParticles(0);
1028 Int_t nprim = mcparticles->GetEntriesFast();
1032 AliAODMCParticle * prim = (AliAODMCParticle *) mcparticles->At(i);
1033 Int_t pdg = prim->GetPdgCode();
1034 if( pdg == 111 || pdg == 221){
1035 Double_t pi0Pt = prim->Pt() ;
1036 //printf("pi0, pt %2.2f\n",pi0Pt);
1037 if(prim->E() == TMath::Abs(prim->Pz())) continue ; //Protection against floating point exception
1038 Double_t pi0Y = 0.5*TMath::Log((prim->E()-prim->Pz())/(prim->E()+prim->Pz())) ;
1039 Double_t phi = TMath::RadToDeg()*prim->Phi() ;
1041 if(TMath::Abs(pi0Y) < 0.5){
1042 fhPrimPi0Pt->Fill(pi0Pt) ;
1044 fhPrimPi0Y ->Fill(pi0Pt, pi0Y) ;
1045 fhPrimPi0Phi->Fill(pi0Pt, phi) ;
1047 else if(pdg == 221){
1048 if(TMath::Abs(pi0Y) < 0.5){
1049 fhPrimEtaPt->Fill(pi0Pt) ;
1051 fhPrimEtaY ->Fill(pi0Pt, pi0Y) ;
1052 fhPrimEtaPhi->Fill(pi0Pt, phi) ;
1054 //Check if both photons hit Calorimeter
1055 if(prim->GetNDaughters()!=2) return; //Only interested in 2 gamma decay
1056 Int_t iphot1=prim->GetDaughter(0) ;
1057 Int_t iphot2=prim->GetDaughter(1) ;
1058 if(iphot1>-1 && iphot1<nprim && iphot2>-1 && iphot2<nprim){
1059 AliAODMCParticle * phot1 = (AliAODMCParticle *) mcparticles->At(iphot1);
1060 AliAODMCParticle * phot2 = (AliAODMCParticle *) mcparticles->At(iphot2);
1061 if(phot1 && phot2 && phot1->GetPdgCode()==22 && phot2->GetPdgCode()==22){
1062 TLorentzVector lv1, lv2;
1063 lv1.SetPxPyPzE(phot1->Px(),phot1->Py(),phot1->Pz(),phot1->E());
1064 lv2.SetPxPyPzE(phot2->Px(),phot2->Py(),phot2->Pz(),phot2->E());
1066 Bool_t inacceptance = kFALSE;
1067 if(fCalorimeter == "PHOS"){
1068 if(GetPHOSGeometry() && GetCaloUtils()->IsPHOSGeoMatrixSet()){
1071 Double_t vtx []={phot1->Xv(),phot1->Yv(),phot1->Zv()};
1072 Double_t vtx2[]={phot2->Xv(),phot2->Yv(),phot2->Zv()};
1073 if(GetPHOSGeometry()->ImpactOnEmc(vtx, phot1->Theta(),phot1->Phi(),mod,z,x) &&
1074 GetPHOSGeometry()->ImpactOnEmc(vtx2,phot2->Theta(),phot2->Phi(),mod,z,x))
1075 inacceptance = kTRUE;
1076 if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
1080 if(GetFiducialCut()->IsInFiducialCut(lv1,fCalorimeter) && GetFiducialCut()->IsInFiducialCut(lv2,fCalorimeter))
1081 inacceptance = kTRUE ;
1082 if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
1086 else if(fCalorimeter == "EMCAL" && GetCaloUtils()->IsEMCALGeoMatrixSet()){
1087 if(GetEMCALGeometry()){
1092 //TVector3 vtx(phot1->Xv(),phot1->Yv(),phot1->Zv());
1093 //TVector3 vimpact(0,0,0);
1095 //GetEMCALGeometry()->ImpactOnEmcal(vtx,phot1->Theta(),phot1->Phi(),absID1,vimpact);
1096 //TVector3 vtx2(phot2->Xv(),phot2->Yv(),phot2->Zv());
1097 //TVector3 vimpact2(0,0,0);
1098 //GetEMCALGeometry()->ImpactOnEmcal(vtx2,phot2->Theta(),phot2->Phi(),absID2,vimpact2);
1100 GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(phot1->Eta(),phot1->Phi(),absID1);
1101 GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(phot2->Eta(),phot2->Phi(),absID2);
1103 // if(TMath::Abs(phot1->Eta()) < 0.7 && phot1->Phi() > 80*TMath::DegToRad() && phot1->Phi() < 120*TMath::DegToRad() )
1104 // printf("Phot1 ccepted? %d\n",absID1);
1105 // if(TMath::Abs(phot2->Eta()) < 0.7 && phot2->Phi() > 80*TMath::DegToRad() && phot2->Phi() < 120*TMath::DegToRad() )
1106 // printf("Phot2 accepted? %d\n",absID2);
1108 if( absID1 >= 0 && absID2 >= 0)
1109 inacceptance = kTRUE;
1111 // if(pdg==111 && inacceptance) printf("2 photons: photon 1: absId %d, pt %2.2f, phi %3.2f, eta %1.2f; photon 2: absId %d, pt %2.2f, phi %3.2f, eta %1.2f\n",
1112 // absID1,phot1->Pt(), phot1->Phi()*TMath::RadToDeg(), phot1->Eta(),
1113 // absID2,phot2->Pt(), phot2->Phi()*TMath::RadToDeg(), phot2->Eta());
1117 if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
1120 if(GetFiducialCut()->IsInFiducialCut(lv1,fCalorimeter) && GetFiducialCut()->IsInFiducialCut(lv2,fCalorimeter))
1121 inacceptance = kTRUE ;
1122 if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
1128 // printf("ACCEPTED pi0: pt %2.2f, phi %3.2f, eta %1.2f\n",pi0Pt,phi,pi0Y);
1129 fhPrimPi0AccPt ->Fill(pi0Pt) ;
1130 fhPrimPi0AccPhi->Fill(pi0Pt, phi) ;
1131 fhPrimPi0AccY ->Fill(pi0Pt, pi0Y) ;
1132 Double_t angle = lv1.Angle(lv2.Vect());
1133 fhPrimPi0OpeningAngle ->Fill(pi0Pt,angle);
1134 fhPrimPi0CosOpeningAngle->Fill(pi0Pt,TMath::Cos(angle));
1137 fhPrimEtaAccPt ->Fill(pi0Pt) ;
1138 fhPrimEtaAccPhi->Fill(pi0Pt, phi) ;
1139 fhPrimEtaAccY ->Fill(pi0Pt, pi0Y) ;
1143 }//Check daughters exist
1144 }// Primary pi0 or eta
1145 }//loop on primaries
1146 }//stack exists and data is MC
1152 //_____________________________________________________________
1153 void AliAnaPi0::FillMCVersusRecDataHistograms(const Int_t index1, const Int_t index2,
1154 const Float_t pt1, const Float_t pt2,
1155 const Int_t ncell1, const Int_t ncell2,
1156 const Double_t mass, const Double_t pt, const Double_t asym,
1157 const Double_t deta, const Double_t dphi){
1158 //Do some MC checks on the origin of the pair, is there any common ancestor and if there is one, who?
1159 //Adjusted for Pythia, need to see what to do for other generators.
1160 //Array of histograms ordered as follows: 0-Photon, 1-electron, 2-pi0, 3-eta, 4-a-proton, 5-a-neutron, 6-stable particles,
1161 // 7-other decays, 8-string, 9-final parton, 10-initial parton, intermediate, 11-colliding proton, 12-unrelated
1164 Int_t ancStatus = 0;
1165 TLorentzVector ancMomentum;
1166 Int_t ancLabel = GetMCAnalysisUtils()->CheckCommonAncestor(index1, index2,
1167 GetReader(), ancPDG, ancStatus,ancMomentum);
1169 if(GetDebug() > 1) printf("AliAnaPi0::FillMCVersusRecDataHistograms() - Common ancestor label %d, pdg %d, name %s, status %d; \n",
1170 ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus);
1173 if(ancPDG==22){//gamma
1174 fhMCOrgMass[0]->Fill(pt,mass);
1175 fhMCOrgAsym[0]->Fill(pt,asym);
1176 fhMCOrgDeltaEta[0]->Fill(pt,deta);
1177 fhMCOrgDeltaPhi[0]->Fill(pt,dphi);
1179 else if(TMath::Abs(ancPDG)==11){//e
1180 fhMCOrgMass[1]->Fill(pt,mass);
1181 fhMCOrgAsym[1]->Fill(pt,asym);
1182 fhMCOrgDeltaEta[1]->Fill(pt,deta);
1183 fhMCOrgDeltaPhi[1]->Fill(pt,dphi);
1185 else if(ancPDG==111){//Pi0
1186 fhMCOrgMass[2]->Fill(pt,mass);
1187 fhMCOrgAsym[2]->Fill(pt,asym);
1188 fhMCOrgDeltaEta[2]->Fill(pt,deta);
1189 fhMCOrgDeltaPhi[2]->Fill(pt,dphi);
1190 if(fMultiCutAnaSim){
1191 for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
1192 for(Int_t icell=0; icell<fNCellNCuts; icell++){
1193 for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
1194 Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
1195 if(pt1 > fPtCuts[ipt] && pt2 > fPtCuts[ipt] &&
1196 asym < fAsymCuts[iasym] &&
1197 ncell1 >= fCellNCuts[icell] && ncell2 >= fCellNCuts[icell]){
1198 fhMCPi0MassPtRec [index]->Fill(pt,mass);
1199 fhMCPi0MassPtTrue[index]->Fill(ancMomentum.Pt(),mass);
1200 if(mass < 0.17 && mass > 0.1) fhMCPi0PtTruePtRec[index]->Fill(ancMomentum.Pt(),pt);
1201 }//pass the different cuts
1202 }// pid bit cut loop
1205 }//Multi cut ana sim
1207 fhMCPi0MassPtTrue[0]->Fill(ancMomentum.Pt(),mass);
1208 if(mass < 0.17 && mass > 0.1) fhMCPi0PtTruePtRec[0]->Fill(ancMomentum.Pt(),pt);
1211 else if(ancPDG==221){//Eta
1212 fhMCOrgMass[3]->Fill(pt,mass);
1213 fhMCOrgAsym[3]->Fill(pt,asym);
1214 fhMCOrgDeltaEta[3]->Fill(pt,deta);
1215 fhMCOrgDeltaPhi[3]->Fill(pt,dphi);
1216 if(fMultiCutAnaSim){
1217 for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
1218 for(Int_t icell=0; icell<fNCellNCuts; icell++){
1219 for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
1220 Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
1221 if(pt1 > fPtCuts[ipt] && pt2 > fPtCuts[ipt] &&
1222 asym < fAsymCuts[iasym] &&
1223 ncell1 >= fCellNCuts[icell] && ncell2 >= fCellNCuts[icell]){
1224 fhMCEtaMassPtRec [index]->Fill(pt,mass);
1225 fhMCEtaMassPtTrue[index]->Fill(ancMomentum.Pt(),mass);
1226 if(mass < 0.65 && mass > 0.45) fhMCEtaPtTruePtRec[index]->Fill(ancMomentum.Pt(),pt);
1227 }//pass the different cuts
1228 }// pid bit cut loop
1231 } //Multi cut ana sim
1233 fhMCEtaMassPtTrue[0]->Fill(ancMomentum.Pt(),mass);
1234 if(mass < 0.65 && mass > 0.45) fhMCEtaPtTruePtRec[0]->Fill(ancMomentum.Pt(),pt);
1237 else if(ancPDG==-2212){//AProton
1238 fhMCOrgMass[4]->Fill(pt,mass);
1239 fhMCOrgAsym[4]->Fill(pt,asym);
1240 fhMCOrgDeltaEta[4]->Fill(pt,deta);
1241 fhMCOrgDeltaPhi[4]->Fill(pt,dphi);
1243 else if(ancPDG==-2112){//ANeutron
1244 fhMCOrgMass[5]->Fill(pt,mass);
1245 fhMCOrgAsym[5]->Fill(pt,asym);
1246 fhMCOrgDeltaEta[5]->Fill(pt,deta);
1247 fhMCOrgDeltaPhi[5]->Fill(pt,dphi);
1249 else if(TMath::Abs(ancPDG)==13){//muons
1250 fhMCOrgMass[6]->Fill(pt,mass);
1251 fhMCOrgAsym[6]->Fill(pt,asym);
1252 fhMCOrgDeltaEta[6]->Fill(pt,deta);
1253 fhMCOrgDeltaPhi[6]->Fill(pt,dphi);
1255 else if (TMath::Abs(ancPDG) > 100 && ancLabel > 7) {
1256 if(ancStatus==1){//Stable particles, converted? not decayed resonances
1257 fhMCOrgMass[6]->Fill(pt,mass);
1258 fhMCOrgAsym[6]->Fill(pt,asym);
1259 fhMCOrgDeltaEta[6]->Fill(pt,deta);
1260 fhMCOrgDeltaPhi[6]->Fill(pt,dphi);
1262 else{//resonances and other decays, more hadron conversions?
1263 fhMCOrgMass[7]->Fill(pt,mass);
1264 fhMCOrgAsym[7]->Fill(pt,asym);
1265 fhMCOrgDeltaEta[7]->Fill(pt,deta);
1266 fhMCOrgDeltaPhi[7]->Fill(pt,dphi);
1269 else {//Partons, colliding protons, strings, intermediate corrections
1270 if(ancStatus==11 || ancStatus==12){//String fragmentation
1271 fhMCOrgMass[8]->Fill(pt,mass);
1272 fhMCOrgAsym[8]->Fill(pt,asym);
1273 fhMCOrgDeltaEta[8]->Fill(pt,deta);
1274 fhMCOrgDeltaPhi[8]->Fill(pt,dphi);
1276 else if (ancStatus==21){
1277 if(ancLabel < 2) {//Colliding protons
1278 fhMCOrgMass[11]->Fill(pt,mass);
1279 fhMCOrgAsym[11]->Fill(pt,asym);
1280 fhMCOrgDeltaEta[11]->Fill(pt,deta);
1281 fhMCOrgDeltaPhi[11]->Fill(pt,dphi);
1282 }//colliding protons
1283 else if(ancLabel < 6){//partonic initial states interactions
1284 fhMCOrgMass[9]->Fill(pt,mass);
1285 fhMCOrgAsym[9]->Fill(pt,asym);
1286 fhMCOrgDeltaEta[9]->Fill(pt,deta);
1287 fhMCOrgDeltaPhi[9]->Fill(pt,dphi);
1289 else if(ancLabel < 8){//Final state partons radiations?
1290 fhMCOrgMass[10]->Fill(pt,mass);
1291 fhMCOrgAsym[10]->Fill(pt,asym);
1292 fhMCOrgDeltaEta[10]->Fill(pt,deta);
1293 fhMCOrgDeltaPhi[10]->Fill(pt,dphi);
1296 printf("AliAnaPi0::FillMCVersusRecDataHistograms() - Check ** Common ancestor label %d, pdg %d, name %s, status %d; \n",
1297 ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus);
1301 printf("AliAnaPi0::FillMCVersusRecDataHistograms() - Check *** Common ancestor label %d, pdg %d, name %s, status %d; \n",
1302 ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus);
1304 }////Partons, colliding protons, strings, intermediate corrections
1306 else { //ancLabel <= -1
1307 //printf("Not related at all label = %d\n",ancLabel);
1308 fhMCOrgMass[12]->Fill(pt,mass);
1309 fhMCOrgAsym[12]->Fill(pt,asym);
1310 fhMCOrgDeltaEta[12]->Fill(pt,deta);
1311 fhMCOrgDeltaPhi[12]->Fill(pt,dphi);
1315 //____________________________________________________________________________________________________________________________________________________
1316 void AliAnaPi0::MakeAnalysisFillHistograms()
1318 //Process one event and extract photons from AOD branch
1319 // filled with AliAnaPhoton and fill histos with invariant mass
1321 //In case of simulated data, fill acceptance histograms
1322 if(IsDataMC())FillAcceptanceHistograms();
1323 if (GetReader()->GetEventNumber()%10000 == 0)
1324 printf("--- Event %d ---\n",GetReader()->GetEventNumber());
1325 //Init some variables
1326 //Int_t iRun = (GetReader()->GetInputEvent())->GetRunNumber() ;
1327 Int_t nPhot = GetInputAODBranch()->GetEntriesFast() ;
1330 Float_t eClusTot = 0;
1331 Float_t eCellTot = 0;
1332 Float_t eDenClus = 0;
1333 Float_t eDenCell = 0;
1335 // Float_t rtmp = 0;
1336 // Float_t rtmpw = 0;
1338 // Float_t rxzw = 0;
1341 // Float_t emax = 0;
1344 printf("AliAnaPi0::MakeAnalysisFillHistograms() - Photon entries %d\n", nPhot);
1346 //If less than photon 2 entries in the list, skip this event
1347 if(nPhot < 2 ) return ;
1349 // Count the number of clusters and cells, in case multiplicity bins dependent on such numbers
1351 if(fCalorimeter=="EMCAL"){
1352 nClus = GetAODEMCAL() ->GetEntriesFast();
1353 nCell = GetEMCALCells()->GetNumberOfCells();
1354 for(Int_t icl=0; icl < nClus; icl++) {
1355 Float_t e1 = ((AliVCluster*)GetAODEMCAL()->At(icl))->E();
1357 // if(e1 > emax) emax = e1;
1358 // ((AliVCluster*)GetAODEMCAL()->At(icl))->GetPosition(pos1);
1359 // for(Int_t icl2=icl+1; icl2 < nClus; icl2++) {
1360 // Float_t e2 = ((AliVCluster*)GetAODEMCAL()->At(icl2))->E();
1361 // ((AliVCluster*)GetAODEMCAL()->At(icl2))->GetPosition(pos2);
1362 // rtmp = TMath::Sqrt((pos1[0]-pos2[0])*(pos1[0]-pos2[0]) + (pos1[2]-pos2[2])*(pos1[2]-pos2[2]));
1363 // rtmpw = TMath::Sqrt((pos1[0]*e1-pos2[0]*e2)*(pos1[0]*e1-pos2[0]*e2) + (pos1[2]*e1-pos2[2]*e2)*(pos1[2]*e1-pos2[2]*e2))/(e1+e2);
1367 // fhClusterPairDist ->Fill(rtmp);
1368 // fhClusterPairDistWeight->Fill(rtmpw);
1369 // //printf("Distance: %f; weighted %f\n ",rtmp,rtmp/(e1+((AliVCluster*)GetAODEMCAL()->At(icl2))->E()));
1371 // }// second cluster loop
1374 for(Int_t jce=0; jce < nCell; jce++) eCellTot += GetEMCALCells()->GetAmplitude(jce);
1377 nClus = GetAODPHOS() ->GetEntriesFast();
1378 nCell = GetPHOSCells()->GetNumberOfCells();
1379 for(Int_t icl=0; icl < nClus; icl++) {
1380 Float_t e1 = ((AliVCluster*)GetAODPHOS()->At(icl))->E();
1382 // ((AliVCluster*)GetAODPHOS()->At(icl))->GetPosition(pos1);
1383 // for(Int_t icl2=icl+1; icl2 < nClus; icl2++) {
1384 // Float_t e2 = ((AliVCluster*)GetAODPHOS()->At(icl2))->E();
1385 // ((AliVCluster*)GetAODPHOS()->At(icl2))->GetPosition(pos2);
1386 // rtmp = TMath::Sqrt((pos1[0]-pos2[0])*(pos1[0]-pos2[0]) + (pos1[2]-pos2[2])*(pos1[2]-pos2[2]));
1387 // rtmpw = TMath::Sqrt((pos1[0]*e1-pos2[0]*e2)*(pos1[0]*e1-pos2[0]*e2) + (pos1[2]*e1-pos2[2]*e2)*(pos1[2]*e1-pos2[2]*e2))/(e1+e2);
1391 // fhClusterPairDist ->Fill(rtmp);
1392 // fhClusterPairDistWeight->Fill(rtmpw);
1393 // }// second cluster loop
1395 for(Int_t jce=0; jce < nCell; jce++) eCellTot += GetPHOSCells()->GetAmplitude(jce);
1398 printf("AliAnaPi0::MakeAnalysisFillHistograms() - # Clusters %d, sum cluster E per SM %f,# Cells %d, sum cell E per SM %f\n", nClus,eClusTot,nCell,eCellTot);
1400 //Fill histograms with "energy density", ncell and nclust will be > 0 since there are at least 2 "photons"
1401 eDenClus = eClusTot/nClus;
1402 eDenCell = eCellTot/nCell;
1403 fhEDensityCluster ->Fill(eDenClus);
1404 fhEDensityCell ->Fill(eDenCell);
1405 fhEDensityCellvsCluster->Fill(eDenClus, eDenCell);
1406 //Fill the average number of cells or clusters per SM
1407 eClusTot /=fNModules;
1408 eCellTot /=fNModules;
1409 fhAverTotECluster ->Fill(eClusTot);
1410 fhAverTotECell ->Fill(eCellTot);
1411 fhAverTotECellvsCluster->Fill(eClusTot, eCellTot);
1412 //printf("Average Cluster: E %f, density %f; Average Cell E %f, density %f\n ",eClusTot,eDenClus,eCellTot,eDenCell);
1414 // //Average weighted pair distance
1418 // fhAverClusterPairDist ->Fill(rxz );
1419 // fhAverClusterPairDistWeight ->Fill(rxzw);
1420 // fhAverClusterPairDistvsAverE ->Fill(rxz ,eDenClus);
1421 // fhAverClusterPairDistWeightvsAverE->Fill(rxzw,eDenClus);
1422 // fhAverClusterPairDistvsN ->Fill(rxz ,nClus);
1423 // fhAverClusterPairDistWeightvsN ->Fill(rxzw,nClus);
1426 // fhMaxEvsClustEDen->Fill(emax,eDenClus);
1427 // fhMaxEvsClustMult->Fill(emax,nPhot);
1429 //printf("Average Distance: %f; weighted %f\n ",rxz,rxzw);
1435 Double_t vert[] = {0.0, 0.0, 0.0} ; //vertex
1436 Int_t evtIndex1 = 0 ;
1437 Int_t currentEvtIndex = -1;
1438 Int_t curCentrBin = 0 ;
1439 Int_t curRPBin = 0 ;
1440 Int_t curZvertBin = 0 ;
1442 //---------------------------------
1443 //First loop on photons/clusters
1444 //---------------------------------
1445 for(Int_t i1=0; i1<nPhot-1; i1++){
1446 AliAODPWG4Particle * p1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i1)) ;
1447 //printf("AliAnaPi0::MakeAnalysisFillHistograms() : cluster1 id %d\n",p1->GetCaloLabel(0));
1449 // get the event index in the mixed buffer where the photon comes from
1450 // in case of mixing with analysis frame, not own mixing
1451 evtIndex1 = GetEventIndex(p1, vert) ;
1452 //printf("charge = %d\n", track->Charge());
1453 if ( evtIndex1 == -1 )
1455 if ( evtIndex1 == -2 )
1457 if(TMath::Abs(vert[2]) > GetZvertexCut()) continue ; //vertex cut
1459 //----------------------------------------------------------------------------
1460 // Get the multiplicity bin. Different cases: centrality (PbPb),
1461 // average cluster multiplicity, average cell multiplicity, track multiplicity
1462 // default is centrality bins
1463 //----------------------------------------------------------------------------
1464 if (evtIndex1 != currentEvtIndex) {
1465 if(fUseTrackMultBins){ // Track multiplicity bins
1466 //printf("track mult %d\n",GetTrackMultiplicity());
1467 curCentrBin = (GetTrackMultiplicity()-1)/5;
1468 if(curCentrBin > fNCentrBin-1) curCentrBin=fNCentrBin-1;
1469 //printf("track mult bin %d\n",curCentrBin);
1471 else if(fUsePhotonMultBins){ // Photon multiplicity bins
1472 //printf("photon mult %d cluster mult %d\n",nPhot, nClus);
1474 if(curRPBin > GetNRPBin() -1) curRPBin=GetNRPBin()-1;
1475 //printf("photon mult bin %d\n",curRPBin);
1477 else if(fUseAverClusterEBins){ // Cluster average energy bins
1478 //Bins for pp, if needed can be done in a more general way
1479 curCentrBin = eClusTot/10 * fNCentrBin;
1480 if(curCentrBin > fNCentrBin-1) curCentrBin=fNCentrBin-1;
1481 //printf("cluster E average %f, bin %d \n",eClusTot,curCentrBin);
1483 else if(fUseAverCellEBins){ // Cell average energy bins
1484 //Bins for pp, if needed can be done in a more general way
1485 curCentrBin = eCellTot/10*fNCentrBin;
1486 if(curCentrBin > fNCentrBin-1) curCentrBin=fNCentrBin-1;
1487 //printf("cell E average %f, bin %d \n",eCellTot,curCentrBin);
1489 else if(fUseAverClusterEDenBins){ // Energy density bins
1490 //Bins for pp, if needed can be done in a more general way
1491 curCentrBin = eDenClus/10*fNCentrBin;
1492 if(curCentrBin > fNCentrBin-1) curCentrBin=fNCentrBin-1;
1493 //printf("cluster Eden average %f, bin %d \n",eDenClus,curCentrBin);
1495 // else if(fUseAverClusterPairRBins){ // Cluster average distance bins
1496 // //Bins for pp, if needed can be done in a more general way
1497 // curCentrBin = rxz/650*fNCentrBin;
1498 // if(curCentrBin > fNCentrBin-1) curCentrBin=fNCentrBin-1;
1499 // //printf("cluster pair R average %f, bin %d \n",rxz,curCentrBin);
1501 // else if(fUseAverClusterPairRWeightBins){ // Cluster average distance bins
1502 // //Bins for pp, if needed can be done in a more general way
1503 // curCentrBin = rxzw/350*fNCentrBin;
1504 // if(curCentrBin > fNCentrBin-1) curCentrBin=fNCentrBin-1;
1505 // //printf("cluster pair rW average %f, bin %d \n",rxzw,curCentrBin);
1507 // else if(fUseEMaxBins){ // Cluster average distance bins
1508 // //Bins for pp, if needed can be done in a more general way
1509 // curCentrBin = emax/20*fNCentrBin;
1510 // if(curCentrBin > fNCentrBin-1) curCentrBin=fNCentrBin-1;
1511 // //printf("cluster pair rW average %f, bin %d \n",rxzw,curCentrBin);
1513 else { //Event centrality
1514 curCentrBin = GetEventCentrality();
1517 //Reaction plane bin
1521 curZvertBin = (Int_t)(0.5*GetNZvertBin()*(vert[2]+GetZvertexCut())/GetZvertexCut()) ;
1523 //Fill event bin info
1524 fhEvents->Fill(curCentrBin+0.5,curZvertBin+0.5,curRPBin+0.5) ;
1525 currentEvtIndex = evtIndex1 ;
1527 printf("AliAnaPi0::MakeAnalysisFillHistograms() - Centrality %d, Vertex Bin %d, RP bin %d \n",curCentrBin,curRPBin,curZvertBin);
1530 //printf("AliAnaPi0::MakeAnalysisFillHistograms(): Photon 1 Evt %d Vertex : %f,%f,%f\n",evtIndex1, GetVertex(evtIndex1)[0] ,GetVertex(evtIndex1)[1],GetVertex(evtIndex1)[2]);
1532 //Get the momentum of this cluster
1533 TLorentzVector photon1(p1->Px(),p1->Py(),p1->Pz(),p1->E());
1535 //Get (Super)Module number of this cluster
1536 module1 = GetModuleNumber(p1);
1538 //---------------------------------
1539 //Second loop on photons/clusters
1540 //---------------------------------
1541 for(Int_t i2=i1+1; i2<nPhot; i2++){
1542 AliAODPWG4Particle * p2 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i2)) ;
1544 //In case of mixing frame, check we are not in the same event as the first cluster
1545 Int_t evtIndex2 = GetEventIndex(p2, vert) ;
1546 if ( evtIndex2 == -1 )
1548 if ( evtIndex2 == -2 )
1550 if (GetMixedEvent() && (evtIndex1 == evtIndex2))
1553 //printf("AliAnaPi0::MakeAnalysisFillHistograms(): Photon 2 Evt %d Vertex : %f,%f,%f\n",evtIndex2, GetVertex(evtIndex2)[0] ,GetVertex(evtIndex2)[1],GetVertex(evtIndex2)[2]);
1555 //Get the momentum of this cluster
1556 TLorentzVector photon2(p2->Px(),p2->Py(),p2->Pz(),p2->E());
1558 module2 = GetModuleNumber(p2);
1560 //---------------------------------
1561 // Get pair kinematics
1562 //---------------------------------
1563 Double_t m = (photon1 + photon2).M() ;
1564 Double_t pt = (photon1 + photon2).Pt();
1565 Double_t deta = photon1.Eta() - photon2.Eta();
1566 Double_t dphi = photon1.Phi() - photon2.Phi();
1567 Double_t a = TMath::Abs(p1->E()-p2->E())/(p1->E()+p2->E()) ;
1570 printf(" E: photon1 %f, photon2 %f; Pair: pT %f, mass %f, a %f\n", p1->E(), p2->E(), (photon1 + photon2).E(),m,a);
1572 //--------------------------------
1573 // Opening angle selection
1574 //--------------------------------
1575 //Check if opening angle is too large or too small compared to what is expected
1576 Double_t angle = photon1.Angle(photon2.Vect());
1577 if(fUseAngleEDepCut && !GetNeutralMesonSelection()->IsAngleInWindow((photon1+photon2).E(),angle+0.05)) {
1579 printf("AliAnaPi0::MakeAnalysisFillHistograms() -Real pair angle %f not in E %f window\n",angle, (photon1+photon2).E());
1583 if(fUseAngleCut && (angle < fAngleCut || angle > fAngleMaxCut)) {
1585 printf("AliAnaPi0::MakeAnalysisFillHistograms() - Real pair cut %f < angle %f < cut %f\n",fAngleCut, angle, fAngleMaxCut);
1589 //-------------------------------------------------------------------------------------------------
1590 //Fill module dependent histograms, put a cut on assymmetry on the first available cut in the array
1591 //-------------------------------------------------------------------------------------------------
1592 if(a < fAsymCuts[0]){
1593 if(module1==module2 && module1 >=0 && module1<fNModules)
1594 fhReMod[module1]->Fill(pt,m) ;
1596 fhReDiffMod[fNModules+2]->Fill(pt,m) ;
1598 if(fCalorimeter=="EMCAL"){
1599 if((module1==0 && module2==2) || (module1==2 && module2==0)) fhReDiffMod[0]->Fill(pt,m) ;
1600 if((module1==1 && module2==3) || (module1==3 && module2==1)) fhReDiffMod[1]->Fill(pt,m) ;
1601 if((module1==0 && module2==1) || (module1==1 && module2==0)) fhReDiffMod[2]->Fill(pt,m) ;
1602 if((module1==2 && module2==3) || (module1==3 && module2==2)) fhReDiffMod[3]->Fill(pt,m) ;
1603 if((module1==0 && module2==3) || (module1==3 && module2==0)) fhReDiffMod[4]->Fill(pt,m) ;
1604 if((module1==2 && module2==1) || (module1==1 && module2==2)) fhReDiffMod[5]->Fill(pt,m) ;
1607 if((module1==0 && module2==1) || (module1==1 && module2==0)) fhReDiffMod[0]->Fill(pt,m) ;
1608 if((module1==0 && module2==2) || (module1==2 && module2==0)) fhReDiffMod[1]->Fill(pt,m) ;
1609 if((module1==1 && module2==2) || (module1==2 && module2==1)) fhReDiffMod[2]->Fill(pt,m) ;
1613 //In case we want only pairs in same (super) module, check their origin.
1615 if(fSameSM && module1!=module2) ok=kFALSE;
1618 //Check if one of the clusters comes from a conversion
1619 if (p1->IsTagged() && p2->IsTagged()) fhReConv2->Fill(pt,m);
1620 else if(p1->IsTagged() || p2->IsTagged()) fhReConv ->Fill(pt,m);
1622 //Fill histograms for different bad channel distance, centrality, assymmetry cut and pid bit
1623 for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
1624 if((p1->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton)) && (p2->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton))){
1625 for(Int_t iasym=0; iasym < fNAsymCuts; iasym++){
1626 if(a < fAsymCuts[iasym]){
1627 Int_t index = ((curCentrBin*fNPIDBits)+ipid)*fNAsymCuts + iasym;
1628 //printf("index %d :(cen %d * nPID %d + ipid %d)*nasym %d + iasym %d\n",index,curCentrBin,fNPIDBits,ipid,fNAsymCuts,iasym);
1629 fhRe1 [index]->Fill(pt,m);
1630 if(fMakeInvPtPlots)fhReInvPt1[index]->Fill(pt,m,1./pt) ;
1631 if(fFillBadDistHisto){
1632 if(p1->DistToBad()>0 && p2->DistToBad()>0){
1633 fhRe2 [index]->Fill(pt,m) ;
1634 if(fMakeInvPtPlots)fhReInvPt2[index]->Fill(pt,m,1./pt) ;
1635 if(p1->DistToBad()>1 && p2->DistToBad()>1){
1636 fhRe3 [index]->Fill(pt,m) ;
1637 if(fMakeInvPtPlots)fhReInvPt3[index]->Fill(pt,m,1./pt) ;
1640 }// Fill bad dist histos
1642 }// asymmetry cut loop
1646 //Fill histograms with opening angle
1647 fhRealOpeningAngle ->Fill(pt,angle);
1648 fhRealCosOpeningAngle->Fill(pt,TMath::Cos(angle));
1650 //Fill histograms with pair assymmetry
1651 fhRePtAsym->Fill(pt,a);
1652 if(m > 0.10 && m < 0.17) fhRePtAsymPi0->Fill(pt,a);
1653 if(m > 0.45 && m < 0.65) fhRePtAsymEta->Fill(pt,a);
1655 //-------------------------------------------------------
1656 //Get the number of cells needed for multi cut analysis.
1657 //-------------------------------------------------------
1660 if(fMultiCutAna || (IsDataMC() && fMultiCutAnaSim)){
1662 AliVEvent * event = GetReader()->GetInputEvent();
1664 for(Int_t iclus = 0; iclus < event->GetNumberOfCaloClusters(); iclus++){
1665 AliVCluster *cluster = event->GetCaloCluster(iclus);
1668 if (fCalorimeter == "EMCAL" && GetReader()->IsEMCALCluster(cluster)) is = kTRUE;
1669 else if(fCalorimeter == "PHOS" && GetReader()->IsPHOSCluster (cluster)) is = kTRUE;
1672 if (p1->GetCaloLabel(0) == cluster->GetID()) ncell1 = cluster->GetNCells();
1673 else if (p2->GetCaloLabel(0) == cluster->GetID()) ncell2 = cluster->GetNCells();
1674 } // PHOS or EMCAL cluster as requested in analysis
1676 if(ncell2 > 0 && ncell1 > 0) break; // No need to continue the iteration
1679 //printf("e 1: %2.2f, e 2: %2.2f, ncells: n1 %d, n2 %d\n", p1->E(), p2->E(),ncell1,ncell2);
1686 //Do some MC checks on the origin of the pair, is there any common ancestor and if there is one, who?
1687 if(IsDataMC()) FillMCVersusRecDataHistograms(p1->GetLabel(), p2->GetLabel(),p1->Pt(), p2->Pt(),ncell1, ncell2, m, pt, a,deta, dphi);
1689 //-----------------------
1690 //Multi cuts analysis
1691 //-----------------------
1693 //Histograms for different PID bits selection
1694 for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
1696 if(p1->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton) &&
1697 p2->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton)) fhRePIDBits[ipid]->Fill(pt,m) ;
1699 //printf("ipt %d, ipid%d, name %s\n",ipt, ipid, fhRePtPIDCuts[ipt*fNPIDBitsBits+ipid]->GetName());
1700 } // pid bit cut loop
1702 //Several pt,ncell and asymmetry cuts
1703 for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
1704 for(Int_t icell=0; icell<fNCellNCuts; icell++){
1705 for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
1706 Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
1707 if(p1->Pt() > fPtCuts[ipt] && p2->Pt() > fPtCuts[ipt] &&
1708 a < fAsymCuts[iasym] &&
1709 ncell1 >= fCellNCuts[icell] && ncell2 >= fCellNCuts[icell]){
1710 fhRePtNCellAsymCuts[index]->Fill(pt,m) ;
1711 //printf("ipt %d, icell%d, iasym %d, name %s\n",ipt, icell, iasym, fhRePtNCellAsymCuts[((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym]->GetName());
1712 if(module1==module2){
1713 if (module1==0) fhRePtNCellAsymCutsSM0[index]->Fill(pt,m) ;
1714 else if(module1==1) fhRePtNCellAsymCutsSM1[index]->Fill(pt,m) ;
1715 else if(module1==2) fhRePtNCellAsymCutsSM2[index]->Fill(pt,m) ;
1716 else if(module1==3) fhRePtNCellAsymCutsSM3[index]->Fill(pt,m) ;
1717 else printf("AliAnaPi0::FillHistograms() - WRONG SM NUMBER\n");
1720 }// pid bit cut loop
1723 for(Int_t iasym = 0; iasym < fNAsymCuts; iasym++){
1724 if(a < fAsymCuts[iasym])fhRePtMult[iasym]->Fill(pt,GetTrackMultiplicity(),m) ;
1726 }// multiple cuts analysis
1728 }// second same event particle
1731 //-------------------------------------------------------------
1733 //-------------------------------------------------------------
1735 //printf("Cen bin %d, RP bin %d, e aver %f, mult %d\n",curCentrBin,curRPBin, eClusTot, nPhot);
1736 //Recover events in with same characteristics as the current event
1737 TList * evMixList=fEventsList[curCentrBin*GetNZvertBin()*GetNRPBin()+curZvertBin*GetNRPBin()+curRPBin] ;
1738 Int_t nMixed = evMixList->GetSize() ;
1739 for(Int_t ii=0; ii<nMixed; ii++){
1740 TClonesArray* ev2= (TClonesArray*) (evMixList->At(ii));
1741 Int_t nPhot2=ev2->GetEntriesFast() ;
1744 printf("AliAnaPi0::MakeAnalysisFillHistograms() - Mixed event %d photon entries %d, centrality bin %d\n", ii, nPhot2, curCentrBin);
1746 //---------------------------------
1747 //First loop on photons/clusters
1748 //---------------------------------
1749 for(Int_t i1=0; i1<nPhot; i1++){
1750 AliAODPWG4Particle * p1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i1)) ;
1751 if(fSameSM && GetModuleNumber(p1)!=module1) continue;
1753 //Get kinematics of cluster and (super) module of this cluster
1754 TLorentzVector photon1(p1->Px(),p1->Py(),p1->Pz(),p1->E());
1755 module1 = GetModuleNumber(p1);
1757 //---------------------------------
1758 //First loop on photons/clusters
1759 //---------------------------------
1760 for(Int_t i2=0; i2<nPhot2; i2++){
1761 AliAODPWG4Particle * p2 = (AliAODPWG4Particle*) (ev2->At(i2)) ;
1763 //Get kinematics of second cluster and calculate those of the pair
1764 TLorentzVector photon2(p2->Px(),p2->Py(),p2->Pz(),p2->E());
1765 m = (photon1+photon2).M() ;
1766 Double_t pt = (photon1 + photon2).Pt();
1767 Double_t a = TMath::Abs(p1->E()-p2->E())/(p1->E()+p2->E()) ;
1769 //Check if opening angle is too large or too small compared to what is expected
1770 Double_t angle = photon1.Angle(photon2.Vect());
1771 if(fUseAngleEDepCut && !GetNeutralMesonSelection()->IsAngleInWindow((photon1+photon2).E(),angle+0.05)){
1773 printf("AliAnaPi0::MakeAnalysisFillHistograms() -Mix pair angle %f not in E %f window\n",angle, (photon1+photon2).E());
1776 if(fUseAngleCut && (angle < fAngleCut || angle > fAngleMaxCut)) {
1778 printf("AliAnaPi0::MakeAnalysisFillHistograms() -Mix pair angle %f < cut %f\n",angle,fAngleCut);
1784 printf("AliAnaPi0::MakeAnalysisFillHistograms() - Mixed Event: pT: photon1 %2.2f, photon2 %2.2f; Pair: pT %2.2f, mass %2.3f, a %f2.3\n",
1785 p1->Pt(), p2->Pt(), pt,m,a);
1787 //In case we want only pairs in same (super) module, check their origin.
1788 module2 = GetModuleNumber(p2);
1790 //-------------------------------------------------------------------------------------------------
1791 //Fill module dependent histograms, put a cut on assymmetry on the first available cut in the array
1792 //-------------------------------------------------------------------------------------------------
1793 if(a < fAsymCuts[0]){
1794 if(module1==module2 && module1 >=0 && module1<fNModules)
1795 fhMiMod[module1]->Fill(pt,m) ;
1797 fhMiDiffMod[fNModules+2]->Fill(pt,m) ;
1799 if(fCalorimeter=="EMCAL"){
1800 if((module1==0 && module2==2) || (module1==2 && module2==0)) fhMiDiffMod[0]->Fill(pt,m) ;
1801 if((module1==1 && module2==3) || (module1==3 && module2==1)) fhMiDiffMod[1]->Fill(pt,m) ;
1802 if((module1==0 && module2==1) || (module1==1 && module2==0)) fhMiDiffMod[2]->Fill(pt,m) ;
1803 if((module1==2 && module2==3) || (module1==3 && module2==2)) fhMiDiffMod[3]->Fill(pt,m) ;
1804 if((module1==0 && module2==3) || (module1==3 && module2==0)) fhMiDiffMod[4]->Fill(pt,m) ;
1805 if((module1==2 && module2==1) || (module1==1 && module2==2)) fhMiDiffMod[5]->Fill(pt,m) ;
1809 if((module1==0 && module2==1) || (module1==1 && module2==0)) fhMiDiffMod[0]->Fill(pt,m) ;
1810 if((module1==0 && module2==2) || (module1==2 && module2==0)) fhMiDiffMod[1]->Fill(pt,m) ;
1811 if((module1==1 && module2==2) || (module1==2 && module2==1)) fhMiDiffMod[2]->Fill(pt,m) ;
1816 if(fSameSM && module1!=module2) ok=kFALSE;
1819 //Check if one of the clusters comes from a conversion
1820 if (p1->IsTagged() && p2->IsTagged()) fhMiConv2->Fill(pt,m);
1821 else if(p1->IsTagged() || p2->IsTagged()) fhMiConv ->Fill(pt,m);
1823 //Fill histograms for different bad channel distance, centrality, assymmetry cut and pid bit
1824 for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
1825 if((p1->IsPIDOK(ipid,AliCaloPID::kPhoton)) && (p2->IsPIDOK(ipid,AliCaloPID::kPhoton))){
1826 for(Int_t iasym=0; iasym < fNAsymCuts; iasym++){
1827 if(a < fAsymCuts[iasym]){
1828 Int_t index = ((curCentrBin*fNPIDBits)+ipid)*fNAsymCuts + iasym;
1829 fhMi1 [index]->Fill(pt,m) ;
1830 if(fMakeInvPtPlots)fhMiInvPt1[index]->Fill(pt,m,1./pt) ;
1831 if(fFillBadDistHisto){
1832 if(p1->DistToBad()>0 && p2->DistToBad()>0){
1833 fhMi2 [index]->Fill(pt,m) ;
1834 if(fMakeInvPtPlots)fhMiInvPt2[index]->Fill(pt,m,1./pt) ;
1835 if(p1->DistToBad()>1 && p2->DistToBad()>1){
1836 fhMi3 [index]->Fill(pt,m) ;
1837 if(fMakeInvPtPlots)fhMiInvPt3[index]->Fill(pt,m,1./pt) ;
1840 }// Fill bad dist histo
1844 }//loop for histograms
1846 //-----------------------
1847 //Multi cuts analysis
1848 //-----------------------
1850 //Several pt,ncell and asymmetry cuts
1851 for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
1852 for(Int_t icell=0; icell<fNCellNCuts; icell++){
1853 for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
1854 Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
1855 if(p1->Pt() > fPtCuts[ipt] && p2->Pt() > fPtCuts[ipt] &&
1856 a < fAsymCuts[iasym] &&
1857 p1->GetBtag() >= fCellNCuts[icell] && p2->GetBtag() >= fCellNCuts[icell]){
1858 fhMiPtNCellAsymCuts[index]->Fill(pt,m) ;
1859 //printf("ipt %d, icell%d, iasym %d, name %s\n",ipt, icell, iasym, fhRePtNCellAsymCuts[((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym]->GetName());
1861 }// pid bit cut loop
1866 //Fill histograms with opening angle
1867 fhMixedOpeningAngle ->Fill(pt,angle);
1868 fhMixedCosOpeningAngle->Fill(pt,TMath::Cos(angle));
1870 }// second cluster loop
1871 }//first cluster loop
1872 }//loop on mixed events
1874 //--------------------------------------------------------
1875 //Add the current event to the list of events for mixing
1876 //--------------------------------------------------------
1877 TClonesArray *currentEvent = new TClonesArray(*GetInputAODBranch());
1878 //Add current event to buffer and Remove redundant events
1879 if(currentEvent->GetEntriesFast()>0){
1880 evMixList->AddFirst(currentEvent) ;
1881 currentEvent=0 ; //Now list of particles belongs to buffer and it will be deleted with buffer
1882 if(evMixList->GetSize()>=fNmaxMixEv)
1884 TClonesArray * tmp = (TClonesArray*) (evMixList->Last()) ;
1885 evMixList->RemoveLast() ;
1890 delete currentEvent ;
1897 //________________________________________________________________________
1898 void AliAnaPi0::ReadHistograms(TList* outputList)
1900 // Needed when Terminate is executed in distributed environment
1901 // Refill analysis histograms of this class with corresponding histograms in output list.
1903 // Histograms of this analsys are kept in the same list as other analysis, recover the position of
1904 // the first one and then add the next.
1905 Int_t index = outputList->IndexOf(outputList->FindObject(GetAddedHistogramsStringToName()+"hRe_cen0_pid0_dist1"));
1907 if(!fhRe1) fhRe1 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
1908 if(!fhRe2) fhRe2 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
1909 if(!fhRe3) fhRe3 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
1910 if(!fhMi1) fhMi1 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
1911 if(!fhMi2) fhMi2 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
1912 if(!fhMi3) fhMi3 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
1913 if(!fhReInvPt1) fhReInvPt1 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
1914 if(!fhReInvPt2) fhReInvPt2 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
1915 if(!fhReInvPt3) fhReInvPt3 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
1916 if(!fhMiInvPt1) fhMiInvPt1 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
1917 if(!fhMiInvPt2) fhMiInvPt2 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
1918 if(!fhMiInvPt3) fhMiInvPt3 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
1919 if(!fhReMod) fhReMod = new TH2D*[fNModules] ;
1920 if(!fhReDiffMod)fhReDiffMod = new TH2D*[fNModules+1] ;
1921 if(!fhMiMod) fhMiMod = new TH2D*[fNModules] ;
1922 if(!fhMiDiffMod)fhMiDiffMod = new TH2D*[fNModules+1] ;
1924 fhReConv = (TH2D*) outputList->At(index++);
1925 fhMiConv = (TH2D*) outputList->At(index++);
1926 fhReConv2 = (TH2D*) outputList->At(index++);
1927 fhMiConv2 = (TH2D*) outputList->At(index++);
1929 for(Int_t ic=0; ic<fNCentrBin; ic++){
1930 for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
1931 for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
1932 Int_t ihisto = ((ic*fNPIDBits)+ipid)*fNAsymCuts + iasym;
1934 fhRe1[ihisto] = (TH2D*) outputList->At(index++);
1935 fhRe2[ihisto] = (TH2D*) outputList->At(index++);
1936 fhRe3[ihisto] = (TH2D*) outputList->At(index++);
1938 fhReInvPt1[ihisto] = (TH2D*) outputList->At(index++);
1939 fhReInvPt2[ihisto] = (TH2D*) outputList->At(index++);
1940 fhReInvPt3[ihisto] = (TH2D*) outputList->At(index++);
1943 fhMi1[ihisto] = (TH2D*) outputList->At(index++);
1944 fhMi2[ihisto] = (TH2D*) outputList->At(index++);
1945 fhMi3[ihisto] = (TH2D*) outputList->At(index++);
1947 fhMiInvPt1[ihisto] = (TH2D*) outputList->At(index++);
1948 fhMiInvPt2[ihisto] = (TH2D*) outputList->At(index++);
1949 fhMiInvPt3[ihisto] = (TH2D*) outputList->At(index++);
1955 fhRePtAsym = (TH2D*)outputList->At(index++);
1956 fhRePtAsymPi0 = (TH2D*)outputList->At(index++);
1957 fhRePtAsymEta = (TH2D*)outputList->At(index++);
1961 if(!fhRePtNCellAsymCuts) fhRePtNCellAsymCuts = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
1962 if(!fhRePIDBits) fhRePIDBits = new TH2D*[fNPIDBits];
1964 for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
1965 fhRePIDBits[ipid] = (TH2D*) outputList->At(index++);
1968 for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
1969 for(Int_t icell=0; icell<fNCellNCuts; icell++){
1970 for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
1971 fhRePtNCellAsymCuts[((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym] = (TH2D*) outputList->At(index++);
1976 if(!fhRePtMult) fhRePtMult = new TH3D*[fNAsymCuts] ;
1977 for(Int_t iasym = 0; iasym < fNAsymCuts; iasym++)
1978 fhRePtMult[iasym] = (TH3D*) outputList->At(index++);
1979 }// multi cut analysis
1981 fhEvents = (TH3D *) outputList->At(index++);
1983 fhRealOpeningAngle = (TH2D*) outputList->At(index++);
1984 fhRealCosOpeningAngle = (TH2D*) outputList->At(index++);
1986 fhMixedOpeningAngle = (TH2D*) outputList->At(index++);
1987 fhMixedCosOpeningAngle = (TH2D*) outputList->At(index++);
1990 //Histograms filled only if MC data is requested
1991 if(IsDataMC() || (GetReader()->GetDataType() == AliCaloTrackReader::kMC) ){
1992 fhPrimPi0Pt = (TH1D*) outputList->At(index++);
1993 fhPrimPi0AccPt = (TH1D*) outputList->At(index++);
1994 fhPrimPi0Y = (TH2D*) outputList->At(index++);
1995 fhPrimPi0AccY = (TH2D*) outputList->At(index++);
1996 fhPrimPi0Phi = (TH2D*) outputList->At(index++);
1997 fhPrimPi0AccPhi = (TH2D*) outputList->At(index++);
1998 fhPrimEtaPt = (TH1D*) outputList->At(index++);
1999 fhPrimEtaAccPt = (TH1D*) outputList->At(index++);
2000 fhPrimEtaY = (TH2D*) outputList->At(index++);
2001 fhPrimEtaAccY = (TH2D*) outputList->At(index++);
2002 fhPrimEtaPhi = (TH2D*) outputList->At(index++);
2003 fhPrimEtaAccPhi = (TH2D*) outputList->At(index++);
2004 for(Int_t i = 0; i<13; i++){
2005 fhMCOrgMass[i] = (TH2D*) outputList->At(index++);
2006 fhMCOrgAsym[i] = (TH2D*) outputList->At(index++);
2007 fhMCOrgDeltaEta[i] = (TH2D*) outputList->At(index++);
2008 fhMCOrgDeltaPhi[i] = (TH2D*) outputList->At(index++);
2011 if(fMultiCutAnaSim){
2012 fhMCPi0MassPtTrue = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
2013 fhMCPi0MassPtRec = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
2014 fhMCPi0PtTruePtRec = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
2015 fhMCEtaMassPtTrue = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
2016 fhMCEtaMassPtRec = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
2017 fhMCEtaPtTruePtRec = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
2018 for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
2019 for(Int_t icell=0; icell<fNCellNCuts; icell++){
2020 for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
2021 Int_t in = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
2022 fhMCPi0MassPtTrue[in] = (TH2D*) outputList->At(index++);
2023 fhMCPi0PtTruePtRec[in] = (TH2D*) outputList->At(index++);
2024 fhMCEtaMassPtTrue[in] = (TH2D*) outputList->At(index++);
2025 fhMCEtaPtTruePtRec[in] = (TH2D*) outputList->At(index++);
2031 fhMCPi0MassPtTrue = new TH2D*[1];
2032 fhMCPi0PtTruePtRec = new TH2D*[1];
2033 fhMCEtaMassPtTrue = new TH2D*[1];
2034 fhMCEtaPtTruePtRec = new TH2D*[1];
2036 fhMCPi0MassPtTrue[0] = (TH2D*) outputList->At(index++);
2037 fhMCPi0PtTruePtRec[0] = (TH2D*) outputList->At(index++);
2038 fhMCEtaMassPtTrue[0] = (TH2D*) outputList->At(index++);
2039 fhMCEtaPtTruePtRec[0] = (TH2D*) outputList->At(index++);
2043 for(Int_t imod=0; imod < fNModules; imod++){
2044 fhReMod[imod] = (TH2D*) outputList->At(index++);
2045 fhReDiffMod[imod] = (TH2D*) outputList->At(index++);
2047 fhMiMod[imod] = (TH2D*) outputList->At(index++);
2048 fhMiDiffMod[imod] = (TH2D*) outputList->At(index++);
2055 //____________________________________________________________________________________________________________________________________________________
2056 void AliAnaPi0::Terminate(TList* outputList)
2058 //Do some calculations and plots from the final histograms.
2060 printf(" *** %s Terminate:\n", GetName()) ;
2062 //Recover histograms from output histograms list, needed for distributed analysis.
2063 ReadHistograms(outputList);
2066 printf("AliAnaPi0::Terminate() - Error: Remote output histograms not imported in AliAnaPi0 object");
2070 printf("AliAnaPi0::Terminate() Mgg Real : %5.3f , RMS : %5.3f \n", fhRe1[0]->GetMean(), fhRe1[0]->GetRMS() ) ;
2072 const Int_t buffersize = 255;
2074 char nameIM[buffersize];
2075 snprintf(nameIM, buffersize,"AliAnaPi0_%s_cPt",fCalorimeter.Data());
2076 TCanvas * cIM = new TCanvas(nameIM, "", 400, 10, 600, 700) ;
2081 TH1D * hIMAllPt = (TH1D*) fhRe1[0]->ProjectionY(Form("IMPtAll_%s",fCalorimeter.Data()));
2082 hIMAllPt->SetLineColor(2);
2083 hIMAllPt->SetTitle("No cut on p_{T, #gamma#gamma} ");
2087 TH1D * hIMPt5 = (TH1D*) fhRe1[0]->ProjectionY(Form("IMPt0-5_%s",fCalorimeter.Data()),0, fhRe1[0]->GetXaxis()->FindBin(5.));
2088 // hRe1Pt5->GetXaxis()->SetRangeUser(0,5);
2089 // TH1D * hIMPt5 = (TH1D*) hRe1Pt5->Project3D(Form("IMPt5_%s_pz",fCalorimeter.Data()));
2090 hIMPt5->SetLineColor(2);
2091 hIMPt5->SetTitle("0 < p_{T, #gamma#gamma} < 5 GeV/c");
2095 TH1D * hIMPt10 = (TH1D*) fhRe1[0]->ProjectionY(Form("IMPt5-10_%s",fCalorimeter.Data()), fhRe1[0]->GetXaxis()->FindBin(5.),fhRe1[0]->GetXaxis()->FindBin(10.));
2096 // hRe1Pt10->GetXaxis()->SetRangeUser(5,10);
2097 // TH1D * hIMPt10 = (TH1D*) hRe1Pt10->Project3D(Form("IMPt10_%s_pz",fCalorimeter.Data()));
2098 hIMPt10->SetLineColor(2);
2099 hIMPt10->SetTitle("5 < p_{T, #gamma#gamma} < 10 GeV/c");
2103 TH1D * hIMPt20 = (TH1D*) fhRe1[0]->ProjectionY(Form("IMPt10-20_%s",fCalorimeter.Data()), fhRe1[0]->GetXaxis()->FindBin(10.),fhRe1[0]->GetXaxis()->FindBin(20.));
2104 // TH3F * hRe1Pt20 = (TH3F*)fhRe1[0]->Clone(Form("IMPt20_%s",fCalorimeter.Data()));
2105 // hRe1Pt20->GetXaxis()->SetRangeUser(10,20);
2106 // TH1D * hIMPt20 = (TH1D*) hRe1Pt20->Project3D(Form("IMPt20_%s_pz",fCalorimeter.Data()));
2107 hIMPt20->SetLineColor(2);
2108 hIMPt20->SetTitle("10 < p_{T, #gamma#gamma} < 20 GeV/c");
2111 char nameIMF[buffersize];
2112 snprintf(nameIMF,buffersize,"AliAnaPi0_%s_Mgg.eps",fCalorimeter.Data());
2113 cIM->Print(nameIMF);
2115 char namePt[buffersize];
2116 snprintf(namePt,buffersize,"AliAnaPi0_%s_cPt",fCalorimeter.Data());
2117 TCanvas * cPt = new TCanvas(namePt, "", 400, 10, 600, 700) ;
2122 TH1D * hPt = (TH1D*) fhRe1[0]->ProjectionX(Form("Pt0_%s",fCalorimeter.Data()),-1,-1);
2123 hPt->SetLineColor(2);
2124 hPt->SetTitle("No cut on M_{#gamma#gamma} ");
2128 TH1D * hPtIM1 = (TH1D*)fhRe1[0]->ProjectionX(Form("Pt1_%s",fCalorimeter.Data()), fhRe1[0]->GetZaxis()->FindBin(0.05),fhRe1[0]->GetZaxis()->FindBin(0.21));
2129 // TH3F * hRe1IM1 = (TH3F*)fhRe1[0]->Clone(Form("Pt1_%s",fCalorimeter.Data()));
2130 // hRe1IM1->GetZaxis()->SetRangeUser(0.05,0.21);
2131 // TH1D * hPtIM1 = (TH1D*) hRe1IM1->Project3D("x");
2132 hPtIM1->SetLineColor(2);
2133 hPtIM1->SetTitle("0.05 < M_{#gamma#gamma} < 0.21 GeV/c^{2}");
2137 TH1D * hPtIM2 = (TH1D*)fhRe1[0]->ProjectionX(Form("Pt2_%s",fCalorimeter.Data()), fhRe1[0]->GetZaxis()->FindBin(0.09),fhRe1[0]->GetZaxis()->FindBin(0.17));
2138 // TH3F * hRe1IM2 = (TH3F*)fhRe1[0]->Clone(Form("Pt2_%s",fCalorimeter.Data()));
2139 // hRe1IM2->GetZaxis()->SetRangeUser(0.09,0.17);
2140 // TH1D * hPtIM2 = (TH1D*) hRe1IM2->Project3D("x");
2141 hPtIM2->SetLineColor(2);
2142 hPtIM2->SetTitle("0.09 < M_{#gamma#gamma} < 0.17 GeV/c^{2}");
2146 TH1D * hPtIM3 = (TH1D*)fhRe1[0]->ProjectionX(Form("Pt3_%s",fCalorimeter.Data()), fhRe1[0]->GetZaxis()->FindBin(0.11),fhRe1[0]->GetZaxis()->FindBin(0.15));
2147 // TH3F * hRe1IM3 = (TH3F*)fhRe1[0]->Clone(Form("Pt3_%s",fCalorimeter.Data()));
2148 // hRe1IM3->GetZaxis()->SetRangeUser(0.11,0.15);
2149 // TH1D * hPtIM3 = (TH1D*) hRe1IM1->Project3D("x");
2150 hPtIM3->SetLineColor(2);
2151 hPtIM3->SetTitle("0.11 < M_{#gamma#gamma} < 0.15 GeV/c^{2}");
2154 char namePtF[buffersize];
2155 snprintf(namePtF,buffersize,"AliAnaPi0_%s_Pt.eps",fCalorimeter.Data());
2156 cPt->Print(namePtF);
2158 char line[buffersize] ;
2159 snprintf(line,buffersize,".!tar -zcf %s_%s.tar.gz *.eps", GetName(),fCalorimeter.Data()) ;
2160 gROOT->ProcessLine(line);
2161 snprintf(line, buffersize,".!rm -fR AliAnaPi0_%s*.eps",fCalorimeter.Data());
2162 gROOT->ProcessLine(line);
2164 printf(" AliAnaPi0::Terminate() - !! All the eps files are in %s_%s.tar.gz !!!\n", GetName(), fCalorimeter.Data());
2167 //____________________________________________________________________________________________________________________________________________________
2168 Int_t AliAnaPi0::GetEventIndex(AliAODPWG4Particle * part, Double_t * vert)
2170 // retieves the event index and checks the vertex
2171 // in the mixed buffer returns -2 if vertex NOK
2172 // for normal events returns 0 if vertex OK and -1 if vertex NOK
2174 Int_t evtIndex = -1 ;
2175 if(GetReader()->GetDataType()!=AliCaloTrackReader::kMC){
2177 if (GetMixedEvent()){
2179 evtIndex = GetMixedEvent()->EventIndexForCaloCluster(part->GetCaloLabel(0)) ;
2180 GetVertex(vert,evtIndex);
2182 if(TMath::Abs(vert[2])> GetZvertexCut())
2183 evtIndex = -2 ; //Event can not be used (vertex, centrality,... cuts not fulfilled)
2184 } else {// Single event
2188 if(TMath::Abs(vert[2])> GetZvertexCut())
2189 evtIndex = -1 ; //Event can not be used (vertex, centrality,... cuts not fulfilled)