]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG4/PartCorrDep/AliAnaPi0.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Correct the way primary photons are in EMCAL acceptance, change some histograms title...
[u/mrichter/AliRoot.git] / PWG4 / PartCorrDep / AliAnaPi0.cxx
CommitLineData
1c5acb87 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15/* $Id: $ */
16
17//_________________________________________________________________________
18// Class to collect two-photon invariant mass distributions for
6175da48 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
1c5acb87 22//
23//-- Author: Dmitri Peressounko (RRC "KI")
24//-- Adapted to PartCorr frame by Lamia Benhabib (SUBATECH)
25//-- and Gustavo Conesa (INFN-Frascati)
26//_________________________________________________________________________
27
28
29// --- ROOT system ---
30#include "TH3.h"
50f39b97 31#include "TH2D.h"
1c5acb87 32//#include "Riostream.h"
6639984f 33#include "TCanvas.h"
34#include "TPad.h"
35#include "TROOT.h"
477d6cee 36#include "TClonesArray.h"
0c1383b5 37#include "TObjString.h"
6175da48 38#include "TDatabasePDG.h"
1c5acb87 39
40//---- AliRoot system ----
41#include "AliAnaPi0.h"
42#include "AliCaloTrackReader.h"
43#include "AliCaloPID.h"
6639984f 44#include "AliStack.h"
ff45398a 45#include "AliFiducialCut.h"
477d6cee 46#include "TParticle.h"
477d6cee 47#include "AliVEvent.h"
6921fa00 48#include "AliESDCaloCluster.h"
49#include "AliESDEvent.h"
50#include "AliAODEvent.h"
50f39b97 51#include "AliNeutralMesonSelection.h"
c8fe2783 52#include "AliMixedEvent.h"
6175da48 53#include "AliAODMCParticle.h"
591cc579 54
1c5acb87 55ClassImp(AliAnaPi0)
56
57//________________________________________________________________________________________________________________________________________________
58AliAnaPi0::AliAnaPi0() : AliAnaPartCorrBaseClass(),
5025c139 59fDoOwnMix(kFALSE),fNCentrBin(0),//fNZvertBin(0),fNrpBin(0),
af7b3903 60fNmaxMixEv(0), fCalorimeter(""),
6175da48 61fNModules(12), fUseAngleCut(kFALSE), fUseAngleEDepCut(kFALSE),fAngleCut(0), fAngleMaxCut(7.),fEventsList(0x0), fMultiCutAna(kFALSE), fMultiCutAnaSim(kFALSE),
62fNPtCuts(0),fNAsymCuts(0), fNCellNCuts(0),fNPIDBits(0), fMakeInvPtPlots(kFALSE), fSameSM(kFALSE),
156549ae 63fUseTrackMultBins(kFALSE),fUsePhotonMultBins(kFALSE),fUseAverCellEBins(kFALSE), fUseAverClusterEBins(kFALSE),
64fUseAverClusterEDenBins(0), //fUseAverClusterPairRBins(0), fUseAverClusterPairRWeightBins(0), fUseEMaxBins(0),
65fFillBadDistHisto(kFALSE),
66fhAverTotECluster(0), fhAverTotECell(0), fhAverTotECellvsCluster(0),
67fhEDensityCluster(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),
71fhReMod(0x0), fhReDiffMod(0x0), fhMiMod(0x0), fhMiDiffMod(0x0),
72fhReConv(0x0), fhMiConv(0x0), fhReConv2(0x0), fhMiConv2(0x0),
73fhRe1(0x0), fhMi1(0x0), fhRe2(0x0), fhMi2(0x0), fhRe3(0x0), fhMi3(0x0),
74fhReInvPt1(0x0), fhMiInvPt1(0x0), fhReInvPt2(0x0), fhMiInvPt2(0x0), fhReInvPt3(0x0), fhMiInvPt3(0x0),
6175da48 75fhRePtNCellAsymCuts(0x0), fhRePtNCellAsymCutsSM0(0x0), fhRePtNCellAsymCutsSM1(0x0), fhRePtNCellAsymCutsSM2(0x0), fhRePtNCellAsymCutsSM3(0x0), fhMiPtNCellAsymCuts(0x0),
76fhRePIDBits(0x0),fhRePtMult(0x0), fhRePtAsym(0x0), fhRePtAsymPi0(0x0),fhRePtAsymEta(0x0),
77fhEvents(0x0), fhRealOpeningAngle(0x0),fhRealCosOpeningAngle(0x0), fhMixedOpeningAngle(0x0),fhMixedCosOpeningAngle(0x0),
78fhPrimPi0Pt(0x0), fhPrimPi0AccPt(0x0), fhPrimPi0Y(0x0), fhPrimPi0AccY(0x0), fhPrimPi0Phi(0x0), fhPrimPi0AccPhi(0x0),
79fhPrimPi0OpeningAngle(0x0),fhPrimPi0CosOpeningAngle(0x0),
80fhPrimEtaPt(0x0), fhPrimEtaAccPt(0x0), fhPrimEtaY(0x0), fhPrimEtaAccY(0x0), fhPrimEtaPhi(0x0), fhPrimEtaAccPhi(0x0),
81fhMCOrgMass(),fhMCOrgAsym(), fhMCOrgDeltaEta(),fhMCOrgDeltaPhi(),
82fhMCPi0MassPtRec(), fhMCPi0MassPtTrue(), fhMCPi0PtTruePtRec(), fhMCEtaMassPtRec(), fhMCEtaMassPtTrue(), fhMCEtaPtTruePtRec()
1c5acb87 83{
84//Default Ctor
85 InitParameters();
86
87}
1c5acb87 88
89//________________________________________________________________________________________________________________________________________________
90AliAnaPi0::~AliAnaPi0() {
477d6cee 91 // Remove event containers
7e7694bb 92
93 if(fDoOwnMix && fEventsList){
477d6cee 94 for(Int_t ic=0; ic<fNCentrBin; ic++){
5025c139 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] ;
7e7694bb 99 }
477d6cee 100 }
101 }
102 delete[] fEventsList;
103 fEventsList=0 ;
104 }
591cc579 105
1c5acb87 106}
107
108//________________________________________________________________________________________________________________________________________________
109void AliAnaPi0::InitParameters()
110{
111//Init parameters when first called the analysis
112//Set default parameters
a3aebfff 113 SetInputAODName("PWG4Particle");
114
115 AddToHistogramsName("AnaPi0_");
6921fa00 116 fNModules = 12; // set maximum to maximum number of EMCAL modules
477d6cee 117 fNCentrBin = 1;
5025c139 118// fNZvertBin = 1;
119// fNrpBin = 1;
477d6cee 120 fNmaxMixEv = 10;
5025c139 121
477d6cee 122 fCalorimeter = "PHOS";
50f39b97 123 fUseAngleCut = kFALSE;
6175da48 124 fUseAngleEDepCut = kFALSE;
125 fAngleCut = 0.;
126 fAngleMaxCut = TMath::Pi();
127
5ae09196 128 fMultiCutAna = kFALSE;
129
130 fNPtCuts = 3;
af7b3903 131 fPtCuts[0] = 0.; fPtCuts[1] = 0.3; fPtCuts[2] = 0.5;
132 for(Int_t i = fNPtCuts; i < 10; i++)fPtCuts[i] = 0.;
5ae09196 133
9c59b5fe 134 fNAsymCuts = 4;
135 fAsymCuts[0] = 1.; fAsymCuts[1] = 0.8; fAsymCuts[2] = 0.6; fAsymCuts[3] = 0.1;
af7b3903 136 for(Int_t i = fNAsymCuts; i < 10; i++)fAsymCuts[i] = 0.;
137
5ae09196 138 fNCellNCuts = 3;
af7b3903 139 fCellNCuts[0] = 0; fCellNCuts[1] = 1; fCellNCuts[2] = 2;
021c573a 140 for(Int_t i = fNCellNCuts; i < 10; i++)fCellNCuts[i] = 0;
af7b3903 141
142 fNPIDBits = 2;
143 fPIDBits[0] = 0; fPIDBits[1] = 2; // fPIDBits[2] = 4; fPIDBits[3] = 6;// check, no cut, dispersion, neutral, dispersion&&neutral
021c573a 144 for(Int_t i = fNPIDBits; i < 10; i++)fPIDBits[i] = 0;
af7b3903 145
1c5acb87 146}
1c5acb87 147
0c1383b5 148
149//________________________________________________________________________________________________________________________________________________
150TObjString * AliAnaPi0::GetAnalysisCuts()
151{
af7b3903 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") ;
157 parList+=onePar ;
158 snprintf(onePar,buffersize,"Number of bins in Centrality: %d \n",fNCentrBin) ;
159 parList+=onePar ;
160 snprintf(onePar,buffersize,"Number of bins in Z vert. pos: %d \n",GetNZvertBin()) ;
161 parList+=onePar ;
162 snprintf(onePar,buffersize,"Number of bins in Reac. Plain: %d \n",GetNRPBin()) ;
163 parList+=onePar ;
164 snprintf(onePar,buffersize,"Depth of event buffer: %d \n",fNmaxMixEv) ;
165 parList+=onePar ;
6175da48 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) ;
168 parList+=onePar ;
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) ;
af7b3903 170 parList+=onePar ;
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]);
173 parList+=onePar ;
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]);
176 parList+=onePar ;
177 snprintf(onePar,buffersize,"Cuts: \n") ;
178 parList+=onePar ;
179 snprintf(onePar,buffersize,"Z vertex position: -%f < z < %f \n",GetZvertexCut(),GetZvertexCut()) ;
180 parList+=onePar ;
181 snprintf(onePar,buffersize,"Calorimeter: %s \n",fCalorimeter.Data()) ;
182 parList+=onePar ;
183 snprintf(onePar,buffersize,"Number of modules: %d \n",fNModules) ;
184 parList+=onePar ;
db2bf6fd 185 if(fMultiCutAna){
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]);
188 parList+=onePar ;
db2bf6fd 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]);
191 parList+=onePar ;
db2bf6fd 192 }
193
af7b3903 194 return new TObjString(parList) ;
0c1383b5 195}
196
2e557d1c 197//________________________________________________________________________________________________________________________________________________
198TList * AliAnaPi0::GetCreateOutputObjects()
199{
477d6cee 200 // Create histograms to be saved in output file and
201 // store them in fOutputContainer
202
203 //create event containers
5025c139 204 fEventsList = new TList*[fNCentrBin*GetNZvertBin()*GetNRPBin()] ;
1c5acb87 205
477d6cee 206 for(Int_t ic=0; ic<fNCentrBin; ic++){
5025c139 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() ;
af7b3903 210 fEventsList[ic*GetNZvertBin()*GetNRPBin()+iz*GetNRPBin()+irp]->SetOwner(kFALSE);
477d6cee 211 }
212 }
213 }
7e7694bb 214
477d6cee 215 TList * outputContainer = new TList() ;
216 outputContainer->SetName(GetName());
6921fa00 217
af7b3903 218 fhReMod = new TH2D*[fNModules] ;
6175da48 219 fhReDiffMod = new TH2D*[fNModules+3] ;
220
221 fhMiMod = new TH2D*[fNModules] ;
222 fhMiDiffMod = new TH2D*[fNModules+3] ;
af7b3903 223
224 fhRe1 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
af7b3903 225 fhMi1 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
6175da48 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] ;
231 }
232 if(fMakeInvPtPlots) {
398c93cc 233 fhReInvPt1 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
398c93cc 234 fhMiInvPt1 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
6175da48 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] ;
240 }
398c93cc 241 }
6175da48 242
5ae09196 243 const Int_t buffersize = 255;
244 char key[buffersize] ;
245 char title[buffersize] ;
477d6cee 246
5a2dbc3c 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();
256
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();
af7b3903 263 Int_t ntrmbins = GetHistoTrackMultiplicityBins();
264 Int_t ntrmmax = GetHistoTrackMultiplicityMax();
265 Int_t ntrmmin = GetHistoTrackMultiplicityMin();
6175da48 266
267 fhAverTotECluster = new TH1F("hAverTotECluster","hAverTotECluster",200,0,50) ;
268 fhAverTotECluster->SetXTitle("E_{cluster, aver. SM} (GeV)");
269 outputContainer->Add(fhAverTotECluster) ;
270
271 fhAverTotECell = new TH1F("hAverTotECell","hAverTotECell",200,0,50) ;
272 fhAverTotECell->SetXTitle("E_{cell, aver. SM} (GeV)");
273 outputContainer->Add(fhAverTotECell) ;
274
156549ae 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) ;
279
280 fhEDensityCluster = new TH1F("hEDensityCluster","hEDensityCluster",200,0,50) ;
281 fhEDensityCluster->SetXTitle("#Sigma E_{cluster} / N_{cluster} (GeV)");
282 outputContainer->Add(fhEDensityCluster) ;
283
284 fhEDensityCell = new TH1F("hEDensityCell","hEDensityCell",200,0,50) ;
285 fhEDensityCell->SetXTitle("#Sigma E_{cell} / N_{cell} (GeV)");
286 outputContainer->Add(fhEDensityCell) ;
287
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) ;
292
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) ;
296//
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) ;
300//
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) ;
304//
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) ;
308//
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) ;
313//
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) ;
318
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) ;
323//
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) ;
328
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) ;
333//
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) ;
338
6175da48 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) ;
343
6175da48 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) ;
348
156549ae 349 if(fDoOwnMix){
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) ;
354
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) ;
359 }
6175da48 360
477d6cee 361 for(Int_t ic=0; ic<fNCentrBin; ic++){
6175da48 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);
398c93cc 366 //Distance to bad module 1
6175da48 367 snprintf(key, buffersize,"hRe_cen%d_pidbit%d_asy%d_dist1",ic,ipid,iasym) ;
398c93cc 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]) ;
6175da48 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]) ;
af7b3903 375
6175da48 376 if(fFillBadDistHisto){
398c93cc 377 //Distance to bad module 2
6175da48 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",
398c93cc 380 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
6175da48 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]) ;
398c93cc 385
386 //Distance to bad module 3
6175da48 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",
398c93cc 389 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
6175da48 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]) ;
398c93cc 394 }
6175da48 395
396 //Inverse pT
397 if(fMakeInvPtPlots){
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]) ;
406
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]) ;
416
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]) ;
425 }
426 }
427 if(fDoOwnMix){
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]) ;
445
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]) ;
454 }
455 //Inverse pT
456 if(fMakeInvPtPlots){
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]) ;
474
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]) ;
483 }
484 }
485 }
486 }
7e7694bb 487 }
1c5acb87 488 }
477d6cee 489
9c59b5fe 490 fhRePtAsym = new TH2D("hRePtAsym","Asymmetry vs pt, for pairs",nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax) ;
af7b3903 491 fhRePtAsym->SetXTitle("p_{T} (GeV/c)");
492 fhRePtAsym->SetYTitle("Asymmetry");
493 outputContainer->Add(fhRePtAsym);
494
9c59b5fe 495 fhRePtAsymPi0 = new TH2D("hRePtAsymPi0","Asymmetry vs pt, for pairs close to #pi^{0} mass",nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax) ;
af7b3903 496 fhRePtAsymPi0->SetXTitle("p_{T} (GeV/c)");
497 fhRePtAsymPi0->SetYTitle("Asymmetry");
498 outputContainer->Add(fhRePtAsymPi0);
499
9c59b5fe 500 fhRePtAsymEta = new TH2D("hRePtAsymEta","Asymmetry vs pt, for pairs close to #eta mass",nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax) ;
af7b3903 501 fhRePtAsymEta->SetXTitle("p_{T} (GeV/c)");
502 fhRePtAsymEta->SetYTitle("Asymmetry");
503 outputContainer->Add(fhRePtAsymEta);
504
5ae09196 505 if(fMultiCutAna){
506
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) ;
af7b3903 512 fhRePIDBits[ipid]->SetXTitle("p_{T} (GeV/c)");
513 fhRePIDBits[ipid]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
5ae09196 514 outputContainer->Add(fhRePIDBits[ipid]) ;
515 }// pid bit loop
516
6175da48 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];
5ae09196 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) ;
af7b3903 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]) ;
5ae09196 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) ;
af7b3903 531 fhRePtNCellAsymCuts[index]->SetXTitle("p_{T} (GeV/c)");
532 fhRePtNCellAsymCuts[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
5ae09196 533 outputContainer->Add(fhRePtNCellAsymCuts[index]) ;
6175da48 534
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]) ;
541
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]) ;
548
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]) ;
555
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]) ;
562
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]) ;
569
5ae09196 570 }
571 }
572 }
821c8090 573
af7b3903 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]) ;
582 }
583
5ae09196 584 }// multi cuts analysis
585
477d6cee 586 fhEvents=new TH3D("hEvents","Number of events",fNCentrBin,0.,1.*fNCentrBin,
5025c139 587 GetNZvertBin(),0.,1.*GetNZvertBin(),GetNRPBin(),0.,1.*GetNRPBin()) ;
477d6cee 588 outputContainer->Add(fhEvents) ;
50f39b97 589
590 fhRealOpeningAngle = new TH2D
6175da48 591 ("hRealOpeningAngle","Angle between all #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,300,0,TMath::Pi());
50f39b97 592 fhRealOpeningAngle->SetYTitle("#theta(rad)");
593 fhRealOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
594 outputContainer->Add(fhRealOpeningAngle) ;
7e7694bb 595
50f39b97 596 fhRealCosOpeningAngle = new TH2D
6175da48 597 ("hRealCosOpeningAngle","Cosinus of angle between all #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,100,0,1);
50f39b97 598 fhRealCosOpeningAngle->SetYTitle("cos (#theta) ");
599 fhRealCosOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
600 outputContainer->Add(fhRealCosOpeningAngle) ;
601
6175da48 602 if(fDoOwnMix){
603
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) ;
609
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) ;
615
616 }
617
477d6cee 618 //Histograms filled only if MC data is requested
0ae57829 619 if(IsDataMC()){
6175da48 620 //Pi0
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) ;
156549ae 623 fhPrimPi0Pt ->SetXTitle("p_{T} (GeV/c)");
624 fhPrimPi0AccPt->SetXTitle("p_{T} (GeV/c)");
6175da48 625 outputContainer->Add(fhPrimPi0Pt) ;
626 outputContainer->Add(fhPrimPi0AccPt) ;
627
156549ae 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)");
6175da48 631 outputContainer->Add(fhPrimPi0Y) ;
632
156549ae 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)");
6175da48 636 outputContainer->Add(fhPrimPi0AccY) ;
477d6cee 637
156549ae 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)");
6175da48 641 outputContainer->Add(fhPrimPi0Phi) ;
477d6cee 642
156549ae 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)");
6175da48 646 outputContainer->Add(fhPrimPi0AccPhi) ;
477d6cee 647
6175da48 648 //Eta
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) ;
156549ae 651 fhPrimEtaPt ->SetXTitle("p_{T} (GeV/c)");
652 fhPrimEtaAccPt->SetXTitle("p_{T} (GeV/c)");
6175da48 653 outputContainer->Add(fhPrimEtaPt) ;
654 outputContainer->Add(fhPrimEtaAccPt) ;
477d6cee 655
156549ae 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)");
6175da48 659 outputContainer->Add(fhPrimEtaY) ;
50f39b97 660
156549ae 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)");
6175da48 664 outputContainer->Add(fhPrimEtaAccY) ;
50f39b97 665
156549ae 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)");
6175da48 669 outputContainer->Add(fhPrimEtaPhi) ;
670
156549ae 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)");
6175da48 674 outputContainer->Add(fhPrimEtaAccPhi) ;
675
50f39b97 676
6175da48 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) ;
682
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) ;
688
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]) ;
694
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]) ;
699
156549ae 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) ;
6175da48 701 fhMCOrgDeltaEta[i]->SetXTitle("p_{T} (GeV/c)");
156549ae 702 fhMCOrgDeltaEta[i]->SetYTitle("#Delta #eta");
6175da48 703 outputContainer->Add(fhMCOrgDeltaEta[i]) ;
704
156549ae 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) ;
6175da48 706 fhMCOrgDeltaPhi[i]->SetXTitle("p_{T} (GeV/c)");
156549ae 707 fhMCOrgDeltaPhi[i]->SetYTitle("#Delta #phi (rad)");
6175da48 708 outputContainer->Add(fhMCOrgDeltaPhi[i]) ;
709
710 }
50f39b97 711
6175da48 712 if(fMultiCutAnaSim){
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;
723
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]) ;
730
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]) ;
737
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]) ;
744
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]) ;
751
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]) ;
758
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]) ;
765 }
766 }
767 }
768 }//multi cut ana
769 else {
770 fhMCPi0MassPtTrue = new TH2D*[1];
771 fhMCPi0PtTruePtRec = new TH2D*[1];
772 fhMCEtaMassPtTrue = new TH2D*[1];
773 fhMCEtaPtTruePtRec = new TH2D*[1];
774
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]) ;
779
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]) ;
784
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]) ;
789
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]) ;
794 }
477d6cee 795 }
50f39b97 796
6175da48 797 TString * pairname = new TString[fNModules+3];
821c8090 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)";
6175da48 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]="";}
821c8090 807 if(fCalorimeter=="PHOS") {
808 pairname[0]="(0-1)";
809 pairname[1]="(0-2)";
af7b3903 810 pairname[2]="(1-2)";
821c8090 811 for(Int_t i = 3 ; i < fNModules ; i++) pairname[i]="";}
812
6921fa00 813 for(Int_t imod=0; imod<fNModules; imod++){
50f39b97 814 //Module dependent invariant mass
5ae09196 815 snprintf(key, buffersize,"hReMod_%d",imod) ;
816 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for Module %d",imod) ;
af7b3903 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})");
50f39b97 820 outputContainer->Add(fhReMod[imod]) ;
821c8090 821
822 snprintf(key, buffersize,"hReDiffMod_%d",imod) ;
823 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for Different Modules: %s",(pairname[imod]).Data()) ;
af7b3903 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})");
821c8090 827 outputContainer->Add(fhReDiffMod[imod]) ;
6175da48 828
829 if(fDoOwnMix){
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]) ;
836
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]) ;
843 }
844
6921fa00 845 }
50f39b97 846
6175da48 847 for (Int_t imod=4; imod<7; imod++) {
848
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]) ;
821c8090 855
6175da48 856 if(fDoOwnMix){
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]) ;
863 }
864 }
821c8090 865
6175da48 866 delete [] pairname;
821c8090 867
eee5fcf1 868// for(Int_t i = 0; i < outputContainer->GetEntries() ; i++){
869//
870// printf("Histogram %d, name: %s\n ",i, outputContainer->At(i)->GetName());
871//
872// }
873
477d6cee 874 return outputContainer;
1c5acb87 875}
876
877//_________________________________________________________________________________________________________________________________________________
878void AliAnaPi0::Print(const Option_t * /*opt*/) const
879{
477d6cee 880 //Print some relevant parameters set for the analysis
a3aebfff 881 printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
477d6cee 882 AliAnaPartCorrBaseClass::Print(" ");
a3aebfff 883
477d6cee 884 printf("Number of bins in Centrality: %d \n",fNCentrBin) ;
5025c139 885 printf("Number of bins in Z vert. pos: %d \n",GetNZvertBin()) ;
886 printf("Number of bins in Reac. Plain: %d \n",GetNRPBin()) ;
477d6cee 887 printf("Depth of event buffer: %d \n",fNmaxMixEv) ;
af7b3903 888 printf("Pair in same Module: %d \n",fSameSM) ;
477d6cee 889 printf("Cuts: \n") ;
5025c139 890 printf("Z vertex position: -%2.3f < z < %2.3f \n",GetZvertexCut(),GetZvertexCut()) ;
50f39b97 891 printf("Number of modules: %d \n",fNModules) ;
6175da48 892 printf("Select pairs with their angle: %d, edep %d, min angle %2.3f, max angle %2.3f \n",fUseAngleCut, fUseAngleEDepCut, fAngleCut, fAngleMaxCut) ;
af7b3903 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]);
896 printf("\n");
897
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]);
901 printf("\n");
902
db2bf6fd 903 if(fMultiCutAna){
904 printf("pT cuts: n = %d, \n",fNPtCuts) ;
905 printf("\tpT > ");
906 for(Int_t i = 0; i < fNPtCuts; i++) printf("%2.2f ",fPtCuts[i]);
907 printf("GeV/c\n");
908
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]);
912 printf("\n");
913
db2bf6fd 914 }
477d6cee 915 printf("------------------------------------------------------\n") ;
1c5acb87 916}
917
5ae09196 918//_____________________________________________________________
919void AliAnaPi0::FillAcceptanceHistograms(){
920 //Fill acceptance histograms if MC data is available
c8fe2783 921
6175da48 922 if(GetReader()->ReadStack()){
5ae09196 923 AliStack * stack = GetMCStack();
6175da48 924 if(stack){
5ae09196 925 for(Int_t i=0 ; i<stack->GetNprimary(); i++){
926 TParticle * prim = stack->Particle(i) ;
6175da48 927 Int_t pdg = prim->GetPdgCode();
928 if( pdg == 111 || pdg == 221){
5ae09196 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() ;
6175da48 934 if(pdg == 111){
935 if(TMath::Abs(pi0Y) < 0.5){
936 fhPrimPi0Pt->Fill(pi0Pt) ;
937 }
156549ae 938 fhPrimPi0Y ->Fill(pi0Pt, pi0Y) ;
939 fhPrimPi0Phi->Fill(pi0Pt, phi) ;
6175da48 940 }
941 else if(pdg == 221){
942 if(TMath::Abs(pi0Y) < 0.5){
943 fhPrimEtaPt->Fill(pi0Pt) ;
944 }
156549ae 945 fhPrimEtaY ->Fill(pi0Pt, pi0Y) ;
946 fhPrimEtaPhi->Fill(pi0Pt, phi) ;
5ae09196 947 }
5ae09196 948 //Check if both photons hit Calorimeter
6175da48 949 if(prim->GetNDaughters()!=2) return; //Only interested in 2 gamma decay
5ae09196 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());
958
959 TLorentzVector lv1, lv2;
960 phot1->Momentum(lv1);
961 phot2->Momentum(lv2);
962
963 Bool_t inacceptance = kFALSE;
964 if(fCalorimeter == "PHOS"){
965 if(GetPHOSGeometry() && GetCaloUtils()->IsPHOSGeoMatrixSet()){
966 Int_t mod ;
967 Double_t x,z ;
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);
971 }
972 else{
973
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);
977 }
978
979 }
980 else if(fCalorimeter == "EMCAL" && GetCaloUtils()->IsEMCALGeoMatrixSet()){
981 if(GetEMCALGeometry()){
156549ae 982
983 Int_t absID1=0;
984 Int_t absID2=0;
985
986 GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(phot1->Eta(),phot1->Phi(),absID1);
987 GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(phot2->Eta(),phot2->Phi(),absID2);
988
989 if( absID1 >= 0 && absID2 >= 0)
5ae09196 990 inacceptance = kTRUE;
156549ae 991
992// if(GetEMCALGeometry()->Impact(phot1) && GetEMCALGeometry()->Impact(phot2))
993// inacceptance = kTRUE;
5ae09196 994 if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
995 }
996 else{
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);
1000 }
1001 }
1002
1003 if(inacceptance){
6175da48 1004 if(pdg==111){
156549ae 1005 fhPrimPi0AccPt ->Fill(pi0Pt) ;
1006 fhPrimPi0AccPhi->Fill(pi0Pt, phi) ;
1007 fhPrimPi0AccY ->Fill(pi0Pt, pi0Y) ;
6175da48 1008 Double_t angle = lv1.Angle(lv2.Vect());
1009 fhPrimPi0OpeningAngle ->Fill(pi0Pt,angle);
1010 fhPrimPi0CosOpeningAngle->Fill(pi0Pt,TMath::Cos(angle));
1011 }
1012 else if(pdg==221){
156549ae 1013 fhPrimEtaAccPt ->Fill(pi0Pt) ;
1014 fhPrimEtaAccPhi->Fill(pi0Pt, phi) ;
1015 fhPrimEtaAccY ->Fill(pi0Pt, pi0Y) ;
6175da48 1016 }
5ae09196 1017 }//Accepted
1018 }// 2 photons
1019 }//Check daughters exist
156549ae 1020 }// Primary pi0 or eta
5ae09196 1021 }//loop on primaries
1022 }//stack exists and data is MC
1023 }//read stack
1024 else if(GetReader()->ReadAODMCParticles()){
6175da48 1025
1026 TClonesArray * mcparticles = GetReader()->GetAODMCParticles(0);
1027 if(mcparticles){
1028 Int_t nprim = mcparticles->GetEntriesFast();
156549ae 1029 for(Int_t i=0
1030 ; i < nprim; i++)
1031 {
6175da48 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() ;
1040 if(pdg == 111){
1041 if(TMath::Abs(pi0Y) < 0.5){
1042 fhPrimPi0Pt->Fill(pi0Pt) ;
1043 }
156549ae 1044 fhPrimPi0Y ->Fill(pi0Pt, pi0Y) ;
1045 fhPrimPi0Phi->Fill(pi0Pt, phi) ;
6175da48 1046 }
1047 else if(pdg == 221){
1048 if(TMath::Abs(pi0Y) < 0.5){
1049 fhPrimEtaPt->Fill(pi0Pt) ;
1050 }
156549ae 1051 fhPrimEtaY ->Fill(pi0Pt, pi0Y) ;
1052 fhPrimEtaPhi->Fill(pi0Pt, phi) ;
6175da48 1053 }
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){
6175da48 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());
1065
1066 Bool_t inacceptance = kFALSE;
1067 if(fCalorimeter == "PHOS"){
1068 if(GetPHOSGeometry() && GetCaloUtils()->IsPHOSGeoMatrixSet()){
1069 Int_t mod ;
1070 Double_t x,z ;
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);
1077 }
1078 else{
1079
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);
1083 }
1084
1085 }
1086 else if(fCalorimeter == "EMCAL" && GetCaloUtils()->IsEMCALGeoMatrixSet()){
1087 if(GetEMCALGeometry()){
156549ae 1088
6175da48 1089 Int_t absID1=0;
6175da48 1090 Int_t absID2=0;
156549ae 1091
1092 //TVector3 vtx(phot1->Xv(),phot1->Yv(),phot1->Zv());
1093 //TVector3 vimpact(0,0,0);
1094
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);
1099
1100 GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(phot1->Eta(),phot1->Phi(),absID1);
1101 GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(phot2->Eta(),phot2->Phi(),absID2);
1102
6175da48 1103// if(TMath::Abs(phot1->Eta()) < 0.7 && phot1->Phi() > 80*TMath::DegToRad() && phot1->Phi() < 120*TMath::DegToRad() )
156549ae 1104// printf("Phot1 ccepted? %d\n",absID1);
6175da48 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);
1107
1108 if( absID1 >= 0 && absID2 >= 0)
1109 inacceptance = kTRUE;
156549ae 1110
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());
1114
1115
1116
6175da48 1117 if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
1118 }
1119 else{
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);
1123 }
1124 }
1125
1126 if(inacceptance){
1127 if(pdg==111){
156549ae 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) ;
6175da48 1132 Double_t angle = lv1.Angle(lv2.Vect());
1133 fhPrimPi0OpeningAngle ->Fill(pi0Pt,angle);
1134 fhPrimPi0CosOpeningAngle->Fill(pi0Pt,TMath::Cos(angle));
1135 }
1136 else if(pdg==221){
156549ae 1137 fhPrimEtaAccPt ->Fill(pi0Pt) ;
1138 fhPrimEtaAccPhi->Fill(pi0Pt, phi) ;
1139 fhPrimEtaAccY ->Fill(pi0Pt, pi0Y) ;
6175da48 1140 }
1141 }//Accepted
1142 }// 2 photons
1143 }//Check daughters exist
156549ae 1144 }// Primary pi0 or eta
6175da48 1145 }//loop on primaries
1146 }//stack exists and data is MC
1147
1148
1149 } // read AOD MC
5ae09196 1150}
1151
6175da48 1152//_____________________________________________________________
1153void 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
1162
1163 Int_t ancPDG = 0;
1164 Int_t ancStatus = 0;
1165 TLorentzVector ancMomentum;
1166 Int_t ancLabel = GetMCAnalysisUtils()->CheckCommonAncestor(index1, index2,
1167 GetReader(), ancPDG, ancStatus,ancMomentum);
1168
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);
1171
1172 if(ancLabel > -1){
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);
1178 }
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);
1184 }
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
1203 }// icell loop
1204 }// pt cut loop
1205 }//Multi cut ana sim
1206 else {
1207 fhMCPi0MassPtTrue[0]->Fill(ancMomentum.Pt(),mass);
1208 if(mass < 0.17 && mass > 0.1) fhMCPi0PtTruePtRec[0]->Fill(ancMomentum.Pt(),pt);
1209 }
1210 }
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);
156549ae 1226 if(mass < 0.65 && mass > 0.45) fhMCEtaPtTruePtRec[index]->Fill(ancMomentum.Pt(),pt);
6175da48 1227 }//pass the different cuts
1228 }// pid bit cut loop
1229 }// icell loop
1230 }// pt cut loop
1231 } //Multi cut ana sim
1232 else {
1233 fhMCEtaMassPtTrue[0]->Fill(ancMomentum.Pt(),mass);
156549ae 1234 if(mass < 0.65 && mass > 0.45) fhMCEtaPtTruePtRec[0]->Fill(ancMomentum.Pt(),pt);
6175da48 1235 }
1236 }
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);
1242 }
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);
1248 }
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);
1254 }
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);
1261 }
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);
1267 }
1268 }
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);
1275 }
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);
1288 }
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);
1294 }
1295 else {
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);
1298 }
1299 }//status 21
1300 else {
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);
1303 }
1304 }////Partons, colliding protons, strings, intermediate corrections
1305 }//ancLabel > -1
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);
1312 }
1313}
1314
1c5acb87 1315//____________________________________________________________________________________________________________________________________________________
6639984f 1316void AliAnaPi0::MakeAnalysisFillHistograms()
1c5acb87 1317{
477d6cee 1318 //Process one event and extract photons from AOD branch
1319 // filled with AliAnaPhoton and fill histos with invariant mass
1320
6175da48 1321 //In case of simulated data, fill acceptance histograms
1322 if(IsDataMC())FillAcceptanceHistograms();
156549ae 1323 if (GetReader()->GetEventNumber()%10000 == 0)
1324 printf("--- Event %d ---\n",GetReader()->GetEventNumber());
6175da48 1325 //Init some variables
1326//Int_t iRun = (GetReader()->GetInputEvent())->GetRunNumber() ;
1327 Int_t nPhot = GetInputAODBranch()->GetEntriesFast() ;
1328 Int_t nClus = 0;
1329 Int_t nCell = 0;
1330 Float_t eClusTot = 0;
1331 Float_t eCellTot = 0;
156549ae 1332 Float_t eDenClus = 0;
1333 Float_t eDenCell = 0;
1334// Int_t ncomb = 0;
1335// Float_t rtmp = 0;
1336// Float_t rtmpw = 0;
1337// Float_t rxz = 0;
1338// Float_t rxzw = 0;
1339// Float_t pos1[3];
1340// Float_t pos2[3];
1341// Float_t emax = 0;
477d6cee 1342
156549ae 1343 if(GetDebug() > 1)
1344 printf("AliAnaPi0::MakeAnalysisFillHistograms() - Photon entries %d\n", nPhot);
1345
1346 //If less than photon 2 entries in the list, skip this event
1347 if(nPhot < 2 ) return ;
1348
6175da48 1349 // Count the number of clusters and cells, in case multiplicity bins dependent on such numbers
1350 // are requested
1351 if(fCalorimeter=="EMCAL"){
1352 nClus = GetAODEMCAL() ->GetEntriesFast();
1353 nCell = GetEMCALCells()->GetNumberOfCells();
156549ae 1354 for(Int_t icl=0; icl < nClus; icl++) {
1355 Float_t e1 = ((AliVCluster*)GetAODEMCAL()->At(icl))->E();
1356 eClusTot += e1;
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);
1364// rxz += rtmp;
1365// rxzw += rtmpw;
1366// ncomb++;
1367// fhClusterPairDist ->Fill(rtmp);
1368// fhClusterPairDistWeight->Fill(rtmpw);
1369// //printf("Distance: %f; weighted %f\n ",rtmp,rtmp/(e1+((AliVCluster*)GetAODEMCAL()->At(icl2))->E()));
1370//
1371// }// second cluster loop
1372 }// first cluster
1373
6175da48 1374 for(Int_t jce=0; jce < nCell; jce++) eCellTot += GetEMCALCells()->GetAmplitude(jce);
6175da48 1375 }
1376 else {
1377 nClus = GetAODPHOS() ->GetEntriesFast();
1378 nCell = GetPHOSCells()->GetNumberOfCells();
156549ae 1379 for(Int_t icl=0; icl < nClus; icl++) {
1380 Float_t e1 = ((AliVCluster*)GetAODPHOS()->At(icl))->E();
1381 eClusTot += e1;
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);
1388// rxz += rtmp;
1389// rxzw += rtmpw;
1390// ncomb++;
1391// fhClusterPairDist ->Fill(rtmp);
1392// fhClusterPairDistWeight->Fill(rtmpw);
1393// }// second cluster loop
1394 }// first cluster
6175da48 1395 for(Int_t jce=0; jce < nCell; jce++) eCellTot += GetPHOSCells()->GetAmplitude(jce);
1396 }
156549ae 1397 if(GetDebug() > 1)
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);
6175da48 1399
156549ae 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);
6175da48 1406 //Fill the average number of cells or clusters per SM
1407 eClusTot /=fNModules;
1408 eCellTot /=fNModules;
156549ae 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);
6175da48 1413
156549ae 1414// //Average weighted pair distance
1415// rxz /= ncomb;
1416// rxzw /= ncomb;
1417//
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);
1424//
1425// //emax
1426// fhMaxEvsClustEDen->Fill(emax,eDenClus);
1427// fhMaxEvsClustMult->Fill(emax,nPhot);
7e7694bb 1428
156549ae 1429 //printf("Average Distance: %f; weighted %f\n ",rxz,rxzw);
1430
6175da48 1431
1432 //Init variables
1433 Int_t module1 = -1;
1434 Int_t module2 = -1;
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 ;
1441
1442 //---------------------------------
1443 //First loop on photons/clusters
1444 //---------------------------------
477d6cee 1445 for(Int_t i1=0; i1<nPhot-1; i1++){
1446 AliAODPWG4Particle * p1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i1)) ;
6175da48 1447 //printf("AliAnaPi0::MakeAnalysisFillHistograms() : cluster1 id %d\n",p1->GetCaloLabel(0));
1448
7e7694bb 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
c8fe2783 1451 evtIndex1 = GetEventIndex(p1, vert) ;
5025c139 1452 //printf("charge = %d\n", track->Charge());
c8fe2783 1453 if ( evtIndex1 == -1 )
1454 return ;
1455 if ( evtIndex1 == -2 )
1456 continue ;
2244659d 1457 if(TMath::Abs(vert[2]) > GetZvertexCut()) continue ; //vertex cut
6175da48 1458
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 //----------------------------------------------------------------------------
c8fe2783 1464 if (evtIndex1 != currentEvtIndex) {
6175da48 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);
1470 }
1471 else if(fUsePhotonMultBins){ // Photon multiplicity bins
1472 //printf("photon mult %d cluster mult %d\n",nPhot, nClus);
156549ae 1473 curRPBin = nPhot-2;
1474 if(curRPBin > GetNRPBin() -1) curRPBin=GetNRPBin()-1;
1475 //printf("photon mult bin %d\n",curRPBin);
6175da48 1476 }
156549ae 1477 else if(fUseAverClusterEBins){ // Cluster average energy bins
6175da48 1478 //Bins for pp, if needed can be done in a more general way
156549ae 1479 curCentrBin = eClusTot/10 * fNCentrBin;
6175da48 1480 if(curCentrBin > fNCentrBin-1) curCentrBin=fNCentrBin-1;
1481 //printf("cluster E average %f, bin %d \n",eClusTot,curCentrBin);
1482 }
156549ae 1483 else if(fUseAverCellEBins){ // Cell average energy bins
6175da48 1484 //Bins for pp, if needed can be done in a more general way
156549ae 1485 curCentrBin = eCellTot/10*fNCentrBin;
6175da48 1486 if(curCentrBin > fNCentrBin-1) curCentrBin=fNCentrBin-1;
1487 //printf("cell E average %f, bin %d \n",eCellTot,curCentrBin);
1488 }
156549ae 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);
1494 }
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);
1500// }
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);
1506// }
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);
1512// }
6175da48 1513 else { //Event centrality
1514 curCentrBin = GetEventCentrality();
1515 }
1516
156549ae 1517 //Reaction plane bin
c8fe2783 1518 curRPBin = 0 ;
156549ae 1519
1520 //Get vertex z bin
5025c139 1521 curZvertBin = (Int_t)(0.5*GetNZvertBin()*(vert[2]+GetZvertexCut())/GetZvertexCut()) ;
6175da48 1522
1523 //Fill event bin info
c8fe2783 1524 fhEvents->Fill(curCentrBin+0.5,curZvertBin+0.5,curRPBin+0.5) ;
1525 currentEvtIndex = evtIndex1 ;
ca468d44 1526 if(GetDebug() > 1)
6175da48 1527 printf("AliAnaPi0::MakeAnalysisFillHistograms() - Centrality %d, Vertex Bin %d, RP bin %d \n",curCentrBin,curRPBin,curZvertBin);
c8fe2783 1528 }
7e7694bb 1529
f8006433 1530 //printf("AliAnaPi0::MakeAnalysisFillHistograms(): Photon 1 Evt %d Vertex : %f,%f,%f\n",evtIndex1, GetVertex(evtIndex1)[0] ,GetVertex(evtIndex1)[1],GetVertex(evtIndex1)[2]);
af7b3903 1531
6175da48 1532 //Get the momentum of this cluster
477d6cee 1533 TLorentzVector photon1(p1->Px(),p1->Py(),p1->Pz(),p1->E());
6175da48 1534
1535 //Get (Super)Module number of this cluster
59b6bd99 1536 module1 = GetModuleNumber(p1);
6175da48 1537
1538 //---------------------------------
1539 //Second loop on photons/clusters
1540 //---------------------------------
477d6cee 1541 for(Int_t i2=i1+1; i2<nPhot; i2++){
1542 AliAODPWG4Particle * p2 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i2)) ;
6175da48 1543
1544 //In case of mixing frame, check we are not in the same event as the first cluster
c8fe2783 1545 Int_t evtIndex2 = GetEventIndex(p2, vert) ;
1546 if ( evtIndex2 == -1 )
1547 return ;
1548 if ( evtIndex2 == -2 )
1549 continue ;
1550 if (GetMixedEvent() && (evtIndex1 == evtIndex2))
7e7694bb 1551 continue ;
6175da48 1552
f8006433 1553 //printf("AliAnaPi0::MakeAnalysisFillHistograms(): Photon 2 Evt %d Vertex : %f,%f,%f\n",evtIndex2, GetVertex(evtIndex2)[0] ,GetVertex(evtIndex2)[1],GetVertex(evtIndex2)[2]);
6175da48 1554
1555 //Get the momentum of this cluster
477d6cee 1556 TLorentzVector photon2(p2->Px(),p2->Py(),p2->Pz(),p2->E());
59b6bd99 1557 //Get module number
6175da48 1558 module2 = GetModuleNumber(p2);
1559
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()) ;
1568
477d6cee 1569 if(GetDebug() > 2)
6175da48 1570 printf(" E: photon1 %f, photon2 %f; Pair: pT %f, mass %f, a %f\n", p1->E(), p2->E(), (photon1 + photon2).E(),m,a);
1571
1572 //--------------------------------
1573 // Opening angle selection
1574 //--------------------------------
50f39b97 1575 //Check if opening angle is too large or too small compared to what is expected
1576 Double_t angle = photon1.Angle(photon2.Vect());
6175da48 1577 if(fUseAngleEDepCut && !GetNeutralMesonSelection()->IsAngleInWindow((photon1+photon2).E(),angle+0.05)) {
1578 if(GetDebug() > 2)
1579 printf("AliAnaPi0::MakeAnalysisFillHistograms() -Real pair angle %f not in E %f window\n",angle, (photon1+photon2).E());
c8fe2783 1580 continue;
6175da48 1581 }
af7b3903 1582
6175da48 1583 if(fUseAngleCut && (angle < fAngleCut || angle > fAngleMaxCut)) {
1584 if(GetDebug() > 2)
1585 printf("AliAnaPi0::MakeAnalysisFillHistograms() - Real pair cut %f < angle %f < cut %f\n",fAngleCut, angle, fAngleMaxCut);
1586 continue;
1587 }
1588
1589 //-------------------------------------------------------------------------------------------------
af7b3903 1590 //Fill module dependent histograms, put a cut on assymmetry on the first available cut in the array
6175da48 1591 //-------------------------------------------------------------------------------------------------
af7b3903 1592 if(a < fAsymCuts[0]){
1593 if(module1==module2 && module1 >=0 && module1<fNModules)
1594 fhReMod[module1]->Fill(pt,m) ;
1595 else
6175da48 1596 fhReDiffMod[fNModules+2]->Fill(pt,m) ;
af7b3903 1597
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) ;
6175da48 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) ;
af7b3903 1605 }
1606 else {
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) ;
1610 }
821c8090 1611 }
7e7694bb 1612
af7b3903 1613 //In case we want only pairs in same (super) module, check their origin.
1614 Bool_t ok = kTRUE;
1615 if(fSameSM && module1!=module2) ok=kFALSE;
1616 if(ok){
6175da48 1617
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);
1621
af7b3903 1622 //Fill histograms for different bad channel distance, centrality, assymmetry cut and pid bit
5ae09196 1623 for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
af7b3903 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;
6175da48 1628 //printf("index %d :(cen %d * nPID %d + ipid %d)*nasym %d + iasym %d\n",index,curCentrBin,fNPIDBits,ipid,fNAsymCuts,iasym);
af7b3903 1629 fhRe1 [index]->Fill(pt,m);
398c93cc 1630 if(fMakeInvPtPlots)fhReInvPt1[index]->Fill(pt,m,1./pt) ;
6175da48 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) ;
1638 }// bad 3
1639 }// bad2
1640 }// Fill bad dist histos
1641 }//assymetry cut
1642 }// asymmetry cut loop
af7b3903 1643 }// bad 1
1644 }// pid bit loop
5ae09196 1645
af7b3903 1646 //Fill histograms with opening angle
1647 fhRealOpeningAngle ->Fill(pt,angle);
1648 fhRealCosOpeningAngle->Fill(pt,TMath::Cos(angle));
1649
1650 //Fill histograms with pair assymmetry
1651 fhRePtAsym->Fill(pt,a);
6175da48 1652 if(m > 0.10 && m < 0.17) fhRePtAsymPi0->Fill(pt,a);
af7b3903 1653 if(m > 0.45 && m < 0.65) fhRePtAsymEta->Fill(pt,a);
1654
6175da48 1655 //-------------------------------------------------------
1656 //Get the number of cells needed for multi cut analysis.
1657 //-------------------------------------------------------
1658 Int_t ncell1 = 0;
1659 Int_t ncell2 = 0;
1660 if(fMultiCutAna || (IsDataMC() && fMultiCutAnaSim)){
1661
af7b3903 1662 AliVEvent * event = GetReader()->GetInputEvent();
1663 if(event){
1664 for(Int_t iclus = 0; iclus < event->GetNumberOfCaloClusters(); iclus++){
1665 AliVCluster *cluster = event->GetCaloCluster(iclus);
5ae09196 1666
af7b3903 1667 Bool_t is = kFALSE;
1668 if (fCalorimeter == "EMCAL" && GetReader()->IsEMCALCluster(cluster)) is = kTRUE;
1669 else if(fCalorimeter == "PHOS" && GetReader()->IsPHOSCluster (cluster)) is = kTRUE;
5ae09196 1670
af7b3903 1671 if(is){
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
1675
1676 if(ncell2 > 0 && ncell1 > 0) break; // No need to continue the iteration
1677
1678 }
1679 //printf("e 1: %2.2f, e 2: %2.2f, ncells: n1 %d, n2 %d\n", p1->E(), p2->E(),ncell1,ncell2);
1680 }
6175da48 1681 }
1682
1683 //---------
1684 // MC data
1685 //---------
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);
1688
1689 //-----------------------
1690 //Multi cuts analysis
1691 //-----------------------
1692 if(fMultiCutAna){
1693 //Histograms for different PID bits selection
1694 for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
1695
1696 if(p1->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton) &&
1697 p2->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton)) fhRePIDBits[ipid]->Fill(pt,m) ;
1698
1699 //printf("ipt %d, ipid%d, name %s\n",ipt, ipid, fhRePtPIDCuts[ipt*fNPIDBitsBits+ipid]->GetName());
1700 } // pid bit cut loop
1701
1702 //Several pt,ncell and asymmetry cuts
af7b3903 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] &&
6175da48 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");
1718 }
1719 }
af7b3903 1720 }// pid bit cut loop
1721 }// icell loop
1722 }// pt cut loop
1723 for(Int_t iasym = 0; iasym < fNAsymCuts; iasym++){
1724 if(a < fAsymCuts[iasym])fhRePtMult[iasym]->Fill(pt,GetTrackMultiplicity(),m) ;
1725 }
af7b3903 1726 }// multiple cuts analysis
1727 }// ok if same sm
7e7694bb 1728 }// second same event particle
1729 }// first cluster
6175da48 1730
1731 //-------------------------------------------------------------
1732 // Mixing
1733 //-------------------------------------------------------------
7e7694bb 1734 if(fDoOwnMix){
156549ae 1735 //printf("Cen bin %d, RP bin %d, e aver %f, mult %d\n",curCentrBin,curRPBin, eClusTot, nPhot);
6175da48 1736 //Recover events in with same characteristics as the current event
5025c139 1737 TList * evMixList=fEventsList[curCentrBin*GetNZvertBin()*GetNRPBin()+curZvertBin*GetNRPBin()+curRPBin] ;
7e7694bb 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() ;
1742 Double_t m = -999;
6175da48 1743 if(GetDebug() > 1)
1744 printf("AliAnaPi0::MakeAnalysisFillHistograms() - Mixed event %d photon entries %d, centrality bin %d\n", ii, nPhot2, curCentrBin);
7e7694bb 1745
6175da48 1746 //---------------------------------
1747 //First loop on photons/clusters
1748 //---------------------------------
7e7694bb 1749 for(Int_t i1=0; i1<nPhot; i1++){
1750 AliAODPWG4Particle * p1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i1)) ;
6175da48 1751 if(fSameSM && GetModuleNumber(p1)!=module1) continue;
1752
1753 //Get kinematics of cluster and (super) module of this cluster
7e7694bb 1754 TLorentzVector photon1(p1->Px(),p1->Py(),p1->Pz(),p1->E());
af7b3903 1755 module1 = GetModuleNumber(p1);
6175da48 1756
1757 //---------------------------------
1758 //First loop on photons/clusters
1759 //---------------------------------
7e7694bb 1760 for(Int_t i2=0; i2<nPhot2; i2++){
1761 AliAODPWG4Particle * p2 = (AliAODPWG4Particle*) (ev2->At(i2)) ;
1762
6175da48 1763 //Get kinematics of second cluster and calculate those of the pair
7e7694bb 1764 TLorentzVector photon2(p2->Px(),p2->Py(),p2->Pz(),p2->E());
6175da48 1765 m = (photon1+photon2).M() ;
7e7694bb 1766 Double_t pt = (photon1 + photon2).Pt();
1767 Double_t a = TMath::Abs(p1->E()-p2->E())/(p1->E()+p2->E()) ;
1768
1769 //Check if opening angle is too large or too small compared to what is expected
1770 Double_t angle = photon1.Angle(photon2.Vect());
6175da48 1771 if(fUseAngleEDepCut && !GetNeutralMesonSelection()->IsAngleInWindow((photon1+photon2).E(),angle+0.05)){
1772 if(GetDebug() > 2)
1773 printf("AliAnaPi0::MakeAnalysisFillHistograms() -Mix pair angle %f not in E %f window\n",angle, (photon1+photon2).E());
1774 continue;
1775 }
1776 if(fUseAngleCut && (angle < fAngleCut || angle > fAngleMaxCut)) {
1777 if(GetDebug() > 2)
1778 printf("AliAnaPi0::MakeAnalysisFillHistograms() -Mix pair angle %f < cut %f\n",angle,fAngleCut);
1779 continue;
1780
1781 }
7e7694bb 1782
1783 if(GetDebug() > 2)
1784 printf("AliAnaPi0::MakeAnalysisFillHistograms() - Mixed Event: pT: photon1 %2.2f, photon2 %2.2f; Pair: pT %2.2f, mass %2.3f, a %f2.3\n",
af7b3903 1785 p1->Pt(), p2->Pt(), pt,m,a);
6175da48 1786
af7b3903 1787 //In case we want only pairs in same (super) module, check their origin.
1788 module2 = GetModuleNumber(p2);
6175da48 1789
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) ;
1796 else
1797 fhMiDiffMod[fNModules+2]->Fill(pt,m) ;
1798
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) ;
1806
1807 }
1808 else {
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) ;
1812 }
1813 }
1814
af7b3903 1815 Bool_t ok = kTRUE;
1816 if(fSameSM && module1!=module2) ok=kFALSE;
1817 if(ok){
6175da48 1818
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);
1822
1823 //Fill histograms for different bad channel distance, centrality, assymmetry cut and pid bit
af7b3903 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) ;
398c93cc 1830 if(fMakeInvPtPlots)fhMiInvPt1[index]->Fill(pt,m,1./pt) ;
6175da48 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) ;
1838 }
af7b3903 1839 }
6175da48 1840 }// Fill bad dist histo
af7b3903 1841 }//Asymmetry cut
1842 }// Asymmetry loop
1843 }//PID cut
1844 }//loop for histograms
6175da48 1845
1846 //-----------------------
1847 //Multi cuts analysis
1848 //-----------------------
1849 if(fMultiCutAna){
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());
1860 }
1861 }// pid bit cut loop
1862 }// icell loop
1863 }// pt cut loop
1864 } // Multi cut ana
1865
1866 //Fill histograms with opening angle
1867 fhMixedOpeningAngle ->Fill(pt,angle);
1868 fhMixedCosOpeningAngle->Fill(pt,TMath::Cos(angle));
af7b3903 1869 }//ok
7e7694bb 1870 }// second cluster loop
1871 }//first cluster loop
1872 }//loop on mixed events
1873
6175da48 1874 //--------------------------------------------------------
1875 //Add the current event to the list of events for mixing
1876 //--------------------------------------------------------
7e7694bb 1877 TClonesArray *currentEvent = new TClonesArray(*GetInputAODBranch());
af7b3903 1878 //Add current event to buffer and Remove redundant events
7e7694bb 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)
1883 {
1884 TClonesArray * tmp = (TClonesArray*) (evMixList->Last()) ;
1885 evMixList->RemoveLast() ;
1886 delete tmp ;
1887 }
1888 }
1889 else{ //empty event
1890 delete currentEvent ;
1891 currentEvent=0 ;
477d6cee 1892 }
7e7694bb 1893 }// DoOwnMix
c8fe2783 1894
1c5acb87 1895}
1896
a5cc4f03 1897//________________________________________________________________________
1898void AliAnaPi0::ReadHistograms(TList* outputList)
1899{
50f39b97 1900 // Needed when Terminate is executed in distributed environment
1901 // Refill analysis histograms of this class with corresponding histograms in output list.
1902
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"));
1906
af7b3903 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] ;
6175da48 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] ;
af7b3903 1919 if(!fhReMod) fhReMod = new TH2D*[fNModules] ;
1920 if(!fhReDiffMod)fhReDiffMod = new TH2D*[fNModules+1] ;
ad30b142 1921 if(!fhMiMod) fhMiMod = new TH2D*[fNModules] ;
1922 if(!fhMiDiffMod)fhMiDiffMod = new TH2D*[fNModules+1] ;
6175da48 1923
1924 fhReConv = (TH2D*) outputList->At(index++);
1925 fhMiConv = (TH2D*) outputList->At(index++);
1926 fhReConv2 = (TH2D*) outputList->At(index++);
1927 fhMiConv2 = (TH2D*) outputList->At(index++);
821c8090 1928
50f39b97 1929 for(Int_t ic=0; ic<fNCentrBin; ic++){
af7b3903 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;
1933
1934 fhRe1[ihisto] = (TH2D*) outputList->At(index++);
1935 fhRe2[ihisto] = (TH2D*) outputList->At(index++);
1936 fhRe3[ihisto] = (TH2D*) outputList->At(index++);
6175da48 1937
1938 fhReInvPt1[ihisto] = (TH2D*) outputList->At(index++);
1939 fhReInvPt2[ihisto] = (TH2D*) outputList->At(index++);
1940 fhReInvPt3[ihisto] = (TH2D*) outputList->At(index++);
5ae09196 1941
af7b3903 1942 if(fDoOwnMix){
1943 fhMi1[ihisto] = (TH2D*) outputList->At(index++);
1944 fhMi2[ihisto] = (TH2D*) outputList->At(index++);
1945 fhMi3[ihisto] = (TH2D*) outputList->At(index++);
6175da48 1946
1947 fhMiInvPt1[ihisto] = (TH2D*) outputList->At(index++);
1948 fhMiInvPt2[ihisto] = (TH2D*) outputList->At(index++);
1949 fhMiInvPt3[ihisto] = (TH2D*) outputList->At(index++);
af7b3903 1950 }//Own mix
1951 }//asymmetry loop
1952 }// pid loop
1953 }// centrality loop
1954
1955 fhRePtAsym = (TH2D*)outputList->At(index++);
1956 fhRePtAsymPi0 = (TH2D*)outputList->At(index++);
1957 fhRePtAsymEta = (TH2D*)outputList->At(index++);
eee5fcf1 1958
5ae09196 1959 if(fMultiCutAna){
1960
eee5fcf1 1961 if(!fhRePtNCellAsymCuts) fhRePtNCellAsymCuts = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
1962 if(!fhRePIDBits) fhRePIDBits = new TH2D*[fNPIDBits];
1963
5ae09196 1964 for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
1965 fhRePIDBits[ipid] = (TH2D*) outputList->At(index++);
1966 }// ipid loop
1967
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++);
1972 }// iasym
1973 }// icell loop
1974 }// ipt loop
af7b3903 1975
1976 if(!fhRePtMult) fhRePtMult = new TH3D*[fNAsymCuts] ;
1977 for(Int_t iasym = 0; iasym < fNAsymCuts; iasym++)
1978 fhRePtMult[iasym] = (TH3D*) outputList->At(index++);
5ae09196 1979 }// multi cut analysis
50f39b97 1980
1981 fhEvents = (TH3D *) outputList->At(index++);
1982
af7b3903 1983 fhRealOpeningAngle = (TH2D*) outputList->At(index++);
1984 fhRealCosOpeningAngle = (TH2D*) outputList->At(index++);
6175da48 1985 if(fDoOwnMix){
1986 fhMixedOpeningAngle = (TH2D*) outputList->At(index++);
1987 fhMixedCosOpeningAngle = (TH2D*) outputList->At(index++);
1988 }
af7b3903 1989
50f39b97 1990 //Histograms filled only if MC data is requested
1991 if(IsDataMC() || (GetReader()->GetDataType() == AliCaloTrackReader::kMC) ){
6175da48 1992 fhPrimPi0Pt = (TH1D*) outputList->At(index++);
1993 fhPrimPi0AccPt = (TH1D*) outputList->At(index++);
156549ae 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++);
6175da48 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++);
2009 }
2010
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++);
2026 }
2027 }
2028 }
2029 }
2030 else{
2031 fhMCPi0MassPtTrue = new TH2D*[1];
2032 fhMCPi0PtTruePtRec = new TH2D*[1];
2033 fhMCEtaMassPtTrue = new TH2D*[1];
2034 fhMCEtaPtTruePtRec = new TH2D*[1];
2035
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++);
2040 }
50f39b97 2041 }
2042
6175da48 2043 for(Int_t imod=0; imod < fNModules; imod++){
2044 fhReMod[imod] = (TH2D*) outputList->At(index++);
2045 fhReDiffMod[imod] = (TH2D*) outputList->At(index++);
2046 if(fDoOwnMix){
2047 fhMiMod[imod] = (TH2D*) outputList->At(index++);
2048 fhMiDiffMod[imod] = (TH2D*) outputList->At(index++);
2049 }
2050 }
eee5fcf1 2051
a5cc4f03 2052}
2053
2054
6639984f 2055//____________________________________________________________________________________________________________________________________________________
a5cc4f03 2056void AliAnaPi0::Terminate(TList* outputList)
6639984f 2057{
2058 //Do some calculations and plots from the final histograms.
477d6cee 2059
fbeaf916 2060 printf(" *** %s Terminate:\n", GetName()) ;
50f39b97 2061
a5cc4f03 2062 //Recover histograms from output histograms list, needed for distributed analysis.
2063 ReadHistograms(outputList);
50f39b97 2064
2e557d1c 2065 if (!fhRe1) {
50f39b97 2066 printf("AliAnaPi0::Terminate() - Error: Remote output histograms not imported in AliAnaPi0 object");
2067 return;
2e557d1c 2068 }
50f39b97 2069
a3aebfff 2070 printf("AliAnaPi0::Terminate() Mgg Real : %5.3f , RMS : %5.3f \n", fhRe1[0]->GetMean(), fhRe1[0]->GetRMS() ) ;
5ae09196 2071
2072 const Int_t buffersize = 255;
2073
2074 char nameIM[buffersize];
2075 snprintf(nameIM, buffersize,"AliAnaPi0_%s_cPt",fCalorimeter.Data());
71dd883b 2076 TCanvas * cIM = new TCanvas(nameIM, "", 400, 10, 600, 700) ;
6639984f 2077 cIM->Divide(2, 2);
50f39b97 2078
6639984f 2079 cIM->cd(1) ;
2080 //gPad->SetLogy();
af7b3903 2081 TH1D * hIMAllPt = (TH1D*) fhRe1[0]->ProjectionY(Form("IMPtAll_%s",fCalorimeter.Data()));
6639984f 2082 hIMAllPt->SetLineColor(2);
2083 hIMAllPt->SetTitle("No cut on p_{T, #gamma#gamma} ");
2084 hIMAllPt->Draw();
2085
2086 cIM->cd(2) ;
af7b3903 2087 TH1D * hIMPt5 = (TH1D*) fhRe1[0]->ProjectionY(Form("IMPt0-5_%s",fCalorimeter.Data()),0, fhRe1[0]->GetXaxis()->FindBin(5.));
2244659d 2088// hRe1Pt5->GetXaxis()->SetRangeUser(0,5);
2089// TH1D * hIMPt5 = (TH1D*) hRe1Pt5->Project3D(Form("IMPt5_%s_pz",fCalorimeter.Data()));
6639984f 2090 hIMPt5->SetLineColor(2);
2091 hIMPt5->SetTitle("0 < p_{T, #gamma#gamma} < 5 GeV/c");
2092 hIMPt5->Draw();
2093
2094 cIM->cd(3) ;
af7b3903 2095 TH1D * hIMPt10 = (TH1D*) fhRe1[0]->ProjectionY(Form("IMPt5-10_%s",fCalorimeter.Data()), fhRe1[0]->GetXaxis()->FindBin(5.),fhRe1[0]->GetXaxis()->FindBin(10.));
2244659d 2096// hRe1Pt10->GetXaxis()->SetRangeUser(5,10);
2097// TH1D * hIMPt10 = (TH1D*) hRe1Pt10->Project3D(Form("IMPt10_%s_pz",fCalorimeter.Data()));
6639984f 2098 hIMPt10->SetLineColor(2);
2099 hIMPt10->SetTitle("5 < p_{T, #gamma#gamma} < 10 GeV/c");
2100 hIMPt10->Draw();
2101
2102 cIM->cd(4) ;
af7b3903 2103 TH1D * hIMPt20 = (TH1D*) fhRe1[0]->ProjectionY(Form("IMPt10-20_%s",fCalorimeter.Data()), fhRe1[0]->GetXaxis()->FindBin(10.),fhRe1[0]->GetXaxis()->FindBin(20.));
2244659d 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()));
6639984f 2107 hIMPt20->SetLineColor(2);
2108 hIMPt20->SetTitle("10 < p_{T, #gamma#gamma} < 20 GeV/c");
2109 hIMPt20->Draw();
2110
5ae09196 2111 char nameIMF[buffersize];
2112 snprintf(nameIMF,buffersize,"AliAnaPi0_%s_Mgg.eps",fCalorimeter.Data());
71dd883b 2113 cIM->Print(nameIMF);
6639984f 2114
5ae09196 2115 char namePt[buffersize];
2116 snprintf(namePt,buffersize,"AliAnaPi0_%s_cPt",fCalorimeter.Data());
71dd883b 2117 TCanvas * cPt = new TCanvas(namePt, "", 400, 10, 600, 700) ;
6639984f 2118 cPt->Divide(2, 2);
2119
2120 cPt->cd(1) ;
2121 //gPad->SetLogy();
af7b3903 2122 TH1D * hPt = (TH1D*) fhRe1[0]->ProjectionX(Form("Pt0_%s",fCalorimeter.Data()),-1,-1);
6639984f 2123 hPt->SetLineColor(2);
2124 hPt->SetTitle("No cut on M_{#gamma#gamma} ");
2125 hPt->Draw();
2126
2127 cPt->cd(2) ;
af7b3903 2128 TH1D * hPtIM1 = (TH1D*)fhRe1[0]->ProjectionX(Form("Pt1_%s",fCalorimeter.Data()), fhRe1[0]->GetZaxis()->FindBin(0.05),fhRe1[0]->GetZaxis()->FindBin(0.21));
2244659d 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");
6639984f 2132 hPtIM1->SetLineColor(2);
2133 hPtIM1->SetTitle("0.05 < M_{#gamma#gamma} < 0.21 GeV/c^{2}");
2134 hPtIM1->Draw();
2135
2136 cPt->cd(3) ;
af7b3903 2137 TH1D * hPtIM2 = (TH1D*)fhRe1[0]->ProjectionX(Form("Pt2_%s",fCalorimeter.Data()), fhRe1[0]->GetZaxis()->FindBin(0.09),fhRe1[0]->GetZaxis()->FindBin(0.17));
2244659d 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");
6639984f 2141 hPtIM2->SetLineColor(2);
2142 hPtIM2->SetTitle("0.09 < M_{#gamma#gamma} < 0.17 GeV/c^{2}");
2143 hPtIM2->Draw();
2144
2145 cPt->cd(4) ;
af7b3903 2146 TH1D * hPtIM3 = (TH1D*)fhRe1[0]->ProjectionX(Form("Pt3_%s",fCalorimeter.Data()), fhRe1[0]->GetZaxis()->FindBin(0.11),fhRe1[0]->GetZaxis()->FindBin(0.15));
2244659d 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");
6639984f 2150 hPtIM3->SetLineColor(2);
2151 hPtIM3->SetTitle("0.11 < M_{#gamma#gamma} < 0.15 GeV/c^{2}");
2152 hPtIM3->Draw();
2153
164a1d84 2154 char namePtF[buffersize];
5ae09196 2155 snprintf(namePtF,buffersize,"AliAnaPi0_%s_Pt.eps",fCalorimeter.Data());
71dd883b 2156 cPt->Print(namePtF);
1c5acb87 2157
5ae09196 2158 char line[buffersize] ;
2159 snprintf(line,buffersize,".!tar -zcf %s_%s.tar.gz *.eps", GetName(),fCalorimeter.Data()) ;
6639984f 2160 gROOT->ProcessLine(line);
5ae09196 2161 snprintf(line, buffersize,".!rm -fR AliAnaPi0_%s*.eps",fCalorimeter.Data());
6639984f 2162 gROOT->ProcessLine(line);
2163
71dd883b 2164 printf(" AliAnaPi0::Terminate() - !! All the eps files are in %s_%s.tar.gz !!!\n", GetName(), fCalorimeter.Data());
1c5acb87 2165
6639984f 2166}
c8fe2783 2167 //____________________________________________________________________________________________________________________________________________________
2168Int_t AliAnaPi0::GetEventIndex(AliAODPWG4Particle * part, Double_t * vert)
2169{
f8006433 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
2173
2174 Int_t evtIndex = -1 ;
2175 if(GetReader()->GetDataType()!=AliCaloTrackReader::kMC){
2176
2177 if (GetMixedEvent()){
2178
2179 evtIndex = GetMixedEvent()->EventIndexForCaloCluster(part->GetCaloLabel(0)) ;
2180 GetVertex(vert,evtIndex);
2181
5025c139 2182 if(TMath::Abs(vert[2])> GetZvertexCut())
f8006433 2183 evtIndex = -2 ; //Event can not be used (vertex, centrality,... cuts not fulfilled)
2184 } else {// Single event
2185
2186 GetVertex(vert);
2187
5025c139 2188 if(TMath::Abs(vert[2])> GetZvertexCut())
f8006433 2189 evtIndex = -1 ; //Event can not be used (vertex, centrality,... cuts not fulfilled)
2190 else
2191 evtIndex = 0 ;
c8fe2783 2192 }
0ae57829 2193 }//No MC reader
f8006433 2194 else {
2195 evtIndex = 0;
2196 vert[0] = 0. ;
2197 vert[1] = 0. ;
2198 vert[2] = 0. ;
2199 }
0ae57829 2200
f8006433 2201 return evtIndex ;
c8fe2783 2202}
f8006433 2203