]> 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
isolation task/cuts modifs (N. Arbor)
[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),
8d230fa8 71fhReMod(0x0), fhReSameSideEMCALMod(0x0), fhReSameSectorEMCALMod(0x0), fhReDiffPHOSMod(0x0),
72fhMiMod(0x0), fhMiSameSideEMCALMod(0x0), fhMiSameSectorEMCALMod(0x0), fhMiDiffPHOSMod(0x0),
156549ae 73fhReConv(0x0), fhMiConv(0x0), fhReConv2(0x0), fhMiConv2(0x0),
74fhRe1(0x0), fhMi1(0x0), fhRe2(0x0), fhMi2(0x0), fhRe3(0x0), fhMi3(0x0),
75fhReInvPt1(0x0), fhMiInvPt1(0x0), fhReInvPt2(0x0), fhMiInvPt2(0x0), fhReInvPt3(0x0), fhMiInvPt3(0x0),
6175da48 76fhRePtNCellAsymCuts(0x0), fhRePtNCellAsymCutsSM0(0x0), fhRePtNCellAsymCutsSM1(0x0), fhRePtNCellAsymCutsSM2(0x0), fhRePtNCellAsymCutsSM3(0x0), fhMiPtNCellAsymCuts(0x0),
77fhRePIDBits(0x0),fhRePtMult(0x0), fhRePtAsym(0x0), fhRePtAsymPi0(0x0),fhRePtAsymEta(0x0),
745f04da 78fhEvents(0x0), fhCentrality(0x0),
79fhRealOpeningAngle(0x0),fhRealCosOpeningAngle(0x0), fhMixedOpeningAngle(0x0),fhMixedCosOpeningAngle(0x0),
6175da48 80fhPrimPi0Pt(0x0), fhPrimPi0AccPt(0x0), fhPrimPi0Y(0x0), fhPrimPi0AccY(0x0), fhPrimPi0Phi(0x0), fhPrimPi0AccPhi(0x0),
08a56f5f 81fhPrimPi0OpeningAngle(0x0), fhPrimPi0CosOpeningAngle(0x0),
6175da48 82fhPrimEtaPt(0x0), fhPrimEtaAccPt(0x0), fhPrimEtaY(0x0), fhPrimEtaAccY(0x0), fhPrimEtaPhi(0x0), fhPrimEtaAccPhi(0x0),
08a56f5f 83fhPrimPi0PtOrigin(0x0), fhPrimEtaPtOrigin(0x0),
6175da48 84fhMCOrgMass(),fhMCOrgAsym(), fhMCOrgDeltaEta(),fhMCOrgDeltaPhi(),
08a56f5f 85fhMCPi0MassPtRec(), fhMCPi0MassPtTrue(), fhMCPi0PtTruePtRec(), fhMCEtaMassPtRec(), fhMCEtaMassPtTrue(), fhMCEtaPtTruePtRec(),
86fhMCPi0PtOrigin(0x0), fhMCEtaPtOrigin(0x0)
1c5acb87 87{
88//Default Ctor
89 InitParameters();
90
91}
1c5acb87 92
93//________________________________________________________________________________________________________________________________________________
94AliAnaPi0::~AliAnaPi0() {
477d6cee 95 // Remove event containers
7e7694bb 96
97 if(fDoOwnMix && fEventsList){
477d6cee 98 for(Int_t ic=0; ic<fNCentrBin; ic++){
5025c139 99 for(Int_t iz=0; iz<GetNZvertBin(); iz++){
100 for(Int_t irp=0; irp<GetNRPBin(); irp++){
101 fEventsList[ic*GetNZvertBin()*GetNRPBin()+iz*GetNRPBin()+irp]->Delete() ;
102 delete fEventsList[ic*GetNZvertBin()*GetNRPBin()+iz*GetNRPBin()+irp] ;
7e7694bb 103 }
477d6cee 104 }
105 }
106 delete[] fEventsList;
107 fEventsList=0 ;
108 }
591cc579 109
1c5acb87 110}
111
112//________________________________________________________________________________________________________________________________________________
113void AliAnaPi0::InitParameters()
114{
115//Init parameters when first called the analysis
116//Set default parameters
a3aebfff 117 SetInputAODName("PWG4Particle");
118
119 AddToHistogramsName("AnaPi0_");
6921fa00 120 fNModules = 12; // set maximum to maximum number of EMCAL modules
477d6cee 121 fNCentrBin = 1;
5025c139 122// fNZvertBin = 1;
123// fNrpBin = 1;
477d6cee 124 fNmaxMixEv = 10;
5025c139 125
477d6cee 126 fCalorimeter = "PHOS";
50f39b97 127 fUseAngleCut = kFALSE;
6175da48 128 fUseAngleEDepCut = kFALSE;
129 fAngleCut = 0.;
130 fAngleMaxCut = TMath::Pi();
131
5ae09196 132 fMultiCutAna = kFALSE;
133
134 fNPtCuts = 3;
af7b3903 135 fPtCuts[0] = 0.; fPtCuts[1] = 0.3; fPtCuts[2] = 0.5;
136 for(Int_t i = fNPtCuts; i < 10; i++)fPtCuts[i] = 0.;
5ae09196 137
9c59b5fe 138 fNAsymCuts = 4;
139 fAsymCuts[0] = 1.; fAsymCuts[1] = 0.8; fAsymCuts[2] = 0.6; fAsymCuts[3] = 0.1;
af7b3903 140 for(Int_t i = fNAsymCuts; i < 10; i++)fAsymCuts[i] = 0.;
141
5ae09196 142 fNCellNCuts = 3;
af7b3903 143 fCellNCuts[0] = 0; fCellNCuts[1] = 1; fCellNCuts[2] = 2;
021c573a 144 for(Int_t i = fNCellNCuts; i < 10; i++)fCellNCuts[i] = 0;
af7b3903 145
146 fNPIDBits = 2;
147 fPIDBits[0] = 0; fPIDBits[1] = 2; // fPIDBits[2] = 4; fPIDBits[3] = 6;// check, no cut, dispersion, neutral, dispersion&&neutral
021c573a 148 for(Int_t i = fNPIDBits; i < 10; i++)fPIDBits[i] = 0;
af7b3903 149
1c5acb87 150}
1c5acb87 151
0c1383b5 152
153//________________________________________________________________________________________________________________________________________________
154TObjString * AliAnaPi0::GetAnalysisCuts()
155{
af7b3903 156 //Save parameters used for analysis
157 TString parList ; //this will be list of parameters used for this analysis.
158 const Int_t buffersize = 255;
159 char onePar[buffersize] ;
160 snprintf(onePar,buffersize,"--- AliAnaPi0 ---\n") ;
161 parList+=onePar ;
162 snprintf(onePar,buffersize,"Number of bins in Centrality: %d \n",fNCentrBin) ;
163 parList+=onePar ;
164 snprintf(onePar,buffersize,"Number of bins in Z vert. pos: %d \n",GetNZvertBin()) ;
165 parList+=onePar ;
166 snprintf(onePar,buffersize,"Number of bins in Reac. Plain: %d \n",GetNRPBin()) ;
167 parList+=onePar ;
168 snprintf(onePar,buffersize,"Depth of event buffer: %d \n",fNmaxMixEv) ;
169 parList+=onePar ;
6175da48 170 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",
171 fSameSM, fUseTrackMultBins, fUsePhotonMultBins, fUseAverClusterEBins, fUseAverCellEBins, fMakeInvPtPlots) ;
172 parList+=onePar ;
173 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 174 parList+=onePar ;
175 snprintf(onePar,buffersize," Asymmetry cuts: n = %d, asymmetry < ",fNAsymCuts) ;
176 for(Int_t i = 0; i < fNAsymCuts; i++) snprintf(onePar,buffersize,"%s %2.2f;",onePar,fAsymCuts[i]);
177 parList+=onePar ;
178 snprintf(onePar,buffersize," PID selection bits: n = %d, PID bit =\n",fNPIDBits) ;
179 for(Int_t i = 0; i < fNPIDBits; i++) snprintf(onePar,buffersize,"%s %d;",onePar,fPIDBits[i]);
180 parList+=onePar ;
181 snprintf(onePar,buffersize,"Cuts: \n") ;
182 parList+=onePar ;
183 snprintf(onePar,buffersize,"Z vertex position: -%f < z < %f \n",GetZvertexCut(),GetZvertexCut()) ;
184 parList+=onePar ;
185 snprintf(onePar,buffersize,"Calorimeter: %s \n",fCalorimeter.Data()) ;
186 parList+=onePar ;
187 snprintf(onePar,buffersize,"Number of modules: %d \n",fNModules) ;
188 parList+=onePar ;
db2bf6fd 189 if(fMultiCutAna){
190 snprintf(onePar, buffersize," pT cuts: n = %d, pt > ",fNPtCuts) ;
191 for(Int_t i = 0; i < fNPtCuts; i++) snprintf(onePar,buffersize,"%s %2.2f;",onePar,fPtCuts[i]);
192 parList+=onePar ;
db2bf6fd 193 snprintf(onePar,buffersize, " N cell in cluster cuts: n = %d, nCell > ",fNCellNCuts) ;
194 for(Int_t i = 0; i < fNCellNCuts; i++) snprintf(onePar,buffersize,"%s %d;",onePar,fCellNCuts[i]);
195 parList+=onePar ;
db2bf6fd 196 }
197
af7b3903 198 return new TObjString(parList) ;
0c1383b5 199}
200
2e557d1c 201//________________________________________________________________________________________________________________________________________________
202TList * AliAnaPi0::GetCreateOutputObjects()
203{
477d6cee 204 // Create histograms to be saved in output file and
205 // store them in fOutputContainer
206
207 //create event containers
5025c139 208 fEventsList = new TList*[fNCentrBin*GetNZvertBin()*GetNRPBin()] ;
1c5acb87 209
477d6cee 210 for(Int_t ic=0; ic<fNCentrBin; ic++){
5025c139 211 for(Int_t iz=0; iz<GetNZvertBin(); iz++){
212 for(Int_t irp=0; irp<GetNRPBin(); irp++){
213 fEventsList[ic*GetNZvertBin()*GetNRPBin()+iz*GetNRPBin()+irp] = new TList() ;
af7b3903 214 fEventsList[ic*GetNZvertBin()*GetNRPBin()+iz*GetNRPBin()+irp]->SetOwner(kFALSE);
477d6cee 215 }
216 }
217 }
7e7694bb 218
477d6cee 219 TList * outputContainer = new TList() ;
220 outputContainer->SetName(GetName());
6921fa00 221
8d230fa8 222 fhReMod = new TH2D*[fNModules] ;
223 fhMiMod = new TH2D*[fNModules] ;
224
225 if(fCalorimeter == "PHOS"){
226 fhReDiffPHOSMod = new TH2D*[fNModules] ;
227 fhMiDiffPHOSMod = new TH2D*[fNModules] ;
228 }
229 else{
230 fhReSameSectorEMCALMod = new TH2D*[fNModules/2] ;
231 fhReSameSideEMCALMod = new TH2D*[fNModules-2] ;
232 fhMiSameSectorEMCALMod = new TH2D*[fNModules/2] ;
233 fhMiSameSideEMCALMod = new TH2D*[fNModules-2] ;
234 }
6175da48 235
af7b3903 236
237 fhRe1 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
af7b3903 238 fhMi1 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
6175da48 239 if(fFillBadDistHisto){
240 fhRe2 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
241 fhRe3 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
242 fhMi2 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
243 fhMi3 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
244 }
245 if(fMakeInvPtPlots) {
398c93cc 246 fhReInvPt1 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
398c93cc 247 fhMiInvPt1 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
6175da48 248 if(fFillBadDistHisto){
249 fhReInvPt2 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
250 fhReInvPt3 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
251 fhMiInvPt2 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
252 fhMiInvPt3 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
253 }
398c93cc 254 }
6175da48 255
5ae09196 256 const Int_t buffersize = 255;
257 char key[buffersize] ;
258 char title[buffersize] ;
477d6cee 259
5a2dbc3c 260 Int_t nptbins = GetHistoPtBins();
261 Int_t nphibins = GetHistoPhiBins();
262 Int_t netabins = GetHistoEtaBins();
263 Float_t ptmax = GetHistoPtMax();
264 Float_t phimax = GetHistoPhiMax();
265 Float_t etamax = GetHistoEtaMax();
266 Float_t ptmin = GetHistoPtMin();
267 Float_t phimin = GetHistoPhiMin();
268 Float_t etamin = GetHistoEtaMin();
269
270 Int_t nmassbins = GetHistoMassBins();
271 Int_t nasymbins = GetHistoAsymmetryBins();
272 Float_t massmax = GetHistoMassMax();
273 Float_t asymmax = GetHistoAsymmetryMax();
274 Float_t massmin = GetHistoMassMin();
275 Float_t asymmin = GetHistoAsymmetryMin();
af7b3903 276 Int_t ntrmbins = GetHistoTrackMultiplicityBins();
277 Int_t ntrmmax = GetHistoTrackMultiplicityMax();
278 Int_t ntrmmin = GetHistoTrackMultiplicityMin();
6175da48 279
280 fhAverTotECluster = new TH1F("hAverTotECluster","hAverTotECluster",200,0,50) ;
281 fhAverTotECluster->SetXTitle("E_{cluster, aver. SM} (GeV)");
282 outputContainer->Add(fhAverTotECluster) ;
283
284 fhAverTotECell = new TH1F("hAverTotECell","hAverTotECell",200,0,50) ;
285 fhAverTotECell->SetXTitle("E_{cell, aver. SM} (GeV)");
286 outputContainer->Add(fhAverTotECell) ;
287
156549ae 288 fhAverTotECellvsCluster = new TH2F("hAverTotECellvsCluster","hAverTotECellvsCluster",200,0,50,200,0,50) ;
289 fhAverTotECellvsCluster->SetYTitle("E_{cell, aver. SM} (GeV)");
290 fhAverTotECellvsCluster->SetXTitle("E_{cluster, aver. SM} (GeV)");
291 outputContainer->Add(fhAverTotECellvsCluster) ;
292
293 fhEDensityCluster = new TH1F("hEDensityCluster","hEDensityCluster",200,0,50) ;
294 fhEDensityCluster->SetXTitle("#Sigma E_{cluster} / N_{cluster} (GeV)");
295 outputContainer->Add(fhEDensityCluster) ;
296
297 fhEDensityCell = new TH1F("hEDensityCell","hEDensityCell",200,0,50) ;
298 fhEDensityCell->SetXTitle("#Sigma E_{cell} / N_{cell} (GeV)");
299 outputContainer->Add(fhEDensityCell) ;
300
301 fhEDensityCellvsCluster = new TH2F("hEDensityCellvsCluster","hEDensityCellvsCluster",200,0,50,200,0,50) ;
302 fhEDensityCellvsCluster->SetYTitle("#Sigma E_{cell} / N_{cell} (GeV)");
303 fhEDensityCellvsCluster->SetXTitle("#Sigma E_{cluster} / N_{cluster} (GeV)");
304 outputContainer->Add(fhEDensityCellvsCluster) ;
305
306// fhClusterPairDist = new TH1F("hClusterPairDist","Distance between clusters",250,0,750) ;
307// fhClusterPairDist->SetXTitle("#sqrt{(x_{1}-x_{2})^2+(z_{1}-z_{2})^2} (cm)");
308// outputContainer->Add(fhClusterPairDist) ;
309//
310// fhClusterPairDistWeight = new TH1F("hClusterPairDistWeighted","Distance between clusters, weighted by pair energy",200,0,400) ;
311// 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)");
312// outputContainer->Add(fhClusterPairDistWeight) ;
313//
314// fhAverClusterPairDist = new TH1F("hAverClusterPairDist","Average distance between clusters",250,0,750) ;
315// fhAverClusterPairDist->SetXTitle("#Sigma (#sqrt{(x_{1}-x_{2})^{2}+(z_{1}-z_{2})^{2}}) / N_{pairs} (cm)");
316// outputContainer->Add(fhAverClusterPairDist) ;
317//
318// fhAverClusterPairDistWeight = new TH1F("hAverClusterPairDistWeighted","Average distance between clusters, weighted by pair energy",200,0,400) ;
319// 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)");
320// outputContainer->Add(fhAverClusterPairDistWeight) ;
321//
322// fhAverClusterPairDistvsAverE = new TH2F("hAverClusterPairDistvsAverE","Average distance between clusters",250,0,750,200,0,50) ;
323// fhAverClusterPairDistvsAverE->SetXTitle("#Sigma (#sqrt{(x_{1}-x_{2})^{2}+(z_{1}-z_{2})^{2}}) / N_{pairs} (cm)");
324// fhAverClusterPairDistvsAverE->SetYTitle("#Sigma E_{cluster} / N_{cluster} (GeV)");
325// outputContainer->Add(fhAverClusterPairDistvsAverE) ;
326//
327// fhAverClusterPairDistWeightvsAverE = new TH2F("hAverClusterPairDistWeightedvsAverE","Average distance between clusters, weighted by pair energy",200,0,400,200,0,50) ;
328// 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)");
329// fhAverClusterPairDistWeightvsAverE->SetYTitle("#Sigma E_{cluster} / N_{cluster} (GeV)");
330// outputContainer->Add(fhAverClusterPairDistWeightvsAverE) ;
331
332// fhAverClusterPairDistvsN = new TH2F("hAverClusterPairDistvsN","Average distance between clusters",250,0,750,200,0,50) ;
333// fhAverClusterPairDistvsN->SetXTitle("#Sigma (#sqrt{(x_{1}-x_{2})^{2}+(z_{1}-z_{2})^{2}}) / N_{pairs} (cm)");
334// fhAverClusterPairDistvsN->SetYTitle("N_{cluster}");
335// outputContainer->Add(fhAverClusterPairDistvsN) ;
336//
337// fhAverClusterPairDistWeightvsN = new TH2F("hAverClusterPairDistWeightedvsN","Average distance between clusters, weighted by pair energy",200,0,400,200,0,50) ;
338// 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)");
339// fhAverClusterPairDistWeightvsN->SetYTitle("N_{cluster}");
340// outputContainer->Add(fhAverClusterPairDistWeightvsN) ;
341
342// fhMaxEvsClustMult = new TH2F("hMaxEvsClustMult","",nptbins,ptmin,ptmax,50,0,50) ;
343// fhMaxEvsClustMult->SetXTitle("E_{max}");
344// fhMaxEvsClustMult->SetYTitle("N_{cluster}");
345// outputContainer->Add(fhMaxEvsClustMult) ;
346//
347// fhMaxEvsClustEDen = new TH2F("hMaxEvsClustEDen","",nptbins,ptmin,ptmax,200,0,50) ;
348// fhMaxEvsClustEDen->SetXTitle("E_{max}");
349// fhMaxEvsClustEDen->SetYTitle("#Sigma E_{cluster} / N_{cluster} (GeV)");
350// outputContainer->Add(fhMaxEvsClustEDen) ;
351
6175da48 352 fhReConv = new TH2D("hReConv","Real Pair with one recombined conversion ",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
353 fhReConv->SetXTitle("p_{T} (GeV/c)");
354 fhReConv->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
355 outputContainer->Add(fhReConv) ;
356
6175da48 357 fhReConv2 = new TH2D("hReConv2","Real Pair with 2 recombined conversion ",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
358 fhReConv2->SetXTitle("p_{T} (GeV/c)");
359 fhReConv2->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
360 outputContainer->Add(fhReConv2) ;
361
156549ae 362 if(fDoOwnMix){
363 fhMiConv = new TH2D("hMiConv","Mixed Pair with one recombined conversion ",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
364 fhMiConv->SetXTitle("p_{T} (GeV/c)");
365 fhMiConv->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
366 outputContainer->Add(fhMiConv) ;
367
368 fhMiConv2 = new TH2D("hMiConv2","Mixed Pair with 2 recombined conversion ",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
369 fhMiConv2->SetXTitle("p_{T} (GeV/c)");
370 fhMiConv2->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
371 outputContainer->Add(fhMiConv2) ;
372 }
6175da48 373
477d6cee 374 for(Int_t ic=0; ic<fNCentrBin; ic++){
6175da48 375 for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
376 for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
377 Int_t index = ((ic*fNPIDBits)+ipid)*fNAsymCuts + iasym;
378 //printf("cen %d, pid %d, asy %d, Index %d\n",ic,ipid,iasym,index);
398c93cc 379 //Distance to bad module 1
6175da48 380 snprintf(key, buffersize,"hRe_cen%d_pidbit%d_asy%d_dist1",ic,ipid,iasym) ;
398c93cc 381 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 1",
382 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
6175da48 383 fhRe1[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
384 fhRe1[index]->SetXTitle("p_{T} (GeV/c)");
385 fhRe1[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
386 //printf("name: %s\n ",fhRe1[index]->GetName());
387 outputContainer->Add(fhRe1[index]) ;
af7b3903 388
6175da48 389 if(fFillBadDistHisto){
398c93cc 390 //Distance to bad module 2
6175da48 391 snprintf(key, buffersize,"hRe_cen%d_pidbit%d_asy%d_dist2",ic,ipid,iasym) ;
392 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 2",
398c93cc 393 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
6175da48 394 fhRe2[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
395 fhRe2[index]->SetXTitle("p_{T} (GeV/c)");
396 fhRe2[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
397 outputContainer->Add(fhRe2[index]) ;
398c93cc 398
399 //Distance to bad module 3
6175da48 400 snprintf(key, buffersize,"hRe_cen%d_pidbit%d_asy%d_dist3",ic,ipid,iasym) ;
401 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 3",
398c93cc 402 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
6175da48 403 fhRe3[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
404 fhRe3[index]->SetXTitle("p_{T} (GeV/c)");
405 fhRe3[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
406 outputContainer->Add(fhRe3[index]) ;
398c93cc 407 }
6175da48 408
409 //Inverse pT
410 if(fMakeInvPtPlots){
411 //Distance to bad module 1
412 snprintf(key, buffersize,"hReInvPt_cen%d_pidbit%d_asy%d_dist1",ic,ipid,iasym) ;
413 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 1",
414 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
415 fhReInvPt1[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
416 fhReInvPt1[index]->SetXTitle("p_{T} (GeV/c)");
417 fhReInvPt1[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
418 outputContainer->Add(fhReInvPt1[index]) ;
419
420 if(fFillBadDistHisto){
421 //Distance to bad module 2
422 snprintf(key, buffersize,"hReInvPt_cen%d_pidbit%d_asy%d_dist2",ic,ipid,iasym) ;
423 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 2",
424 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
425 fhReInvPt2[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
426 fhReInvPt2[index]->SetXTitle("p_{T} (GeV/c)");
427 fhReInvPt2[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
428 outputContainer->Add(fhReInvPt2[index]) ;
429
430 //Distance to bad module 3
431 snprintf(key, buffersize,"hReInvPt_cen%d_pidbit%d_asy%d_dist3",ic,ipid,iasym) ;
432 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 3",
433 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
434 fhReInvPt3[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
435 fhReInvPt3[index]->SetXTitle("p_{T} (GeV/c)");
436 fhReInvPt3[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
437 outputContainer->Add(fhReInvPt3[index]) ;
438 }
439 }
440 if(fDoOwnMix){
441 //Distance to bad module 1
442 snprintf(key, buffersize,"hMi_cen%d_pidbit%d_asy%d_dist1",ic,ipid,iasym) ;
443 snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 1",
444 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
445 fhMi1[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
446 fhMi1[index]->SetXTitle("p_{T} (GeV/c)");
447 fhMi1[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
448 outputContainer->Add(fhMi1[index]) ;
449 if(fFillBadDistHisto){
450 //Distance to bad module 2
451 snprintf(key, buffersize,"hMi_cen%d_pidbit%d_asy%d_dist2",ic,ipid,iasym) ;
452 snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 2",
453 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
454 fhMi2[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
455 fhMi2[index]->SetXTitle("p_{T} (GeV/c)");
456 fhMi2[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
457 outputContainer->Add(fhMi2[index]) ;
458
459 //Distance to bad module 3
460 snprintf(key, buffersize,"hMi_cen%d_pidbit%d_asy%d_dist3",ic,ipid,iasym) ;
461 snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 3",
462 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
463 fhMi3[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
464 fhMi3[index]->SetXTitle("p_{T} (GeV/c)");
465 fhMi3[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
466 outputContainer->Add(fhMi3[index]) ;
467 }
468 //Inverse pT
469 if(fMakeInvPtPlots){
470 //Distance to bad module 1
471 snprintf(key, buffersize,"hMiInvPt_cen%d_pidbit%d_asy%d_dist1",ic,ipid,iasym) ;
472 snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 1",
473 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
474 fhMiInvPt1[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
475 fhMiInvPt1[index]->SetXTitle("p_{T} (GeV/c)");
476 fhMiInvPt1[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
477 outputContainer->Add(fhMiInvPt1[index]) ;
478 if(fFillBadDistHisto){
479 //Distance to bad module 2
480 snprintf(key, buffersize,"hMiInvPt_cen%d_pidbit%d_asy%d_dist2",ic,ipid,iasym) ;
481 snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 2",
482 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
483 fhMiInvPt2[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
484 fhMiInvPt2[index]->SetXTitle("p_{T} (GeV/c)");
485 fhMiInvPt2[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
486 outputContainer->Add(fhMiInvPt2[index]) ;
487
488 //Distance to bad module 3
489 snprintf(key, buffersize,"hMiInvPt_cen%d_pidbit%d_asy%d_dist3",ic,ipid,iasym) ;
490 snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f,dist bad 3",
491 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
492 fhMiInvPt3[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
493 fhMiInvPt3[index]->SetXTitle("p_{T} (GeV/c)");
494 fhMiInvPt3[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
495 outputContainer->Add(fhMiInvPt3[index]) ;
496 }
497 }
498 }
499 }
7e7694bb 500 }
1c5acb87 501 }
477d6cee 502
9c59b5fe 503 fhRePtAsym = new TH2D("hRePtAsym","Asymmetry vs pt, for pairs",nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax) ;
af7b3903 504 fhRePtAsym->SetXTitle("p_{T} (GeV/c)");
505 fhRePtAsym->SetYTitle("Asymmetry");
506 outputContainer->Add(fhRePtAsym);
507
9c59b5fe 508 fhRePtAsymPi0 = new TH2D("hRePtAsymPi0","Asymmetry vs pt, for pairs close to #pi^{0} mass",nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax) ;
af7b3903 509 fhRePtAsymPi0->SetXTitle("p_{T} (GeV/c)");
510 fhRePtAsymPi0->SetYTitle("Asymmetry");
511 outputContainer->Add(fhRePtAsymPi0);
512
9c59b5fe 513 fhRePtAsymEta = new TH2D("hRePtAsymEta","Asymmetry vs pt, for pairs close to #eta mass",nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax) ;
af7b3903 514 fhRePtAsymEta->SetXTitle("p_{T} (GeV/c)");
515 fhRePtAsymEta->SetYTitle("Asymmetry");
516 outputContainer->Add(fhRePtAsymEta);
517
5ae09196 518 if(fMultiCutAna){
519
520 fhRePIDBits = new TH2D*[fNPIDBits];
521 for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
522 snprintf(key, buffersize,"hRe_pidbit%d",ipid) ;
523 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for PIDBit=%d",fPIDBits[ipid]) ;
524 fhRePIDBits[ipid] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
af7b3903 525 fhRePIDBits[ipid]->SetXTitle("p_{T} (GeV/c)");
526 fhRePIDBits[ipid]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
5ae09196 527 outputContainer->Add(fhRePIDBits[ipid]) ;
528 }// pid bit loop
529
6175da48 530 fhRePtNCellAsymCuts = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
531 fhRePtNCellAsymCutsSM0 = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
532 fhRePtNCellAsymCutsSM1 = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
533 fhRePtNCellAsymCutsSM2 = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
534 fhRePtNCellAsymCutsSM3 = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
535 fhMiPtNCellAsymCuts = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
5ae09196 536 for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
537 for(Int_t icell=0; icell<fNCellNCuts; icell++){
538 for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
539 snprintf(key, buffersize,"hRe_pt%d_cell%d_asym%d",ipt,icell,iasym) ;
af7b3903 540 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 541 Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
542 //printf("ipt %d, icell %d, iassym %d, index %d\n",ipt, icell, iasym, index);
543 fhRePtNCellAsymCuts[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
af7b3903 544 fhRePtNCellAsymCuts[index]->SetXTitle("p_{T} (GeV/c)");
545 fhRePtNCellAsymCuts[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
5ae09196 546 outputContainer->Add(fhRePtNCellAsymCuts[index]) ;
6175da48 547
548 snprintf(key, buffersize,"hRe_pt%d_cell%d_asym%d_SM0",ipt,icell,iasym) ;
549 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]) ;
550 fhRePtNCellAsymCutsSM0[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
551 fhRePtNCellAsymCutsSM0[index]->SetXTitle("p_{T} (GeV/c)");
552 fhRePtNCellAsymCutsSM0[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
553 outputContainer->Add(fhRePtNCellAsymCutsSM0[index]) ;
554
555 snprintf(key, buffersize,"hRe_pt%d_cell%d_asym%d_SM1",ipt,icell,iasym) ;
556 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]) ;
557 fhRePtNCellAsymCutsSM1[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
558 fhRePtNCellAsymCutsSM1[index]->SetXTitle("p_{T} (GeV/c)");
559 fhRePtNCellAsymCutsSM1[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
560 outputContainer->Add(fhRePtNCellAsymCutsSM1[index]) ;
561
562 snprintf(key, buffersize,"hRe_pt%d_cell%d_asym%d_SM2",ipt,icell,iasym) ;
563 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]) ;
564 fhRePtNCellAsymCutsSM2[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
565 fhRePtNCellAsymCutsSM2[index]->SetXTitle("p_{T} (GeV/c)");
566 fhRePtNCellAsymCutsSM2[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
567 outputContainer->Add(fhRePtNCellAsymCutsSM2[index]) ;
568
569 snprintf(key, buffersize,"hRe_pt%d_cell%d_asym%d_SM3",ipt,icell,iasym) ;
570 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]) ;
571 fhRePtNCellAsymCutsSM3[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
572 fhRePtNCellAsymCutsSM3[index]->SetXTitle("p_{T} (GeV/c)");
573 fhRePtNCellAsymCutsSM3[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
574 outputContainer->Add(fhRePtNCellAsymCutsSM3[index]) ;
575
576 snprintf(key, buffersize,"hMi_pt%d_cell%d_asym%d",ipt,icell,iasym) ;
577 snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for pt >%2.2f, ncell>%d and asym >%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]) ;
578 fhMiPtNCellAsymCuts[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
579 fhMiPtNCellAsymCuts[index]->SetXTitle("p_{T} (GeV/c)");
580 fhMiPtNCellAsymCuts[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
581 outputContainer->Add(fhMiPtNCellAsymCuts[index]) ;
582
5ae09196 583 }
584 }
585 }
821c8090 586
af7b3903 587 fhRePtMult = new TH3D*[fNAsymCuts] ;
588 for(Int_t iasym = 0; iasym<fNAsymCuts; iasym++){
589 fhRePtMult[iasym] = new TH3D(Form("hRePtMult_asym%d",iasym),Form("(p_{T},C,M)_{#gamma#gamma}, A<%1.2f",fAsymCuts[iasym]),
590 nptbins,ptmin,ptmax,ntrmbins,ntrmmin,ntrmmax,nmassbins,massmin,massmax);
591 fhRePtMult[iasym]->SetXTitle("p_{T} (GeV/c)");
592 fhRePtMult[iasym]->SetYTitle("Track multiplicity");
593 fhRePtMult[iasym]->SetZTitle("m_{#gamma,#gamma} (GeV/c^{2})");
594 outputContainer->Add(fhRePtMult[iasym]) ;
595 }
596
5ae09196 597 }// multi cuts analysis
598
477d6cee 599 fhEvents=new TH3D("hEvents","Number of events",fNCentrBin,0.,1.*fNCentrBin,
5025c139 600 GetNZvertBin(),0.,1.*GetNZvertBin(),GetNRPBin(),0.,1.*GetNRPBin()) ;
745f04da 601
602 fhEvents->SetXTitle("Centrality bin");
603 fhEvents->SetYTitle("Z vertex bin bin");
604 fhEvents->SetZTitle("RP bin");
477d6cee 605 outputContainer->Add(fhEvents) ;
50f39b97 606
745f04da 607 fhCentrality=new TH1D("hCentralityBin","Number of events in centrality bin",fNCentrBin*10,0.,1.*fNCentrBin) ;
608 fhCentrality->SetXTitle("Centrality bin");
609 outputContainer->Add(fhCentrality) ;
610
50f39b97 611 fhRealOpeningAngle = new TH2D
6175da48 612 ("hRealOpeningAngle","Angle between all #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,300,0,TMath::Pi());
50f39b97 613 fhRealOpeningAngle->SetYTitle("#theta(rad)");
614 fhRealOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
615 outputContainer->Add(fhRealOpeningAngle) ;
7e7694bb 616
50f39b97 617 fhRealCosOpeningAngle = new TH2D
6175da48 618 ("hRealCosOpeningAngle","Cosinus of angle between all #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,100,0,1);
50f39b97 619 fhRealCosOpeningAngle->SetYTitle("cos (#theta) ");
620 fhRealCosOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
621 outputContainer->Add(fhRealCosOpeningAngle) ;
622
6175da48 623 if(fDoOwnMix){
624
625 fhMixedOpeningAngle = new TH2D
626 ("hMixedOpeningAngle","Angle between all #gamma pair vs E_{#pi^{0}}, Mixed pairs",nptbins,ptmin,ptmax,300,0,TMath::Pi());
627 fhMixedOpeningAngle->SetYTitle("#theta(rad)");
628 fhMixedOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
629 outputContainer->Add(fhMixedOpeningAngle) ;
630
631 fhMixedCosOpeningAngle = new TH2D
632 ("hMixedCosOpeningAngle","Cosinus of angle between all #gamma pair vs E_{#pi^{0}}, Mixed pairs",nptbins,ptmin,ptmax,100,0,1);
633 fhMixedCosOpeningAngle->SetYTitle("cos (#theta) ");
634 fhMixedCosOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
635 outputContainer->Add(fhMixedCosOpeningAngle) ;
636
637 }
638
477d6cee 639 //Histograms filled only if MC data is requested
0ae57829 640 if(IsDataMC()){
6175da48 641 //Pi0
642 fhPrimPi0Pt = new TH1D("hPrimPi0Pt","Primary pi0 pt",nptbins,ptmin,ptmax) ;
643 fhPrimPi0AccPt = new TH1D("hPrimPi0AccPt","Primary pi0 pt with both photons in acceptance",nptbins,ptmin,ptmax) ;
156549ae 644 fhPrimPi0Pt ->SetXTitle("p_{T} (GeV/c)");
645 fhPrimPi0AccPt->SetXTitle("p_{T} (GeV/c)");
6175da48 646 outputContainer->Add(fhPrimPi0Pt) ;
647 outputContainer->Add(fhPrimPi0AccPt) ;
648
156549ae 649 fhPrimPi0Y = new TH2D("hPrimPi0Rapidity","Rapidity of primary pi0",nptbins,ptmin,ptmax, netabins,etamin,etamax) ;
650 fhPrimPi0Y ->SetYTitle("Rapidity");
651 fhPrimPi0Y ->SetXTitle("p_{T} (GeV/c)");
6175da48 652 outputContainer->Add(fhPrimPi0Y) ;
653
156549ae 654 fhPrimPi0AccY = new TH2D("hPrimPi0AccRapidity","Rapidity of primary pi0",nptbins,ptmin,ptmax, netabins,etamin,etamax) ;
655 fhPrimPi0AccY->SetYTitle("Rapidity");
656 fhPrimPi0AccY->SetXTitle("p_{T} (GeV/c)");
6175da48 657 outputContainer->Add(fhPrimPi0AccY) ;
477d6cee 658
156549ae 659 fhPrimPi0Phi = new TH2D("hPrimPi0Phi","Azimuthal of primary pi0",nptbins,ptmin,ptmax, nphibins,phimin*TMath::RadToDeg(),phimax*TMath::RadToDeg()) ;
660 fhPrimPi0Phi->SetYTitle("#phi (deg)");
661 fhPrimPi0Phi->SetXTitle("p_{T} (GeV/c)");
6175da48 662 outputContainer->Add(fhPrimPi0Phi) ;
477d6cee 663
156549ae 664 fhPrimPi0AccPhi = new TH2D("hPrimPi0AccPhi","Azimuthal of primary pi0 with accepted daughters",nptbins,ptmin,ptmax, nphibins,phimin*TMath::RadToDeg(),phimax*TMath::RadToDeg()) ;
665 fhPrimPi0AccPhi->SetYTitle("#phi (deg)");
666 fhPrimPi0AccPhi->SetXTitle("p_{T} (GeV/c)");
6175da48 667 outputContainer->Add(fhPrimPi0AccPhi) ;
477d6cee 668
6175da48 669 //Eta
670 fhPrimEtaPt = new TH1D("hPrimEtaPt","Primary eta pt",nptbins,ptmin,ptmax) ;
671 fhPrimEtaAccPt = new TH1D("hPrimEtaAccPt","Primary eta pt with both photons in acceptance",nptbins,ptmin,ptmax) ;
156549ae 672 fhPrimEtaPt ->SetXTitle("p_{T} (GeV/c)");
673 fhPrimEtaAccPt->SetXTitle("p_{T} (GeV/c)");
6175da48 674 outputContainer->Add(fhPrimEtaPt) ;
675 outputContainer->Add(fhPrimEtaAccPt) ;
477d6cee 676
156549ae 677 fhPrimEtaY = new TH2D("hPrimEtaRapidity","Rapidity of primary eta",nptbins,ptmin,ptmax,netabins,etamin,etamax) ;
678 fhPrimEtaY->SetYTitle("Rapidity");
679 fhPrimEtaY->SetXTitle("p_{T} (GeV/c)");
6175da48 680 outputContainer->Add(fhPrimEtaY) ;
50f39b97 681
156549ae 682 fhPrimEtaAccY = new TH2D("hPrimEtaAccRapidity","Rapidity of primary eta",nptbins,ptmin,ptmax, netabins,etamin,etamax) ;
683 fhPrimEtaAccY->SetYTitle("Rapidity");
684 fhPrimEtaAccY->SetXTitle("p_{T} (GeV/c)");
6175da48 685 outputContainer->Add(fhPrimEtaAccY) ;
50f39b97 686
156549ae 687 fhPrimEtaPhi = new TH2D("hPrimEtaPhi","Azimuthal of primary eta",nptbins,ptmin,ptmax, nphibins,phimin*TMath::RadToDeg(),phimax*TMath::RadToDeg()) ;
688 fhPrimEtaPhi->SetYTitle("#phi (deg)");
689 fhPrimEtaPhi->SetXTitle("p_{T} (GeV/c)");
6175da48 690 outputContainer->Add(fhPrimEtaPhi) ;
691
156549ae 692 fhPrimEtaAccPhi = new TH2D("hPrimEtaAccPhi","Azimuthal of primary eta with accepted daughters",nptbins,ptmin,ptmax, nphibins,phimin*TMath::RadToDeg(),phimax*TMath::RadToDeg()) ;
693 fhPrimEtaAccPhi->SetYTitle("#phi (deg)");
694 fhPrimEtaAccPhi->SetXTitle("p_{T} (GeV/c)");
6175da48 695 outputContainer->Add(fhPrimEtaAccPhi) ;
696
50f39b97 697
08a56f5f 698 //Prim origin
699 //Pi0
700 fhPrimPi0PtOrigin = new TH2D("hPrimPi0PtOrigin","Primary pi0 pt vs origin",nptbins,ptmin,ptmax,11,0,11) ;
701 fhPrimPi0PtOrigin->SetXTitle("p_{T} (GeV/c)");
702 fhPrimPi0PtOrigin->SetYTitle("Origin");
703 fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(1 ,"Status 21");
704 fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(2 ,"Quark");
705 fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(3 ,"qq Resonances ");
706 fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(4 ,"Resonances");
707 fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(5 ,"#rho");
708 fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(6 ,"#omega");
709 fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(7 ,"K");
710 fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(8 ,"Other");
711 fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(9 ,"#eta");
712 fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(10 ,"#eta prime");
713 outputContainer->Add(fhPrimPi0PtOrigin) ;
714
715 fhMCPi0PtOrigin = new TH2D("hMCPi0PtOrigin","Reconstructed pair from generated pi0 pt vs origin",nptbins,ptmin,ptmax,11,0,11) ;
716 fhMCPi0PtOrigin->SetXTitle("p_{T} (GeV/c)");
717 fhMCPi0PtOrigin->SetYTitle("Origin");
718 fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(1 ,"Status 21");
719 fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(2 ,"Quark");
720 fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(3 ,"qq Resonances");
721 fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(4 ,"Resonances");
722 fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(5 ,"#rho");
723 fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(6 ,"#omega");
724 fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(7 ,"K");
725 fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(8 ,"Other");
726 fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(9 ,"#eta");
727 fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(10 ,"#eta prime");
728 outputContainer->Add(fhMCPi0PtOrigin) ;
729
730 //Eta
731 fhPrimEtaPtOrigin = new TH2D("hPrimEtaPtOrigin","Primary pi0 pt vs origin",nptbins,ptmin,ptmax,7,0,7) ;
732 fhPrimEtaPtOrigin->SetXTitle("p_{T} (GeV/c)");
733 fhPrimEtaPtOrigin->SetYTitle("Origin");
734 fhPrimEtaPtOrigin->GetYaxis()->SetBinLabel(1 ,"Status 21");
735 fhPrimEtaPtOrigin->GetYaxis()->SetBinLabel(2 ,"Quark");
736 fhPrimEtaPtOrigin->GetYaxis()->SetBinLabel(3 ,"qq Resonances");
737 fhPrimEtaPtOrigin->GetYaxis()->SetBinLabel(4 ,"Resonances");
738 fhPrimEtaPtOrigin->GetYaxis()->SetBinLabel(5 ,"Other");
739 fhPrimEtaPtOrigin->GetYaxis()->SetBinLabel(6 ,"#eta prime ");
740
741 outputContainer->Add(fhPrimEtaPtOrigin) ;
742
743 fhMCEtaPtOrigin = new TH2D("hMCEtaPtOrigin","Reconstructed pair from generated pi0 pt vs origin",nptbins,ptmin,ptmax,7,0,7) ;
744 fhMCEtaPtOrigin->SetXTitle("p_{T} (GeV/c)");
745 fhMCEtaPtOrigin->SetYTitle("Origin");
746 fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(1 ,"Status 21");
747 fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(2 ,"Quark");
748 fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(3 ,"qq Resonances");
749 fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(4 ,"Resonances");
750 fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(5 ,"Other");
751 fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(6 ,"#eta prime");
752
753 outputContainer->Add(fhMCEtaPtOrigin) ;
754
755
6175da48 756 fhPrimPi0OpeningAngle = new TH2D
757 ("hPrimPi0OpeningAngle","Angle between all primary #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,100,0,0.5);
758 fhPrimPi0OpeningAngle->SetYTitle("#theta(rad)");
759 fhPrimPi0OpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
760 outputContainer->Add(fhPrimPi0OpeningAngle) ;
761
762 fhPrimPi0CosOpeningAngle = new TH2D
763 ("hPrimPi0CosOpeningAngle","Cosinus of angle between all primary #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,100,-1,1);
764 fhPrimPi0CosOpeningAngle->SetYTitle("cos (#theta) ");
765 fhPrimPi0CosOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
766 outputContainer->Add(fhPrimPi0CosOpeningAngle) ;
767
768 for(Int_t i = 0; i<13; i++){
769 fhMCOrgMass[i] = new TH2D(Form("hMCOrgMass_%d",i),Form("mass vs pt, origin %d",i),nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
770 fhMCOrgMass[i]->SetXTitle("p_{T} (GeV/c)");
771 fhMCOrgMass[i]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
772 outputContainer->Add(fhMCOrgMass[i]) ;
773
774 fhMCOrgAsym[i]= new TH2D(Form("hMCOrgAsym_%d",i),Form("asymmetry vs pt, origin %d",i),nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax) ;
775 fhMCOrgAsym[i]->SetXTitle("p_{T} (GeV/c)");
776 fhMCOrgAsym[i]->SetYTitle("A");
777 outputContainer->Add(fhMCOrgAsym[i]) ;
778
156549ae 779 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 780 fhMCOrgDeltaEta[i]->SetXTitle("p_{T} (GeV/c)");
156549ae 781 fhMCOrgDeltaEta[i]->SetYTitle("#Delta #eta");
6175da48 782 outputContainer->Add(fhMCOrgDeltaEta[i]) ;
783
156549ae 784 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 785 fhMCOrgDeltaPhi[i]->SetXTitle("p_{T} (GeV/c)");
156549ae 786 fhMCOrgDeltaPhi[i]->SetYTitle("#Delta #phi (rad)");
6175da48 787 outputContainer->Add(fhMCOrgDeltaPhi[i]) ;
788
789 }
50f39b97 790
6175da48 791 if(fMultiCutAnaSim){
792 fhMCPi0MassPtTrue = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
793 fhMCPi0MassPtRec = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
794 fhMCPi0PtTruePtRec = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
795 fhMCEtaMassPtRec = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
796 fhMCEtaMassPtTrue = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
797 fhMCEtaPtTruePtRec = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
798 for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
799 for(Int_t icell=0; icell<fNCellNCuts; icell++){
800 for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
801 Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
802
803 fhMCPi0MassPtRec[index] = new TH2D(Form("hMCPi0MassPtRec_pt%d_cell%d_asym%d",ipt,icell,iasym),
804 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]),
805 nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
806 fhMCPi0MassPtRec[index]->SetXTitle("p_{T, reconstructed} (GeV/c)");
807 fhMCPi0MassPtRec[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
808 outputContainer->Add(fhMCPi0MassPtRec[index]) ;
809
810 fhMCPi0MassPtTrue[index] = new TH2D(Form("hMCPi0MassPtTrue_pt%d_cell%d_asym%d",ipt,icell,iasym),
811 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]),
812 nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
813 fhMCPi0MassPtTrue[index]->SetXTitle("p_{T, generated} (GeV/c)");
814 fhMCPi0MassPtTrue[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
815 outputContainer->Add(fhMCPi0MassPtTrue[index]) ;
816
817 fhMCPi0PtTruePtRec[index] = new TH2D(Form("hMCPi0PtTruePtRec_pt%d_cell%d_asym%d",ipt,icell,iasym),
818 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]),
819 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax) ;
820 fhMCPi0PtTruePtRec[index]->SetXTitle("p_{T, generated} (GeV/c)");
821 fhMCPi0PtTruePtRec[index]->SetYTitle("p_{T, reconstructed} (GeV/c)");
822 outputContainer->Add(fhMCPi0PtTruePtRec[index]) ;
823
824 fhMCEtaMassPtRec[index] = new TH2D(Form("hMCEtaMassPtRec_pt%d_cell%d_asym%d",ipt,icell,iasym),
825 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]),
826 nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
827 fhMCEtaMassPtRec[index]->SetXTitle("p_{T, generated} (GeV/c)");
828 fhMCEtaMassPtRec[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
829 outputContainer->Add(fhMCEtaMassPtRec[index]) ;
830
831 fhMCEtaMassPtTrue[index] = new TH2D(Form("hMCEtaMassPtTrue_pt%d_cell%d_asym%d",ipt,icell,iasym),
832 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]),
833 nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
834 fhMCEtaMassPtTrue[index]->SetXTitle("p_{T, generated} (GeV/c)");
835 fhMCEtaMassPtTrue[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
836 outputContainer->Add(fhMCEtaMassPtTrue[index]) ;
837
838 fhMCEtaPtTruePtRec[index] = new TH2D(Form("hMCEtaPtTruePtRec_pt%d_cell%d_asym%d",ipt,icell,iasym),
839 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]),
840 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax) ;
841 fhMCEtaPtTruePtRec[index]->SetXTitle("p_{T, generated} (GeV/c)");
842 fhMCEtaPtTruePtRec[index]->SetYTitle("p_{T, reconstructed} (GeV/c)");
843 outputContainer->Add(fhMCEtaPtTruePtRec[index]) ;
844 }
845 }
846 }
847 }//multi cut ana
848 else {
849 fhMCPi0MassPtTrue = new TH2D*[1];
850 fhMCPi0PtTruePtRec = new TH2D*[1];
851 fhMCEtaMassPtTrue = new TH2D*[1];
852 fhMCEtaPtTruePtRec = new TH2D*[1];
853
854 fhMCPi0MassPtTrue[0] = new TH2D("hMCPi0MassPtTrue","Reconstructed Mass vs generated p_T of true #pi^{0} cluster pairs",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
855 fhMCPi0MassPtTrue[0]->SetXTitle("p_{T, generated} (GeV/c)");
856 fhMCPi0MassPtTrue[0]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
857 outputContainer->Add(fhMCPi0MassPtTrue[0]) ;
858
859 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) ;
860 fhMCPi0PtTruePtRec[0]->SetXTitle("p_{T, generated} (GeV/c)");
861 fhMCPi0PtTruePtRec[0]->SetYTitle("p_{T, reconstructed} (GeV/c)");
862 outputContainer->Add(fhMCPi0PtTruePtRec[0]) ;
863
864 fhMCEtaMassPtTrue[0] = new TH2D("hMCEtaMassPtTrue","Reconstructed Mass vs generated p_T of true #eta cluster pairs",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
865 fhMCEtaMassPtTrue[0]->SetXTitle("p_{T, generated} (GeV/c)");
866 fhMCEtaMassPtTrue[0]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
867 outputContainer->Add(fhMCEtaMassPtTrue[0]) ;
868
869 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) ;
870 fhMCEtaPtTruePtRec[0]->SetXTitle("p_{T, generated} (GeV/c)");
871 fhMCEtaPtTruePtRec[0]->SetYTitle("p_{T, reconstructed} (GeV/c)");
872 outputContainer->Add(fhMCEtaPtTruePtRec[0]) ;
873 }
477d6cee 874 }
50f39b97 875
8d230fa8 876 TString pairnamePHOS[] = {"(0-1)","(0-2)","(1-2)","(0-3)","(0-4)","(1-3)","(1-4)","(2-3)","(2-4)","(3-4)"};
6921fa00 877 for(Int_t imod=0; imod<fNModules; imod++){
50f39b97 878 //Module dependent invariant mass
5ae09196 879 snprintf(key, buffersize,"hReMod_%d",imod) ;
880 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for Module %d",imod) ;
af7b3903 881 fhReMod[imod] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
882 fhReMod[imod]->SetXTitle("p_{T} (GeV/c)");
883 fhReMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
50f39b97 884 outputContainer->Add(fhReMod[imod]) ;
8d230fa8 885 if(fCalorimeter=="PHOS"){
886 snprintf(key, buffersize,"hReDiffPHOSMod_%d",imod) ;
887 snprintf(title, buffersize,"Real pairs PHOS, clusters in different Modules: %s",(pairnamePHOS[imod]).Data()) ;
888 fhReDiffPHOSMod[imod] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
889 fhReDiffPHOSMod[imod]->SetXTitle("p_{T} (GeV/c)");
890 fhReDiffPHOSMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
891 outputContainer->Add(fhReDiffPHOSMod[imod]) ;
892 }
893 else{//EMCAL
894 if(imod<fNModules/2){
895 snprintf(key, buffersize,"hReSameSectorEMCAL_%d",imod) ;
896 snprintf(title, buffersize,"Real pairs EMCAL, clusters in same sector, SM(%d,%d)",imod*2,imod*2+1) ;
897 fhReSameSectorEMCALMod[imod] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
898 fhReSameSectorEMCALMod[imod]->SetXTitle("p_{T} (GeV/c)");
899 fhReSameSectorEMCALMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
900 outputContainer->Add(fhReSameSectorEMCALMod[imod]) ;
901 }
902 if(imod<fNModules-2){
903 snprintf(key, buffersize,"hReSameSideEMCAL_%d",imod) ;
904 snprintf(title, buffersize,"Real pairs EMCAL, clusters in same side SM(%d,%d)",imod, imod+2) ;
905 fhReSameSideEMCALMod[imod] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
906 fhReSameSideEMCALMod[imod]->SetXTitle("p_{T} (GeV/c)");
907 fhReSameSideEMCALMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
908 outputContainer->Add(fhReSameSideEMCALMod[imod]) ;
909 }
910 }//EMCAL
6175da48 911
912 if(fDoOwnMix){
913 snprintf(key, buffersize,"hMiMod_%d",imod) ;
914 snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for Module %d",imod) ;
915 fhMiMod[imod] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
916 fhMiMod[imod]->SetXTitle("p_{T} (GeV/c)");
917 fhMiMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
918 outputContainer->Add(fhMiMod[imod]) ;
919
8d230fa8 920 if(fCalorimeter=="PHOS"){
921 snprintf(key, buffersize,"hMiDiffPHOSMod_%d",imod) ;
922 snprintf(title, buffersize,"Mixed pairs PHOS, clusters in different Modules: %s",(pairnamePHOS[imod]).Data()) ;
923 fhMiDiffPHOSMod[imod] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
924 fhMiDiffPHOSMod[imod]->SetXTitle("p_{T} (GeV/c)");
925 fhMiDiffPHOSMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
926 outputContainer->Add(fhMiDiffPHOSMod[imod]) ;
927 }//PHOS
928 else{//EMCAL
929 if(imod<fNModules/2){
930 snprintf(key, buffersize,"hMiSameSectorEMCALMod_%d",imod) ;
931 snprintf(title, buffersize,"Mixed pairs EMCAL, clusters in same sector, SM(%d,%d)",imod*2,imod*2+1) ;
932 fhMiSameSectorEMCALMod[imod] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
933 fhMiSameSectorEMCALMod[imod]->SetXTitle("p_{T} (GeV/c)");
934 fhMiSameSectorEMCALMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
935 outputContainer->Add(fhMiSameSectorEMCALMod[imod]) ;
936 }
937 if(imod<fNModules-2){
938 snprintf(key, buffersize,"hMiSameSideEMCALMod_%d",imod) ;
939 snprintf(title, buffersize,"Mixed pairs EMCAL, clusters in same side SM(%d,%d)",imod, imod+2) ;
940 fhMiSameSideEMCALMod[imod] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
941 fhMiSameSideEMCALMod[imod]->SetXTitle("p_{T} (GeV/c)");
942 fhMiSameSideEMCALMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
943 outputContainer->Add(fhMiSameSideEMCALMod[imod]) ;
944 }
945 }//EMCAL
946 }
6175da48 947 }
8d230fa8 948
eee5fcf1 949// for(Int_t i = 0; i < outputContainer->GetEntries() ; i++){
950//
951// printf("Histogram %d, name: %s\n ",i, outputContainer->At(i)->GetName());
952//
953// }
954
477d6cee 955 return outputContainer;
1c5acb87 956}
957
958//_________________________________________________________________________________________________________________________________________________
959void AliAnaPi0::Print(const Option_t * /*opt*/) const
960{
477d6cee 961 //Print some relevant parameters set for the analysis
a3aebfff 962 printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
477d6cee 963 AliAnaPartCorrBaseClass::Print(" ");
a3aebfff 964
477d6cee 965 printf("Number of bins in Centrality: %d \n",fNCentrBin) ;
5025c139 966 printf("Number of bins in Z vert. pos: %d \n",GetNZvertBin()) ;
967 printf("Number of bins in Reac. Plain: %d \n",GetNRPBin()) ;
477d6cee 968 printf("Depth of event buffer: %d \n",fNmaxMixEv) ;
af7b3903 969 printf("Pair in same Module: %d \n",fSameSM) ;
477d6cee 970 printf("Cuts: \n") ;
41121cfe 971// printf("Z vertex position: -%2.3f < z < %2.3f \n",GetZvertexCut(),GetZvertexCut()) ; //It crashes here, why?
50f39b97 972 printf("Number of modules: %d \n",fNModules) ;
6175da48 973 printf("Select pairs with their angle: %d, edep %d, min angle %2.3f, max angle %2.3f \n",fUseAngleCut, fUseAngleEDepCut, fAngleCut, fAngleMaxCut) ;
af7b3903 974 printf("Asymmetry cuts: n = %d, \n",fNAsymCuts) ;
975 printf("\tasymmetry < ");
976 for(Int_t i = 0; i < fNAsymCuts; i++) printf("%2.2f ",fAsymCuts[i]);
977 printf("\n");
978
979 printf("PID selection bits: n = %d, \n",fNPIDBits) ;
980 printf("\tPID bit = ");
981 for(Int_t i = 0; i < fNPIDBits; i++) printf("%d ",fPIDBits[i]);
982 printf("\n");
983
db2bf6fd 984 if(fMultiCutAna){
985 printf("pT cuts: n = %d, \n",fNPtCuts) ;
986 printf("\tpT > ");
987 for(Int_t i = 0; i < fNPtCuts; i++) printf("%2.2f ",fPtCuts[i]);
988 printf("GeV/c\n");
989
990 printf("N cell in cluster cuts: n = %d, \n",fNCellNCuts) ;
991 printf("\tnCell > ");
992 for(Int_t i = 0; i < fNCellNCuts; i++) printf("%d ",fCellNCuts[i]);
993 printf("\n");
994
db2bf6fd 995 }
477d6cee 996 printf("------------------------------------------------------\n") ;
1c5acb87 997}
998
5ae09196 999//_____________________________________________________________
1000void AliAnaPi0::FillAcceptanceHistograms(){
1001 //Fill acceptance histograms if MC data is available
c8fe2783 1002
6175da48 1003 if(GetReader()->ReadStack()){
5ae09196 1004 AliStack * stack = GetMCStack();
6175da48 1005 if(stack){
5ae09196 1006 for(Int_t i=0 ; i<stack->GetNprimary(); i++){
1007 TParticle * prim = stack->Particle(i) ;
6175da48 1008 Int_t pdg = prim->GetPdgCode();
1009 if( pdg == 111 || pdg == 221){
5ae09196 1010 Double_t pi0Pt = prim->Pt() ;
1011 //printf("pi0, pt %2.2f\n",pi0Pt);
1012 if(prim->Energy() == TMath::Abs(prim->Pz())) continue ; //Protection against floating point exception
1013 Double_t pi0Y = 0.5*TMath::Log((prim->Energy()-prim->Pz())/(prim->Energy()+prim->Pz())) ;
1014 Double_t phi = TMath::RadToDeg()*prim->Phi() ;
6175da48 1015 if(pdg == 111){
6eb99ccd 1016 if(TMath::Abs(pi0Y) < 1.0){
6175da48 1017 fhPrimPi0Pt->Fill(pi0Pt) ;
1018 }
156549ae 1019 fhPrimPi0Y ->Fill(pi0Pt, pi0Y) ;
1020 fhPrimPi0Phi->Fill(pi0Pt, phi) ;
6175da48 1021 }
1022 else if(pdg == 221){
6eb99ccd 1023 if(TMath::Abs(pi0Y) < 1.0){
6175da48 1024 fhPrimEtaPt->Fill(pi0Pt) ;
1025 }
156549ae 1026 fhPrimEtaY ->Fill(pi0Pt, pi0Y) ;
1027 fhPrimEtaPhi->Fill(pi0Pt, phi) ;
5ae09196 1028 }
08a56f5f 1029
1030 //Origin of meson
1031 Int_t momindex = prim->GetFirstMother();
41121cfe 1032 if(momindex < 0) continue;
08a56f5f 1033 TParticle* mother = stack->Particle(momindex);
1034 Int_t mompdg = TMath::Abs(mother->GetPdgCode());
1035 Int_t momstatus = mother->GetStatusCode();
1036 if(pdg == 111){
1037 if (momstatus == 21)fhPrimPi0PtOrigin->Fill(pi0Pt,0.5);//parton
1038 else if(mompdg < 22 ) fhPrimPi0PtOrigin->Fill(pi0Pt,1.5);//quark
1039 else if(mompdg > 2100 && mompdg < 2210) fhPrimPi0PtOrigin->Fill(pi0Pt,2.5);// resonances
1040 else if(mompdg == 221) fhPrimPi0PtOrigin->Fill(pi0Pt,8.5);//eta
1041 else if(mompdg == 331) fhPrimPi0PtOrigin->Fill(pi0Pt,9.5);//eta prime
1042 else if(mompdg == 213) fhPrimPi0PtOrigin->Fill(pi0Pt,4.5);//rho
1043 else if(mompdg == 223) fhPrimPi0PtOrigin->Fill(pi0Pt,5.5);//omega
1044 else if(mompdg >= 310 && mompdg <= 323) fhPrimPi0PtOrigin->Fill(pi0Pt,6.5);//k0S, k+-,k*
1045 else if(mompdg == 130) fhPrimPi0PtOrigin->Fill(pi0Pt,6.5);//k0L
1046 else if(momstatus == 11 || momstatus == 12 ) fhPrimPi0PtOrigin->Fill(pi0Pt,3.5);//resonances
1047 else fhPrimPi0PtOrigin->Fill(pi0Pt,7.5);//other?
1048 }//pi0
1049 else {
1050 if (momstatus == 21 ) fhPrimEtaPtOrigin->Fill(pi0Pt,0.5);//parton
1051 else if(mompdg < 22 ) fhPrimEtaPtOrigin->Fill(pi0Pt,1.5);//quark
1052 else if(mompdg > 2100 && mompdg < 2210) fhPrimEtaPtOrigin->Fill(pi0Pt,2.5);//qq resonances
1053 else if(mompdg == 331) fhPrimEtaPtOrigin->Fill(pi0Pt,5.5);//eta prime
1054 else if(momstatus == 11 || momstatus == 12 ) fhPrimEtaPtOrigin->Fill(pi0Pt,3.5);//resonances
1055 else fhPrimEtaPtOrigin->Fill(pi0Pt,4.5);//stable, conversions?
1056 //printf("Other Meson pdg %d, Mother %s, pdg %d, status %d\n",pdg, TDatabasePDG::Instance()->GetParticle(mompdg)->GetName(),mompdg, momstatus );
1057 }
1058
1059
5ae09196 1060 //Check if both photons hit Calorimeter
6175da48 1061 if(prim->GetNDaughters()!=2) return; //Only interested in 2 gamma decay
5ae09196 1062 Int_t iphot1=prim->GetFirstDaughter() ;
1063 Int_t iphot2=prim->GetLastDaughter() ;
1064 if(iphot1>-1 && iphot1<stack->GetNtrack() && iphot2>-1 && iphot2<stack->GetNtrack()){
1065 TParticle * phot1 = stack->Particle(iphot1) ;
1066 TParticle * phot2 = stack->Particle(iphot2) ;
1067 if(phot1 && phot2 && phot1->GetPdgCode()==22 && phot2->GetPdgCode()==22){
1068 //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",
1069 // phot1->Pt(), phot1->Phi()*180./3.1415, phot1->Eta(), phot2->Pt(), phot2->Phi()*180./3.1415, phot2->Eta());
1070
1071 TLorentzVector lv1, lv2;
1072 phot1->Momentum(lv1);
1073 phot2->Momentum(lv2);
1074
1075 Bool_t inacceptance = kFALSE;
1076 if(fCalorimeter == "PHOS"){
1077 if(GetPHOSGeometry() && GetCaloUtils()->IsPHOSGeoMatrixSet()){
1078 Int_t mod ;
1079 Double_t x,z ;
1080 if(GetPHOSGeometry()->ImpactOnEmc(phot1,mod,z,x) && GetPHOSGeometry()->ImpactOnEmc(phot2,mod,z,x))
1081 inacceptance = kTRUE;
1082 if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
1083 }
1084 else{
1085
1086 if(GetFiducialCut()->IsInFiducialCut(lv1,fCalorimeter) && GetFiducialCut()->IsInFiducialCut(lv2,fCalorimeter))
1087 inacceptance = kTRUE ;
1088 if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
1089 }
1090
1091 }
1092 else if(fCalorimeter == "EMCAL" && GetCaloUtils()->IsEMCALGeoMatrixSet()){
1093 if(GetEMCALGeometry()){
156549ae 1094
1095 Int_t absID1=0;
1096 Int_t absID2=0;
1097
1098 GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(phot1->Eta(),phot1->Phi(),absID1);
1099 GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(phot2->Eta(),phot2->Phi(),absID2);
1100
1101 if( absID1 >= 0 && absID2 >= 0)
5ae09196 1102 inacceptance = kTRUE;
156549ae 1103
1104// if(GetEMCALGeometry()->Impact(phot1) && GetEMCALGeometry()->Impact(phot2))
1105// inacceptance = kTRUE;
5ae09196 1106 if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
1107 }
1108 else{
1109 if(GetFiducialCut()->IsInFiducialCut(lv1,fCalorimeter) && GetFiducialCut()->IsInFiducialCut(lv2,fCalorimeter))
1110 inacceptance = kTRUE ;
1111 if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
1112 }
1113 }
1114
1115 if(inacceptance){
6175da48 1116 if(pdg==111){
156549ae 1117 fhPrimPi0AccPt ->Fill(pi0Pt) ;
1118 fhPrimPi0AccPhi->Fill(pi0Pt, phi) ;
1119 fhPrimPi0AccY ->Fill(pi0Pt, pi0Y) ;
6175da48 1120 Double_t angle = lv1.Angle(lv2.Vect());
1121 fhPrimPi0OpeningAngle ->Fill(pi0Pt,angle);
1122 fhPrimPi0CosOpeningAngle->Fill(pi0Pt,TMath::Cos(angle));
1123 }
1124 else if(pdg==221){
156549ae 1125 fhPrimEtaAccPt ->Fill(pi0Pt) ;
1126 fhPrimEtaAccPhi->Fill(pi0Pt, phi) ;
1127 fhPrimEtaAccY ->Fill(pi0Pt, pi0Y) ;
6175da48 1128 }
5ae09196 1129 }//Accepted
1130 }// 2 photons
1131 }//Check daughters exist
156549ae 1132 }// Primary pi0 or eta
5ae09196 1133 }//loop on primaries
1134 }//stack exists and data is MC
1135 }//read stack
1136 else if(GetReader()->ReadAODMCParticles()){
6175da48 1137
1138 TClonesArray * mcparticles = GetReader()->GetAODMCParticles(0);
1139 if(mcparticles){
1140 Int_t nprim = mcparticles->GetEntriesFast();
08a56f5f 1141 for(Int_t i=0; i < nprim; i++)
156549ae 1142 {
08a56f5f 1143 AliAODMCParticle * prim = (AliAODMCParticle *) mcparticles->At(i);
1144 // Only generator particles
1145 if( prim->GetStatus() == 0) break;
1146
6175da48 1147 Int_t pdg = prim->GetPdgCode();
1148 if( pdg == 111 || pdg == 221){
1149 Double_t pi0Pt = prim->Pt() ;
1150 //printf("pi0, pt %2.2f\n",pi0Pt);
1151 if(prim->E() == TMath::Abs(prim->Pz())) continue ; //Protection against floating point exception
1152 Double_t pi0Y = 0.5*TMath::Log((prim->E()-prim->Pz())/(prim->E()+prim->Pz())) ;
1153 Double_t phi = TMath::RadToDeg()*prim->Phi() ;
1154 if(pdg == 111){
1155 if(TMath::Abs(pi0Y) < 0.5){
1156 fhPrimPi0Pt->Fill(pi0Pt) ;
1157 }
156549ae 1158 fhPrimPi0Y ->Fill(pi0Pt, pi0Y) ;
1159 fhPrimPi0Phi->Fill(pi0Pt, phi) ;
6175da48 1160 }
1161 else if(pdg == 221){
1162 if(TMath::Abs(pi0Y) < 0.5){
1163 fhPrimEtaPt->Fill(pi0Pt) ;
1164 }
156549ae 1165 fhPrimEtaY ->Fill(pi0Pt, pi0Y) ;
1166 fhPrimEtaPhi->Fill(pi0Pt, phi) ;
6175da48 1167 }
08a56f5f 1168
1169 //Origin of meson
1170 Int_t momindex = prim->GetMother();
41121cfe 1171 if(momindex < 0) continue;
08a56f5f 1172 AliAODMCParticle* mother = (AliAODMCParticle *) mcparticles->At(momindex);
1173 Int_t mompdg = TMath::Abs(mother->GetPdgCode());
1174 Int_t momstatus = mother->GetStatus();
1175 if(pdg == 111){
1176 if (momstatus == 21) fhPrimPi0PtOrigin->Fill(pi0Pt,0.5);//parton
1177 else if(mompdg < 22 ) fhPrimPi0PtOrigin->Fill(pi0Pt,1.5);//quark
1178 else if(mompdg > 2100 && mompdg < 2210) fhPrimPi0PtOrigin->Fill(pi0Pt,2.5);// resonances
1179 else if(mompdg == 221) fhPrimPi0PtOrigin->Fill(pi0Pt,8.5);//eta
1180 else if(mompdg == 331) fhPrimPi0PtOrigin->Fill(pi0Pt,9.5);//eta prime
1181 else if(mompdg == 213) fhPrimPi0PtOrigin->Fill(pi0Pt,4.5);//rho
1182 else if(mompdg == 223) fhPrimPi0PtOrigin->Fill(pi0Pt,5.5);//omega
1183 else if(mompdg >= 310 && mompdg <= 323) fhPrimPi0PtOrigin->Fill(pi0Pt,6.5);//k0S, k+-,k*
1184 else if(mompdg == 130) fhPrimPi0PtOrigin->Fill(pi0Pt,6.5);//k0L
1185 else if(momstatus == 11 || momstatus == 12 ) fhPrimPi0PtOrigin->Fill(pi0Pt,3.5);//resonances
1186 else fhPrimPi0PtOrigin->Fill(pi0Pt,7.5);//other?
1187 }//pi0
1188 else {
1189 if (momstatus == 21 ) fhPrimEtaPtOrigin->Fill(pi0Pt,0.5);//parton
1190 else if(mompdg < 22 ) fhPrimEtaPtOrigin->Fill(pi0Pt,1.5);//quark
1191 else if(mompdg > 2100 && mompdg < 2210) fhPrimEtaPtOrigin->Fill(pi0Pt,2.5);//qq resonances
1192 else if(mompdg == 331) fhPrimEtaPtOrigin->Fill(pi0Pt,5.5);//eta prime
1193 else if(momstatus == 11 || momstatus == 12 ) fhPrimEtaPtOrigin->Fill(pi0Pt,3.5);//resonances
1194 else fhPrimEtaPtOrigin->Fill(pi0Pt,4.5);//stable, conversions?
1195 //printf("Other Meson pdg %d, Mother %s, pdg %d, status %d\n",pdg, TDatabasePDG::Instance()->GetParticle(mompdg)->GetName(),mompdg, momstatus );
1196 }
1197
6175da48 1198 //Check if both photons hit Calorimeter
1199 if(prim->GetNDaughters()!=2) return; //Only interested in 2 gamma decay
1200 Int_t iphot1=prim->GetDaughter(0) ;
1201 Int_t iphot2=prim->GetDaughter(1) ;
1202 if(iphot1>-1 && iphot1<nprim && iphot2>-1 && iphot2<nprim){
1203 AliAODMCParticle * phot1 = (AliAODMCParticle *) mcparticles->At(iphot1);
1204 AliAODMCParticle * phot2 = (AliAODMCParticle *) mcparticles->At(iphot2);
1205 if(phot1 && phot2 && phot1->GetPdgCode()==22 && phot2->GetPdgCode()==22){
6175da48 1206 TLorentzVector lv1, lv2;
1207 lv1.SetPxPyPzE(phot1->Px(),phot1->Py(),phot1->Pz(),phot1->E());
1208 lv2.SetPxPyPzE(phot2->Px(),phot2->Py(),phot2->Pz(),phot2->E());
1209
1210 Bool_t inacceptance = kFALSE;
1211 if(fCalorimeter == "PHOS"){
1212 if(GetPHOSGeometry() && GetCaloUtils()->IsPHOSGeoMatrixSet()){
1213 Int_t mod ;
1214 Double_t x,z ;
1215 Double_t vtx []={phot1->Xv(),phot1->Yv(),phot1->Zv()};
1216 Double_t vtx2[]={phot2->Xv(),phot2->Yv(),phot2->Zv()};
1217 if(GetPHOSGeometry()->ImpactOnEmc(vtx, phot1->Theta(),phot1->Phi(),mod,z,x) &&
1218 GetPHOSGeometry()->ImpactOnEmc(vtx2,phot2->Theta(),phot2->Phi(),mod,z,x))
1219 inacceptance = kTRUE;
1220 if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
1221 }
1222 else{
1223
1224 if(GetFiducialCut()->IsInFiducialCut(lv1,fCalorimeter) && GetFiducialCut()->IsInFiducialCut(lv2,fCalorimeter))
1225 inacceptance = kTRUE ;
1226 if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
1227 }
1228
1229 }
1230 else if(fCalorimeter == "EMCAL" && GetCaloUtils()->IsEMCALGeoMatrixSet()){
1231 if(GetEMCALGeometry()){
156549ae 1232
6175da48 1233 Int_t absID1=0;
6175da48 1234 Int_t absID2=0;
156549ae 1235
1236 //TVector3 vtx(phot1->Xv(),phot1->Yv(),phot1->Zv());
1237 //TVector3 vimpact(0,0,0);
1238
1239 //GetEMCALGeometry()->ImpactOnEmcal(vtx,phot1->Theta(),phot1->Phi(),absID1,vimpact);
1240 //TVector3 vtx2(phot2->Xv(),phot2->Yv(),phot2->Zv());
1241 //TVector3 vimpact2(0,0,0);
1242 //GetEMCALGeometry()->ImpactOnEmcal(vtx2,phot2->Theta(),phot2->Phi(),absID2,vimpact2);
1243
1244 GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(phot1->Eta(),phot1->Phi(),absID1);
1245 GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(phot2->Eta(),phot2->Phi(),absID2);
1246
6175da48 1247// if(TMath::Abs(phot1->Eta()) < 0.7 && phot1->Phi() > 80*TMath::DegToRad() && phot1->Phi() < 120*TMath::DegToRad() )
156549ae 1248// printf("Phot1 ccepted? %d\n",absID1);
6175da48 1249// if(TMath::Abs(phot2->Eta()) < 0.7 && phot2->Phi() > 80*TMath::DegToRad() && phot2->Phi() < 120*TMath::DegToRad() )
1250// printf("Phot2 accepted? %d\n",absID2);
1251
1252 if( absID1 >= 0 && absID2 >= 0)
1253 inacceptance = kTRUE;
156549ae 1254
1255// 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",
1256// absID1,phot1->Pt(), phot1->Phi()*TMath::RadToDeg(), phot1->Eta(),
1257// absID2,phot2->Pt(), phot2->Phi()*TMath::RadToDeg(), phot2->Eta());
1258
1259
1260
6175da48 1261 if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
1262 }
1263 else{
1264 if(GetFiducialCut()->IsInFiducialCut(lv1,fCalorimeter) && GetFiducialCut()->IsInFiducialCut(lv2,fCalorimeter))
1265 inacceptance = kTRUE ;
1266 if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
1267 }
1268 }
1269
1270 if(inacceptance){
1271 if(pdg==111){
156549ae 1272 // printf("ACCEPTED pi0: pt %2.2f, phi %3.2f, eta %1.2f\n",pi0Pt,phi,pi0Y);
1273 fhPrimPi0AccPt ->Fill(pi0Pt) ;
1274 fhPrimPi0AccPhi->Fill(pi0Pt, phi) ;
1275 fhPrimPi0AccY ->Fill(pi0Pt, pi0Y) ;
6175da48 1276 Double_t angle = lv1.Angle(lv2.Vect());
1277 fhPrimPi0OpeningAngle ->Fill(pi0Pt,angle);
1278 fhPrimPi0CosOpeningAngle->Fill(pi0Pt,TMath::Cos(angle));
1279 }
1280 else if(pdg==221){
156549ae 1281 fhPrimEtaAccPt ->Fill(pi0Pt) ;
1282 fhPrimEtaAccPhi->Fill(pi0Pt, phi) ;
1283 fhPrimEtaAccY ->Fill(pi0Pt, pi0Y) ;
6175da48 1284 }
1285 }//Accepted
1286 }// 2 photons
1287 }//Check daughters exist
156549ae 1288 }// Primary pi0 or eta
6175da48 1289 }//loop on primaries
1290 }//stack exists and data is MC
1291
1292
1293 } // read AOD MC
5ae09196 1294}
1295
6175da48 1296//_____________________________________________________________
1297void AliAnaPi0::FillMCVersusRecDataHistograms(const Int_t index1, const Int_t index2,
1298 const Float_t pt1, const Float_t pt2,
1299 const Int_t ncell1, const Int_t ncell2,
1300 const Double_t mass, const Double_t pt, const Double_t asym,
1301 const Double_t deta, const Double_t dphi){
1302 //Do some MC checks on the origin of the pair, is there any common ancestor and if there is one, who?
1303 //Adjusted for Pythia, need to see what to do for other generators.
1304 //Array of histograms ordered as follows: 0-Photon, 1-electron, 2-pi0, 3-eta, 4-a-proton, 5-a-neutron, 6-stable particles,
1305 // 7-other decays, 8-string, 9-final parton, 10-initial parton, intermediate, 11-colliding proton, 12-unrelated
1306
1307 Int_t ancPDG = 0;
1308 Int_t ancStatus = 0;
1309 TLorentzVector ancMomentum;
1310 Int_t ancLabel = GetMCAnalysisUtils()->CheckCommonAncestor(index1, index2,
1311 GetReader(), ancPDG, ancStatus,ancMomentum);
1312
08a56f5f 1313 Int_t momindex = -1;
1314 Int_t mompdg = -1;
1315 Int_t momstatus = -1;
6175da48 1316 if(GetDebug() > 1) printf("AliAnaPi0::FillMCVersusRecDataHistograms() - Common ancestor label %d, pdg %d, name %s, status %d; \n",
1317 ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus);
1318
1319 if(ancLabel > -1){
1320 if(ancPDG==22){//gamma
1321 fhMCOrgMass[0]->Fill(pt,mass);
1322 fhMCOrgAsym[0]->Fill(pt,asym);
1323 fhMCOrgDeltaEta[0]->Fill(pt,deta);
1324 fhMCOrgDeltaPhi[0]->Fill(pt,dphi);
1325 }
1326 else if(TMath::Abs(ancPDG)==11){//e
1327 fhMCOrgMass[1]->Fill(pt,mass);
1328 fhMCOrgAsym[1]->Fill(pt,asym);
1329 fhMCOrgDeltaEta[1]->Fill(pt,deta);
1330 fhMCOrgDeltaPhi[1]->Fill(pt,dphi);
1331 }
1332 else if(ancPDG==111){//Pi0
1333 fhMCOrgMass[2]->Fill(pt,mass);
1334 fhMCOrgAsym[2]->Fill(pt,asym);
1335 fhMCOrgDeltaEta[2]->Fill(pt,deta);
1336 fhMCOrgDeltaPhi[2]->Fill(pt,dphi);
1337 if(fMultiCutAnaSim){
1338 for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
1339 for(Int_t icell=0; icell<fNCellNCuts; icell++){
1340 for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
1341 Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
1342 if(pt1 > fPtCuts[ipt] && pt2 > fPtCuts[ipt] &&
1343 asym < fAsymCuts[iasym] &&
1344 ncell1 >= fCellNCuts[icell] && ncell2 >= fCellNCuts[icell]){
1345 fhMCPi0MassPtRec [index]->Fill(pt,mass);
1346 fhMCPi0MassPtTrue[index]->Fill(ancMomentum.Pt(),mass);
1347 if(mass < 0.17 && mass > 0.1) fhMCPi0PtTruePtRec[index]->Fill(ancMomentum.Pt(),pt);
1348 }//pass the different cuts
1349 }// pid bit cut loop
1350 }// icell loop
1351 }// pt cut loop
1352 }//Multi cut ana sim
1353 else {
1354 fhMCPi0MassPtTrue[0]->Fill(ancMomentum.Pt(),mass);
08a56f5f 1355 if(mass < 0.17 && mass > 0.1) {
1356 fhMCPi0PtTruePtRec[0]->Fill(ancMomentum.Pt(),pt);
1357
1358 if(GetReader()->ReadStack()){
1359 TParticle* ancestor = GetMCStack()->Particle(ancLabel);
1360 momindex = ancestor->GetFirstMother();
41121cfe 1361 if(momindex < 0) return;
08a56f5f 1362 TParticle* mother = GetMCStack()->Particle(momindex);
1363 mompdg = TMath::Abs(mother->GetPdgCode());
1364 momstatus = mother->GetStatusCode();
1365 }
1366 else {
1367 TClonesArray * mcparticles = GetReader()->GetAODMCParticles(0);
1368 AliAODMCParticle* ancestor = (AliAODMCParticle *) mcparticles->At(ancLabel);
1369 momindex = ancestor->GetMother();
41121cfe 1370 if(momindex < 0) return;
08a56f5f 1371 AliAODMCParticle* mother = (AliAODMCParticle *) mcparticles->At(momindex);
1372 mompdg = TMath::Abs(mother->GetPdgCode());
1373 momstatus = mother->GetStatus();
1374 }
1375
1376 if (momstatus == 21) fhMCPi0PtOrigin->Fill(pt,0.5);//parton
1377 else if(mompdg < 22 ) fhMCPi0PtOrigin->Fill(pt,1.5);//quark
1378 else if(mompdg > 2100 && mompdg < 2210) fhMCPi0PtOrigin->Fill(pt,2.5);// resonances
1379 else if(mompdg == 221) fhMCPi0PtOrigin->Fill(pt,8.5);//eta
1380 else if(mompdg == 331) fhMCPi0PtOrigin->Fill(pt,9.5);//eta prime
1381 else if(mompdg == 213) fhMCPi0PtOrigin->Fill(pt,4.5);//rho
1382 else if(mompdg == 223) fhMCPi0PtOrigin->Fill(pt,5.5);//omega
1383 else if(mompdg >= 310 && mompdg <= 323) fhMCPi0PtOrigin->Fill(pt,6.5);//k0S, k+-,k*
1384 else if(mompdg == 130) fhMCPi0PtOrigin->Fill(pt,6.5);//k0L
1385 else if(momstatus == 11 || momstatus == 12 ) fhMCPi0PtOrigin->Fill(pt,3.5);//resonances
1386 else fhMCPi0PtOrigin->Fill(pt,7.5);//other?
1387
1388 }//pi0 mass region
1389
6175da48 1390 }
1391 }
1392 else if(ancPDG==221){//Eta
1393 fhMCOrgMass[3]->Fill(pt,mass);
1394 fhMCOrgAsym[3]->Fill(pt,asym);
1395 fhMCOrgDeltaEta[3]->Fill(pt,deta);
1396 fhMCOrgDeltaPhi[3]->Fill(pt,dphi);
1397 if(fMultiCutAnaSim){
1398 for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
1399 for(Int_t icell=0; icell<fNCellNCuts; icell++){
1400 for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
1401 Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
1402 if(pt1 > fPtCuts[ipt] && pt2 > fPtCuts[ipt] &&
1403 asym < fAsymCuts[iasym] &&
1404 ncell1 >= fCellNCuts[icell] && ncell2 >= fCellNCuts[icell]){
1405 fhMCEtaMassPtRec [index]->Fill(pt,mass);
1406 fhMCEtaMassPtTrue[index]->Fill(ancMomentum.Pt(),mass);
156549ae 1407 if(mass < 0.65 && mass > 0.45) fhMCEtaPtTruePtRec[index]->Fill(ancMomentum.Pt(),pt);
6175da48 1408 }//pass the different cuts
1409 }// pid bit cut loop
1410 }// icell loop
1411 }// pt cut loop
1412 } //Multi cut ana sim
1413 else {
1414 fhMCEtaMassPtTrue[0]->Fill(ancMomentum.Pt(),mass);
156549ae 1415 if(mass < 0.65 && mass > 0.45) fhMCEtaPtTruePtRec[0]->Fill(ancMomentum.Pt(),pt);
08a56f5f 1416
1417 if(GetReader()->ReadStack()){
1418 TParticle* ancestor = GetMCStack()->Particle(ancLabel);
1419 momindex = ancestor->GetFirstMother();
41121cfe 1420 if(momindex < 0) return;
08a56f5f 1421 TParticle* mother = GetMCStack()->Particle(momindex);
1422 mompdg = TMath::Abs(mother->GetPdgCode());
1423 momstatus = mother->GetStatusCode();
1424 }
1425 else {
1426 TClonesArray * mcparticles = GetReader()->GetAODMCParticles(0);
1427 AliAODMCParticle* ancestor = (AliAODMCParticle *) mcparticles->At(ancLabel);
1428 momindex = ancestor->GetMother();
41121cfe 1429 if(momindex < 0) return;
08a56f5f 1430 AliAODMCParticle* mother = (AliAODMCParticle *) mcparticles->At(momindex);
1431 mompdg = TMath::Abs(mother->GetPdgCode());
1432 momstatus = mother->GetStatus();
1433 }
1434
1435 if (momstatus == 21 ) fhMCEtaPtOrigin->Fill(pt,0.5);//parton
1436 else if(mompdg < 22 ) fhMCEtaPtOrigin->Fill(pt,1.5);//quark
1437 else if(mompdg > 2100 && mompdg < 2210) fhMCEtaPtOrigin->Fill(pt,2.5);//qq resonances
1438 else if(mompdg == 331) fhMCEtaPtOrigin->Fill(pt,5.5);//eta prime
1439 else if(momstatus == 11 || momstatus == 12 ) fhMCEtaPtOrigin->Fill(pt,3.5);//resonances
1440 else fhMCEtaPtOrigin->Fill(pt,4.5);//stable, conversions?
1441 //printf("Other Meson pdg %d, Mother %s, pdg %d, status %d\n",pdg, TDatabasePDG::Instance()->GetParticle(mompdg)->GetName(),mompdg, momstatus );
1442 }// eta mass region
6175da48 1443 }
1444 else if(ancPDG==-2212){//AProton
1445 fhMCOrgMass[4]->Fill(pt,mass);
1446 fhMCOrgAsym[4]->Fill(pt,asym);
1447 fhMCOrgDeltaEta[4]->Fill(pt,deta);
1448 fhMCOrgDeltaPhi[4]->Fill(pt,dphi);
1449 }
1450 else if(ancPDG==-2112){//ANeutron
1451 fhMCOrgMass[5]->Fill(pt,mass);
1452 fhMCOrgAsym[5]->Fill(pt,asym);
1453 fhMCOrgDeltaEta[5]->Fill(pt,deta);
1454 fhMCOrgDeltaPhi[5]->Fill(pt,dphi);
1455 }
1456 else if(TMath::Abs(ancPDG)==13){//muons
1457 fhMCOrgMass[6]->Fill(pt,mass);
1458 fhMCOrgAsym[6]->Fill(pt,asym);
1459 fhMCOrgDeltaEta[6]->Fill(pt,deta);
1460 fhMCOrgDeltaPhi[6]->Fill(pt,dphi);
1461 }
1462 else if (TMath::Abs(ancPDG) > 100 && ancLabel > 7) {
1463 if(ancStatus==1){//Stable particles, converted? not decayed resonances
1464 fhMCOrgMass[6]->Fill(pt,mass);
1465 fhMCOrgAsym[6]->Fill(pt,asym);
1466 fhMCOrgDeltaEta[6]->Fill(pt,deta);
1467 fhMCOrgDeltaPhi[6]->Fill(pt,dphi);
1468 }
1469 else{//resonances and other decays, more hadron conversions?
1470 fhMCOrgMass[7]->Fill(pt,mass);
1471 fhMCOrgAsym[7]->Fill(pt,asym);
1472 fhMCOrgDeltaEta[7]->Fill(pt,deta);
1473 fhMCOrgDeltaPhi[7]->Fill(pt,dphi);
1474 }
1475 }
1476 else {//Partons, colliding protons, strings, intermediate corrections
1477 if(ancStatus==11 || ancStatus==12){//String fragmentation
1478 fhMCOrgMass[8]->Fill(pt,mass);
1479 fhMCOrgAsym[8]->Fill(pt,asym);
1480 fhMCOrgDeltaEta[8]->Fill(pt,deta);
1481 fhMCOrgDeltaPhi[8]->Fill(pt,dphi);
1482 }
1483 else if (ancStatus==21){
1484 if(ancLabel < 2) {//Colliding protons
1485 fhMCOrgMass[11]->Fill(pt,mass);
1486 fhMCOrgAsym[11]->Fill(pt,asym);
1487 fhMCOrgDeltaEta[11]->Fill(pt,deta);
1488 fhMCOrgDeltaPhi[11]->Fill(pt,dphi);
1489 }//colliding protons
1490 else if(ancLabel < 6){//partonic initial states interactions
1491 fhMCOrgMass[9]->Fill(pt,mass);
1492 fhMCOrgAsym[9]->Fill(pt,asym);
1493 fhMCOrgDeltaEta[9]->Fill(pt,deta);
1494 fhMCOrgDeltaPhi[9]->Fill(pt,dphi);
1495 }
1496 else if(ancLabel < 8){//Final state partons radiations?
1497 fhMCOrgMass[10]->Fill(pt,mass);
1498 fhMCOrgAsym[10]->Fill(pt,asym);
1499 fhMCOrgDeltaEta[10]->Fill(pt,deta);
1500 fhMCOrgDeltaPhi[10]->Fill(pt,dphi);
1501 }
1502 else {
1503 printf("AliAnaPi0::FillMCVersusRecDataHistograms() - Check ** Common ancestor label %d, pdg %d, name %s, status %d; \n",
1504 ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus);
1505 }
1506 }//status 21
1507 else {
1508 printf("AliAnaPi0::FillMCVersusRecDataHistograms() - Check *** Common ancestor label %d, pdg %d, name %s, status %d; \n",
1509 ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus);
1510 }
1511 }////Partons, colliding protons, strings, intermediate corrections
1512 }//ancLabel > -1
1513 else { //ancLabel <= -1
1514 //printf("Not related at all label = %d\n",ancLabel);
1515 fhMCOrgMass[12]->Fill(pt,mass);
1516 fhMCOrgAsym[12]->Fill(pt,asym);
1517 fhMCOrgDeltaEta[12]->Fill(pt,deta);
1518 fhMCOrgDeltaPhi[12]->Fill(pt,dphi);
1519 }
1520}
1521
1c5acb87 1522//____________________________________________________________________________________________________________________________________________________
6639984f 1523void AliAnaPi0::MakeAnalysisFillHistograms()
1c5acb87 1524{
477d6cee 1525 //Process one event and extract photons from AOD branch
1526 // filled with AliAnaPhoton and fill histos with invariant mass
1527
6175da48 1528 //In case of simulated data, fill acceptance histograms
1529 if(IsDataMC())FillAcceptanceHistograms();
08a56f5f 1530
1531 //if (GetReader()->GetEventNumber()%10000 == 0)
1532 // printf("--- Event %d ---\n",GetReader()->GetEventNumber());
1533
6175da48 1534 //Init some variables
1535//Int_t iRun = (GetReader()->GetInputEvent())->GetRunNumber() ;
1536 Int_t nPhot = GetInputAODBranch()->GetEntriesFast() ;
1537 Int_t nClus = 0;
1538 Int_t nCell = 0;
1539 Float_t eClusTot = 0;
1540 Float_t eCellTot = 0;
156549ae 1541 Float_t eDenClus = 0;
1542 Float_t eDenCell = 0;
1543// Int_t ncomb = 0;
1544// Float_t rtmp = 0;
1545// Float_t rtmpw = 0;
1546// Float_t rxz = 0;
1547// Float_t rxzw = 0;
1548// Float_t pos1[3];
1549// Float_t pos2[3];
1550// Float_t emax = 0;
477d6cee 1551
156549ae 1552 if(GetDebug() > 1)
1553 printf("AliAnaPi0::MakeAnalysisFillHistograms() - Photon entries %d\n", nPhot);
1554
1555 //If less than photon 2 entries in the list, skip this event
1556 if(nPhot < 2 ) return ;
1557
6175da48 1558 // Count the number of clusters and cells, in case multiplicity bins dependent on such numbers
1559 // are requested
1560 if(fCalorimeter=="EMCAL"){
be518ab0 1561 nClus = GetEMCALClusters() ->GetEntriesFast();
6175da48 1562 nCell = GetEMCALCells()->GetNumberOfCells();
156549ae 1563 for(Int_t icl=0; icl < nClus; icl++) {
be518ab0 1564 Float_t e1 = ((AliVCluster*)GetEMCALClusters()->At(icl))->E();
156549ae 1565 eClusTot += e1;
1566// if(e1 > emax) emax = e1;
be518ab0 1567// ((AliVCluster*)GetEMCALClusters()->At(icl))->GetPosition(pos1);
156549ae 1568// for(Int_t icl2=icl+1; icl2 < nClus; icl2++) {
be518ab0 1569// Float_t e2 = ((AliVCluster*)GetEMCALClusters()->At(icl2))->E();
1570// ((AliVCluster*)GetEMCALClusters()->At(icl2))->GetPosition(pos2);
156549ae 1571// rtmp = TMath::Sqrt((pos1[0]-pos2[0])*(pos1[0]-pos2[0]) + (pos1[2]-pos2[2])*(pos1[2]-pos2[2]));
1572// 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);
1573// rxz += rtmp;
1574// rxzw += rtmpw;
1575// ncomb++;
1576// fhClusterPairDist ->Fill(rtmp);
1577// fhClusterPairDistWeight->Fill(rtmpw);
be518ab0 1578// //printf("Distance: %f; weighted %f\n ",rtmp,rtmp/(e1+((AliVCluster*)GetEMCALClusters()->At(icl2))->E()));
156549ae 1579//
1580// }// second cluster loop
1581 }// first cluster
1582
6175da48 1583 for(Int_t jce=0; jce < nCell; jce++) eCellTot += GetEMCALCells()->GetAmplitude(jce);
6175da48 1584 }
1585 else {
be518ab0 1586 nClus = GetPHOSClusters()->GetEntriesFast();
1587 nCell = GetPHOSCells() ->GetNumberOfCells();
156549ae 1588 for(Int_t icl=0; icl < nClus; icl++) {
be518ab0 1589 Float_t e1 = ((AliVCluster*)GetPHOSClusters()->At(icl))->E();
156549ae 1590 eClusTot += e1;
be518ab0 1591// ((AliVCluster*)GetPHOSClusters()->At(icl))->GetPosition(pos1);
156549ae 1592// for(Int_t icl2=icl+1; icl2 < nClus; icl2++) {
be518ab0 1593// Float_t e2 = ((AliVCluster*)GetPHOSClusters()->At(icl2))->E();
1594// ((AliVCluster*)GetPHOSClusters()->At(icl2))->GetPosition(pos2);
156549ae 1595// rtmp = TMath::Sqrt((pos1[0]-pos2[0])*(pos1[0]-pos2[0]) + (pos1[2]-pos2[2])*(pos1[2]-pos2[2]));
1596// 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);
1597// rxz += rtmp;
1598// rxzw += rtmpw;
1599// ncomb++;
1600// fhClusterPairDist ->Fill(rtmp);
1601// fhClusterPairDistWeight->Fill(rtmpw);
1602// }// second cluster loop
1603 }// first cluster
6175da48 1604 for(Int_t jce=0; jce < nCell; jce++) eCellTot += GetPHOSCells()->GetAmplitude(jce);
1605 }
156549ae 1606 if(GetDebug() > 1)
1607 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 1608
156549ae 1609 //Fill histograms with "energy density", ncell and nclust will be > 0 since there are at least 2 "photons"
1610 eDenClus = eClusTot/nClus;
1611 eDenCell = eCellTot/nCell;
1612 fhEDensityCluster ->Fill(eDenClus);
1613 fhEDensityCell ->Fill(eDenCell);
1614 fhEDensityCellvsCluster->Fill(eDenClus, eDenCell);
6175da48 1615 //Fill the average number of cells or clusters per SM
1616 eClusTot /=fNModules;
1617 eCellTot /=fNModules;
156549ae 1618 fhAverTotECluster ->Fill(eClusTot);
1619 fhAverTotECell ->Fill(eCellTot);
1620 fhAverTotECellvsCluster->Fill(eClusTot, eCellTot);
1621 //printf("Average Cluster: E %f, density %f; Average Cell E %f, density %f\n ",eClusTot,eDenClus,eCellTot,eDenCell);
6175da48 1622
156549ae 1623// //Average weighted pair distance
1624// rxz /= ncomb;
1625// rxzw /= ncomb;
1626//
1627// fhAverClusterPairDist ->Fill(rxz );
1628// fhAverClusterPairDistWeight ->Fill(rxzw);
1629// fhAverClusterPairDistvsAverE ->Fill(rxz ,eDenClus);
1630// fhAverClusterPairDistWeightvsAverE->Fill(rxzw,eDenClus);
1631// fhAverClusterPairDistvsN ->Fill(rxz ,nClus);
1632// fhAverClusterPairDistWeightvsN ->Fill(rxzw,nClus);
1633//
1634// //emax
1635// fhMaxEvsClustEDen->Fill(emax,eDenClus);
1636// fhMaxEvsClustMult->Fill(emax,nPhot);
7e7694bb 1637
156549ae 1638 //printf("Average Distance: %f; weighted %f\n ",rxz,rxzw);
1639
6175da48 1640
1641 //Init variables
1642 Int_t module1 = -1;
1643 Int_t module2 = -1;
1644 Double_t vert[] = {0.0, 0.0, 0.0} ; //vertex
1645 Int_t evtIndex1 = 0 ;
1646 Int_t currentEvtIndex = -1;
1647 Int_t curCentrBin = 0 ;
1648 Int_t curRPBin = 0 ;
1649 Int_t curZvertBin = 0 ;
1650
1651 //---------------------------------
1652 //First loop on photons/clusters
1653 //---------------------------------
477d6cee 1654 for(Int_t i1=0; i1<nPhot-1; i1++){
1655 AliAODPWG4Particle * p1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i1)) ;
6175da48 1656 //printf("AliAnaPi0::MakeAnalysisFillHistograms() : cluster1 id %d\n",p1->GetCaloLabel(0));
1657
7e7694bb 1658 // get the event index in the mixed buffer where the photon comes from
1659 // in case of mixing with analysis frame, not own mixing
c8fe2783 1660 evtIndex1 = GetEventIndex(p1, vert) ;
5025c139 1661 //printf("charge = %d\n", track->Charge());
c8fe2783 1662 if ( evtIndex1 == -1 )
1663 return ;
1664 if ( evtIndex1 == -2 )
1665 continue ;
41121cfe 1666
1667 //printf("z vertex %f < %f\n",vert[2],GetZvertexCut());
2244659d 1668 if(TMath::Abs(vert[2]) > GetZvertexCut()) continue ; //vertex cut
41121cfe 1669
6175da48 1670
1671 //----------------------------------------------------------------------------
1672 // Get the multiplicity bin. Different cases: centrality (PbPb),
1673 // average cluster multiplicity, average cell multiplicity, track multiplicity
1674 // default is centrality bins
1675 //----------------------------------------------------------------------------
c8fe2783 1676 if (evtIndex1 != currentEvtIndex) {
6175da48 1677 if(fUseTrackMultBins){ // Track multiplicity bins
1678 //printf("track mult %d\n",GetTrackMultiplicity());
1679 curCentrBin = (GetTrackMultiplicity()-1)/5;
1680 if(curCentrBin > fNCentrBin-1) curCentrBin=fNCentrBin-1;
1681 //printf("track mult bin %d\n",curCentrBin);
1682 }
1683 else if(fUsePhotonMultBins){ // Photon multiplicity bins
1684 //printf("photon mult %d cluster mult %d\n",nPhot, nClus);
156549ae 1685 curRPBin = nPhot-2;
1686 if(curRPBin > GetNRPBin() -1) curRPBin=GetNRPBin()-1;
1687 //printf("photon mult bin %d\n",curRPBin);
6175da48 1688 }
156549ae 1689 else if(fUseAverClusterEBins){ // Cluster average energy bins
6175da48 1690 //Bins for pp, if needed can be done in a more general way
08a56f5f 1691 curCentrBin = (Int_t) eClusTot/10 * fNCentrBin;
6175da48 1692 if(curCentrBin > fNCentrBin-1) curCentrBin=fNCentrBin-1;
1693 //printf("cluster E average %f, bin %d \n",eClusTot,curCentrBin);
1694 }
156549ae 1695 else if(fUseAverCellEBins){ // Cell average energy bins
6175da48 1696 //Bins for pp, if needed can be done in a more general way
08a56f5f 1697 curCentrBin = (Int_t) eCellTot/10*fNCentrBin;
6175da48 1698 if(curCentrBin > fNCentrBin-1) curCentrBin=fNCentrBin-1;
1699 //printf("cell E average %f, bin %d \n",eCellTot,curCentrBin);
1700 }
156549ae 1701 else if(fUseAverClusterEDenBins){ // Energy density bins
1702 //Bins for pp, if needed can be done in a more general way
08a56f5f 1703 curCentrBin = (Int_t) eDenClus/10*fNCentrBin;
156549ae 1704 if(curCentrBin > fNCentrBin-1) curCentrBin=fNCentrBin-1;
1705 //printf("cluster Eden average %f, bin %d \n",eDenClus,curCentrBin);
1706 }
1707// else if(fUseAverClusterPairRBins){ // Cluster average distance bins
1708// //Bins for pp, if needed can be done in a more general way
1709// curCentrBin = rxz/650*fNCentrBin;
1710// if(curCentrBin > fNCentrBin-1) curCentrBin=fNCentrBin-1;
1711// //printf("cluster pair R average %f, bin %d \n",rxz,curCentrBin);
1712// }
1713// else if(fUseAverClusterPairRWeightBins){ // Cluster average distance bins
1714// //Bins for pp, if needed can be done in a more general way
1715// curCentrBin = rxzw/350*fNCentrBin;
1716// if(curCentrBin > fNCentrBin-1) curCentrBin=fNCentrBin-1;
1717// //printf("cluster pair rW average %f, bin %d \n",rxzw,curCentrBin);
1718// }
1719// else if(fUseEMaxBins){ // Cluster average distance bins
1720// //Bins for pp, if needed can be done in a more general way
1721// curCentrBin = emax/20*fNCentrBin;
1722// if(curCentrBin > fNCentrBin-1) curCentrBin=fNCentrBin-1;
1723// //printf("cluster pair rW average %f, bin %d \n",rxzw,curCentrBin);
1724// }
6175da48 1725 else { //Event centrality
41121cfe 1726 curCentrBin = GetEventCentrality();
6eb99ccd 1727 //printf("curCentrBin %d\n",curCentrBin);
6175da48 1728 }
6eb99ccd 1729
1730 if (curCentrBin < 0 || curCentrBin >= fNCentrBin){
aba0ecc2 1731 if(GetDebug() > 0)
1732 printf("AliAnaPi0::MakeAnalysisFillHistograms() - Centrality bin <%d> not expected, n bins <%d> , return\n",curCentrBin,fNCentrBin);
6eb99ccd 1733 return;
1734 }
1735
156549ae 1736 //Reaction plane bin
c8fe2783 1737 curRPBin = 0 ;
156549ae 1738
1739 //Get vertex z bin
5025c139 1740 curZvertBin = (Int_t)(0.5*GetNZvertBin()*(vert[2]+GetZvertexCut())/GetZvertexCut()) ;
6175da48 1741
1742 //Fill event bin info
745f04da 1743 fhEvents ->Fill(curCentrBin+0.5,curZvertBin+0.5,curRPBin+0.5) ;
1744 fhCentrality->Fill(curCentrBin);
c8fe2783 1745 currentEvtIndex = evtIndex1 ;
ca468d44 1746 if(GetDebug() > 1)
6175da48 1747 printf("AliAnaPi0::MakeAnalysisFillHistograms() - Centrality %d, Vertex Bin %d, RP bin %d \n",curCentrBin,curRPBin,curZvertBin);
c8fe2783 1748 }
7e7694bb 1749
f8006433 1750 //printf("AliAnaPi0::MakeAnalysisFillHistograms(): Photon 1 Evt %d Vertex : %f,%f,%f\n",evtIndex1, GetVertex(evtIndex1)[0] ,GetVertex(evtIndex1)[1],GetVertex(evtIndex1)[2]);
af7b3903 1751
6175da48 1752 //Get the momentum of this cluster
477d6cee 1753 TLorentzVector photon1(p1->Px(),p1->Py(),p1->Pz(),p1->E());
6175da48 1754
1755 //Get (Super)Module number of this cluster
59b6bd99 1756 module1 = GetModuleNumber(p1);
6175da48 1757
1758 //---------------------------------
1759 //Second loop on photons/clusters
1760 //---------------------------------
477d6cee 1761 for(Int_t i2=i1+1; i2<nPhot; i2++){
1762 AliAODPWG4Particle * p2 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i2)) ;
6175da48 1763
1764 //In case of mixing frame, check we are not in the same event as the first cluster
c8fe2783 1765 Int_t evtIndex2 = GetEventIndex(p2, vert) ;
1766 if ( evtIndex2 == -1 )
1767 return ;
1768 if ( evtIndex2 == -2 )
1769 continue ;
1770 if (GetMixedEvent() && (evtIndex1 == evtIndex2))
7e7694bb 1771 continue ;
6175da48 1772
f8006433 1773 //printf("AliAnaPi0::MakeAnalysisFillHistograms(): Photon 2 Evt %d Vertex : %f,%f,%f\n",evtIndex2, GetVertex(evtIndex2)[0] ,GetVertex(evtIndex2)[1],GetVertex(evtIndex2)[2]);
6175da48 1774
1775 //Get the momentum of this cluster
477d6cee 1776 TLorentzVector photon2(p2->Px(),p2->Py(),p2->Pz(),p2->E());
59b6bd99 1777 //Get module number
6175da48 1778 module2 = GetModuleNumber(p2);
1779
1780 //---------------------------------
1781 // Get pair kinematics
1782 //---------------------------------
1783 Double_t m = (photon1 + photon2).M() ;
1784 Double_t pt = (photon1 + photon2).Pt();
1785 Double_t deta = photon1.Eta() - photon2.Eta();
1786 Double_t dphi = photon1.Phi() - photon2.Phi();
1787 Double_t a = TMath::Abs(p1->E()-p2->E())/(p1->E()+p2->E()) ;
1788
477d6cee 1789 if(GetDebug() > 2)
6175da48 1790 printf(" E: photon1 %f, photon2 %f; Pair: pT %f, mass %f, a %f\n", p1->E(), p2->E(), (photon1 + photon2).E(),m,a);
1791
1792 //--------------------------------
1793 // Opening angle selection
1794 //--------------------------------
50f39b97 1795 //Check if opening angle is too large or too small compared to what is expected
1796 Double_t angle = photon1.Angle(photon2.Vect());
6175da48 1797 if(fUseAngleEDepCut && !GetNeutralMesonSelection()->IsAngleInWindow((photon1+photon2).E(),angle+0.05)) {
1798 if(GetDebug() > 2)
1799 printf("AliAnaPi0::MakeAnalysisFillHistograms() -Real pair angle %f not in E %f window\n",angle, (photon1+photon2).E());
c8fe2783 1800 continue;
6175da48 1801 }
af7b3903 1802
6175da48 1803 if(fUseAngleCut && (angle < fAngleCut || angle > fAngleMaxCut)) {
1804 if(GetDebug() > 2)
1805 printf("AliAnaPi0::MakeAnalysisFillHistograms() - Real pair cut %f < angle %f < cut %f\n",fAngleCut, angle, fAngleMaxCut);
1806 continue;
1807 }
1808
1809 //-------------------------------------------------------------------------------------------------
af7b3903 1810 //Fill module dependent histograms, put a cut on assymmetry on the first available cut in the array
6175da48 1811 //-------------------------------------------------------------------------------------------------
af7b3903 1812 if(a < fAsymCuts[0]){
1813 if(module1==module2 && module1 >=0 && module1<fNModules)
1814 fhReMod[module1]->Fill(pt,m) ;
8d230fa8 1815
af7b3903 1816 if(fCalorimeter=="EMCAL"){
8d230fa8 1817
1818 // Same sector
1819 Int_t j=0;
1820 for(Int_t i = 0; i < fNModules/2; i++){
1821 j=2*i;
1822 if((module1==j && module2==j+1) || (module1==j+1 && module2==j)) fhReSameSectorEMCALMod[i]->Fill(pt,m) ;
1823 }
1824
1825 // Same side
1826 for(Int_t i = 0; i < fNModules-2; i++){
1827 if((module1==i && module2==i+2) || (module1==i+2 && module2==i)) fhReSameSideEMCALMod[i]->Fill(pt,m);
1828 }
1829 }//EMCAL
1830 else {//PHOS
1831 if((module1==0 && module2==1) || (module1==1 && module2==0)) fhReDiffPHOSMod[0]->Fill(pt,m) ;
1832 if((module1==0 && module2==2) || (module1==2 && module2==0)) fhReDiffPHOSMod[1]->Fill(pt,m) ;
1833 if((module1==1 && module2==2) || (module1==2 && module2==1)) fhReDiffPHOSMod[2]->Fill(pt,m) ;
1834 }//PHOS
821c8090 1835 }
7e7694bb 1836
af7b3903 1837 //In case we want only pairs in same (super) module, check their origin.
1838 Bool_t ok = kTRUE;
1839 if(fSameSM && module1!=module2) ok=kFALSE;
1840 if(ok){
6175da48 1841
1842 //Check if one of the clusters comes from a conversion
1843 if (p1->IsTagged() && p2->IsTagged()) fhReConv2->Fill(pt,m);
1844 else if(p1->IsTagged() || p2->IsTagged()) fhReConv ->Fill(pt,m);
1845
af7b3903 1846 //Fill histograms for different bad channel distance, centrality, assymmetry cut and pid bit
5ae09196 1847 for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
af7b3903 1848 if((p1->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton)) && (p2->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton))){
1849 for(Int_t iasym=0; iasym < fNAsymCuts; iasym++){
1850 if(a < fAsymCuts[iasym]){
1851 Int_t index = ((curCentrBin*fNPIDBits)+ipid)*fNAsymCuts + iasym;
6175da48 1852 //printf("index %d :(cen %d * nPID %d + ipid %d)*nasym %d + iasym %d\n",index,curCentrBin,fNPIDBits,ipid,fNAsymCuts,iasym);
af7b3903 1853 fhRe1 [index]->Fill(pt,m);
398c93cc 1854 if(fMakeInvPtPlots)fhReInvPt1[index]->Fill(pt,m,1./pt) ;
6175da48 1855 if(fFillBadDistHisto){
1856 if(p1->DistToBad()>0 && p2->DistToBad()>0){
1857 fhRe2 [index]->Fill(pt,m) ;
1858 if(fMakeInvPtPlots)fhReInvPt2[index]->Fill(pt,m,1./pt) ;
1859 if(p1->DistToBad()>1 && p2->DistToBad()>1){
1860 fhRe3 [index]->Fill(pt,m) ;
1861 if(fMakeInvPtPlots)fhReInvPt3[index]->Fill(pt,m,1./pt) ;
1862 }// bad 3
1863 }// bad2
1864 }// Fill bad dist histos
1865 }//assymetry cut
1866 }// asymmetry cut loop
af7b3903 1867 }// bad 1
1868 }// pid bit loop
5ae09196 1869
af7b3903 1870 //Fill histograms with opening angle
1871 fhRealOpeningAngle ->Fill(pt,angle);
1872 fhRealCosOpeningAngle->Fill(pt,TMath::Cos(angle));
1873
1874 //Fill histograms with pair assymmetry
1875 fhRePtAsym->Fill(pt,a);
6175da48 1876 if(m > 0.10 && m < 0.17) fhRePtAsymPi0->Fill(pt,a);
af7b3903 1877 if(m > 0.45 && m < 0.65) fhRePtAsymEta->Fill(pt,a);
1878
6175da48 1879 //-------------------------------------------------------
1880 //Get the number of cells needed for multi cut analysis.
1881 //-------------------------------------------------------
1882 Int_t ncell1 = 0;
1883 Int_t ncell2 = 0;
1884 if(fMultiCutAna || (IsDataMC() && fMultiCutAnaSim)){
1885
af7b3903 1886 AliVEvent * event = GetReader()->GetInputEvent();
1887 if(event){
1888 for(Int_t iclus = 0; iclus < event->GetNumberOfCaloClusters(); iclus++){
1889 AliVCluster *cluster = event->GetCaloCluster(iclus);
5ae09196 1890
af7b3903 1891 Bool_t is = kFALSE;
1892 if (fCalorimeter == "EMCAL" && GetReader()->IsEMCALCluster(cluster)) is = kTRUE;
1893 else if(fCalorimeter == "PHOS" && GetReader()->IsPHOSCluster (cluster)) is = kTRUE;
5ae09196 1894
af7b3903 1895 if(is){
1896 if (p1->GetCaloLabel(0) == cluster->GetID()) ncell1 = cluster->GetNCells();
1897 else if (p2->GetCaloLabel(0) == cluster->GetID()) ncell2 = cluster->GetNCells();
1898 } // PHOS or EMCAL cluster as requested in analysis
1899
1900 if(ncell2 > 0 && ncell1 > 0) break; // No need to continue the iteration
1901
1902 }
1903 //printf("e 1: %2.2f, e 2: %2.2f, ncells: n1 %d, n2 %d\n", p1->E(), p2->E(),ncell1,ncell2);
1904 }
6175da48 1905 }
1906
1907 //---------
1908 // MC data
1909 //---------
1910 //Do some MC checks on the origin of the pair, is there any common ancestor and if there is one, who?
1911 if(IsDataMC()) FillMCVersusRecDataHistograms(p1->GetLabel(), p2->GetLabel(),p1->Pt(), p2->Pt(),ncell1, ncell2, m, pt, a,deta, dphi);
1912
1913 //-----------------------
1914 //Multi cuts analysis
1915 //-----------------------
1916 if(fMultiCutAna){
1917 //Histograms for different PID bits selection
1918 for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
1919
1920 if(p1->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton) &&
1921 p2->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton)) fhRePIDBits[ipid]->Fill(pt,m) ;
1922
1923 //printf("ipt %d, ipid%d, name %s\n",ipt, ipid, fhRePtPIDCuts[ipt*fNPIDBitsBits+ipid]->GetName());
1924 } // pid bit cut loop
1925
1926 //Several pt,ncell and asymmetry cuts
af7b3903 1927 for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
1928 for(Int_t icell=0; icell<fNCellNCuts; icell++){
1929 for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
1930 Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
1931 if(p1->Pt() > fPtCuts[ipt] && p2->Pt() > fPtCuts[ipt] &&
1932 a < fAsymCuts[iasym] &&
6175da48 1933 ncell1 >= fCellNCuts[icell] && ncell2 >= fCellNCuts[icell]){
1934 fhRePtNCellAsymCuts[index]->Fill(pt,m) ;
1935 //printf("ipt %d, icell%d, iasym %d, name %s\n",ipt, icell, iasym, fhRePtNCellAsymCuts[((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym]->GetName());
1936 if(module1==module2){
1937 if (module1==0) fhRePtNCellAsymCutsSM0[index]->Fill(pt,m) ;
1938 else if(module1==1) fhRePtNCellAsymCutsSM1[index]->Fill(pt,m) ;
1939 else if(module1==2) fhRePtNCellAsymCutsSM2[index]->Fill(pt,m) ;
1940 else if(module1==3) fhRePtNCellAsymCutsSM3[index]->Fill(pt,m) ;
1941 else printf("AliAnaPi0::FillHistograms() - WRONG SM NUMBER\n");
1942 }
1943 }
af7b3903 1944 }// pid bit cut loop
1945 }// icell loop
1946 }// pt cut loop
1947 for(Int_t iasym = 0; iasym < fNAsymCuts; iasym++){
1948 if(a < fAsymCuts[iasym])fhRePtMult[iasym]->Fill(pt,GetTrackMultiplicity(),m) ;
1949 }
af7b3903 1950 }// multiple cuts analysis
1951 }// ok if same sm
7e7694bb 1952 }// second same event particle
1953 }// first cluster
6175da48 1954
1955 //-------------------------------------------------------------
1956 // Mixing
1957 //-------------------------------------------------------------
7e7694bb 1958 if(fDoOwnMix){
156549ae 1959 //printf("Cen bin %d, RP bin %d, e aver %f, mult %d\n",curCentrBin,curRPBin, eClusTot, nPhot);
6175da48 1960 //Recover events in with same characteristics as the current event
5025c139 1961 TList * evMixList=fEventsList[curCentrBin*GetNZvertBin()*GetNRPBin()+curZvertBin*GetNRPBin()+curRPBin] ;
7e7694bb 1962 Int_t nMixed = evMixList->GetSize() ;
1963 for(Int_t ii=0; ii<nMixed; ii++){
1964 TClonesArray* ev2= (TClonesArray*) (evMixList->At(ii));
1965 Int_t nPhot2=ev2->GetEntriesFast() ;
1966 Double_t m = -999;
6175da48 1967 if(GetDebug() > 1)
1968 printf("AliAnaPi0::MakeAnalysisFillHistograms() - Mixed event %d photon entries %d, centrality bin %d\n", ii, nPhot2, curCentrBin);
7e7694bb 1969
6175da48 1970 //---------------------------------
1971 //First loop on photons/clusters
1972 //---------------------------------
7e7694bb 1973 for(Int_t i1=0; i1<nPhot; i1++){
1974 AliAODPWG4Particle * p1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i1)) ;
6175da48 1975 if(fSameSM && GetModuleNumber(p1)!=module1) continue;
1976
1977 //Get kinematics of cluster and (super) module of this cluster
7e7694bb 1978 TLorentzVector photon1(p1->Px(),p1->Py(),p1->Pz(),p1->E());
af7b3903 1979 module1 = GetModuleNumber(p1);
6175da48 1980
1981 //---------------------------------
1982 //First loop on photons/clusters
1983 //---------------------------------
7e7694bb 1984 for(Int_t i2=0; i2<nPhot2; i2++){
1985 AliAODPWG4Particle * p2 = (AliAODPWG4Particle*) (ev2->At(i2)) ;
1986
6175da48 1987 //Get kinematics of second cluster and calculate those of the pair
7e7694bb 1988 TLorentzVector photon2(p2->Px(),p2->Py(),p2->Pz(),p2->E());
6175da48 1989 m = (photon1+photon2).M() ;
7e7694bb 1990 Double_t pt = (photon1 + photon2).Pt();
1991 Double_t a = TMath::Abs(p1->E()-p2->E())/(p1->E()+p2->E()) ;
1992
1993 //Check if opening angle is too large or too small compared to what is expected
1994 Double_t angle = photon1.Angle(photon2.Vect());
6175da48 1995 if(fUseAngleEDepCut && !GetNeutralMesonSelection()->IsAngleInWindow((photon1+photon2).E(),angle+0.05)){
1996 if(GetDebug() > 2)
1997 printf("AliAnaPi0::MakeAnalysisFillHistograms() -Mix pair angle %f not in E %f window\n",angle, (photon1+photon2).E());
1998 continue;
1999 }
2000 if(fUseAngleCut && (angle < fAngleCut || angle > fAngleMaxCut)) {
2001 if(GetDebug() > 2)
2002 printf("AliAnaPi0::MakeAnalysisFillHistograms() -Mix pair angle %f < cut %f\n",angle,fAngleCut);
2003 continue;
2004
2005 }
7e7694bb 2006
2007 if(GetDebug() > 2)
2008 printf("AliAnaPi0::MakeAnalysisFillHistograms() - Mixed Event: pT: photon1 %2.2f, photon2 %2.2f; Pair: pT %2.2f, mass %2.3f, a %f2.3\n",
af7b3903 2009 p1->Pt(), p2->Pt(), pt,m,a);
6175da48 2010
af7b3903 2011 //In case we want only pairs in same (super) module, check their origin.
2012 module2 = GetModuleNumber(p2);
6175da48 2013
2014 //-------------------------------------------------------------------------------------------------
2015 //Fill module dependent histograms, put a cut on assymmetry on the first available cut in the array
2016 //-------------------------------------------------------------------------------------------------
2017 if(a < fAsymCuts[0]){
2018 if(module1==module2 && module1 >=0 && module1<fNModules)
2019 fhMiMod[module1]->Fill(pt,m) ;
6175da48 2020
8d230fa8 2021 if(fCalorimeter=="EMCAL"){
2022
2023 // Same sector
2024 Int_t j=0;
2025 for(Int_t i = 0; i < fNModules/2; i++){
2026 j=2*i;
2027 if((module1==j && module2==j+1) || (module1==j+1 && module2==j)) fhMiSameSectorEMCALMod[i]->Fill(pt,m) ;
2028 }
2029
2030 // Same side
2031 for(Int_t i = 0; i < fNModules-2; i++){
2032 if((module1==i && module2==i+2) || (module1==i+2 && module2==i)) fhMiSameSideEMCALMod[i]->Fill(pt,m);
2033 }
2034 }//EMCAL
2035 else {//PHOS
2036 if((module1==0 && module2==1) || (module1==1 && module2==0)) fhMiDiffPHOSMod[0]->Fill(pt,m) ;
2037 if((module1==0 && module2==2) || (module1==2 && module2==0)) fhMiDiffPHOSMod[1]->Fill(pt,m) ;
2038 if((module1==1 && module2==2) || (module1==2 && module2==1)) fhMiDiffPHOSMod[2]->Fill(pt,m) ;
2039 }//PHOS
2040
2041
6175da48 2042 }
2043
af7b3903 2044 Bool_t ok = kTRUE;
2045 if(fSameSM && module1!=module2) ok=kFALSE;
2046 if(ok){
6175da48 2047
2048 //Check if one of the clusters comes from a conversion
2049 if (p1->IsTagged() && p2->IsTagged()) fhMiConv2->Fill(pt,m);
2050 else if(p1->IsTagged() || p2->IsTagged()) fhMiConv ->Fill(pt,m);
2051
2052 //Fill histograms for different bad channel distance, centrality, assymmetry cut and pid bit
af7b3903 2053 for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
2054 if((p1->IsPIDOK(ipid,AliCaloPID::kPhoton)) && (p2->IsPIDOK(ipid,AliCaloPID::kPhoton))){
2055 for(Int_t iasym=0; iasym < fNAsymCuts; iasym++){
2056 if(a < fAsymCuts[iasym]){
2057 Int_t index = ((curCentrBin*fNPIDBits)+ipid)*fNAsymCuts + iasym;
2058 fhMi1 [index]->Fill(pt,m) ;
398c93cc 2059 if(fMakeInvPtPlots)fhMiInvPt1[index]->Fill(pt,m,1./pt) ;
6175da48 2060 if(fFillBadDistHisto){
2061 if(p1->DistToBad()>0 && p2->DistToBad()>0){
2062 fhMi2 [index]->Fill(pt,m) ;
2063 if(fMakeInvPtPlots)fhMiInvPt2[index]->Fill(pt,m,1./pt) ;
2064 if(p1->DistToBad()>1 && p2->DistToBad()>1){
2065 fhMi3 [index]->Fill(pt,m) ;
2066 if(fMakeInvPtPlots)fhMiInvPt3[index]->Fill(pt,m,1./pt) ;
2067 }
af7b3903 2068 }
6175da48 2069 }// Fill bad dist histo
af7b3903 2070 }//Asymmetry cut
2071 }// Asymmetry loop
2072 }//PID cut
2073 }//loop for histograms
6175da48 2074
2075 //-----------------------
2076 //Multi cuts analysis
2077 //-----------------------
2078 if(fMultiCutAna){
2079 //Several pt,ncell and asymmetry cuts
2080 for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
2081 for(Int_t icell=0; icell<fNCellNCuts; icell++){
2082 for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
2083 Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
2084 if(p1->Pt() > fPtCuts[ipt] && p2->Pt() > fPtCuts[ipt] &&
2085 a < fAsymCuts[iasym] &&
2086 p1->GetBtag() >= fCellNCuts[icell] && p2->GetBtag() >= fCellNCuts[icell]){
2087 fhMiPtNCellAsymCuts[index]->Fill(pt,m) ;
2088 //printf("ipt %d, icell%d, iasym %d, name %s\n",ipt, icell, iasym, fhRePtNCellAsymCuts[((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym]->GetName());
2089 }
2090 }// pid bit cut loop
2091 }// icell loop
2092 }// pt cut loop
2093 } // Multi cut ana
2094
2095 //Fill histograms with opening angle
2096 fhMixedOpeningAngle ->Fill(pt,angle);
2097 fhMixedCosOpeningAngle->Fill(pt,TMath::Cos(angle));
af7b3903 2098 }//ok
7e7694bb 2099 }// second cluster loop
2100 }//first cluster loop
2101 }//loop on mixed events
2102
6175da48 2103 //--------------------------------------------------------
2104 //Add the current event to the list of events for mixing
2105 //--------------------------------------------------------
7e7694bb 2106 TClonesArray *currentEvent = new TClonesArray(*GetInputAODBranch());
af7b3903 2107 //Add current event to buffer and Remove redundant events
7e7694bb 2108 if(currentEvent->GetEntriesFast()>0){
2109 evMixList->AddFirst(currentEvent) ;
2110 currentEvent=0 ; //Now list of particles belongs to buffer and it will be deleted with buffer
2111 if(evMixList->GetSize()>=fNmaxMixEv)
2112 {
2113 TClonesArray * tmp = (TClonesArray*) (evMixList->Last()) ;
2114 evMixList->RemoveLast() ;
2115 delete tmp ;
2116 }
2117 }
2118 else{ //empty event
2119 delete currentEvent ;
2120 currentEvent=0 ;
477d6cee 2121 }
7e7694bb 2122 }// DoOwnMix
c8fe2783 2123
1c5acb87 2124}
2125
a5cc4f03 2126//________________________________________________________________________
2127void AliAnaPi0::ReadHistograms(TList* outputList)
2128{
50f39b97 2129 // Needed when Terminate is executed in distributed environment
2130 // Refill analysis histograms of this class with corresponding histograms in output list.
2131
2132 // Histograms of this analsys are kept in the same list as other analysis, recover the position of
2133 // the first one and then add the next.
2134 Int_t index = outputList->IndexOf(outputList->FindObject(GetAddedHistogramsStringToName()+"hRe_cen0_pid0_dist1"));
2135
af7b3903 2136 if(!fhRe1) fhRe1 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
2137 if(!fhRe2) fhRe2 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
2138 if(!fhRe3) fhRe3 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
2139 if(!fhMi1) fhMi1 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
2140 if(!fhMi2) fhMi2 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
2141 if(!fhMi3) fhMi3 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
6175da48 2142 if(!fhReInvPt1) fhReInvPt1 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
2143 if(!fhReInvPt2) fhReInvPt2 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
2144 if(!fhReInvPt3) fhReInvPt3 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
2145 if(!fhMiInvPt1) fhMiInvPt1 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
2146 if(!fhMiInvPt2) fhMiInvPt2 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
2147 if(!fhMiInvPt3) fhMiInvPt3 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
af7b3903 2148 if(!fhReMod) fhReMod = new TH2D*[fNModules] ;
8d230fa8 2149 if(!fhReDiffPHOSMod) fhReDiffPHOSMod = new TH2D*[fNModules] ;
2150 if(!fhReSameSectorEMCALMod)fhReSameSectorEMCALMod = new TH2D*[fNModules/2] ;
2151 if(!fhReSameSideEMCALMod) fhReSameSideEMCALMod = new TH2D*[fNModules-2] ;
ad30b142 2152 if(!fhMiMod) fhMiMod = new TH2D*[fNModules] ;
8d230fa8 2153 if(!fhMiDiffPHOSMod) fhMiDiffPHOSMod = new TH2D*[fNModules] ;
2154 if(!fhMiSameSectorEMCALMod)fhMiSameSectorEMCALMod = new TH2D*[fNModules/2] ;
2155 if(!fhMiSameSideEMCALMod) fhMiSameSideEMCALMod = new TH2D*[fNModules-2] ;
2156
6175da48 2157 fhReConv = (TH2D*) outputList->At(index++);
2158 fhMiConv = (TH2D*) outputList->At(index++);
2159 fhReConv2 = (TH2D*) outputList->At(index++);
2160 fhMiConv2 = (TH2D*) outputList->At(index++);
821c8090 2161
50f39b97 2162 for(Int_t ic=0; ic<fNCentrBin; ic++){
af7b3903 2163 for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
2164 for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
2165 Int_t ihisto = ((ic*fNPIDBits)+ipid)*fNAsymCuts + iasym;
2166
2167 fhRe1[ihisto] = (TH2D*) outputList->At(index++);
2168 fhRe2[ihisto] = (TH2D*) outputList->At(index++);
2169 fhRe3[ihisto] = (TH2D*) outputList->At(index++);
6175da48 2170
2171 fhReInvPt1[ihisto] = (TH2D*) outputList->At(index++);
2172 fhReInvPt2[ihisto] = (TH2D*) outputList->At(index++);
2173 fhReInvPt3[ihisto] = (TH2D*) outputList->At(index++);
5ae09196 2174
af7b3903 2175 if(fDoOwnMix){
2176 fhMi1[ihisto] = (TH2D*) outputList->At(index++);
2177 fhMi2[ihisto] = (TH2D*) outputList->At(index++);
2178 fhMi3[ihisto] = (TH2D*) outputList->At(index++);
6175da48 2179
2180 fhMiInvPt1[ihisto] = (TH2D*) outputList->At(index++);
2181 fhMiInvPt2[ihisto] = (TH2D*) outputList->At(index++);
2182 fhMiInvPt3[ihisto] = (TH2D*) outputList->At(index++);
af7b3903 2183 }//Own mix
2184 }//asymmetry loop
2185 }// pid loop
2186 }// centrality loop
2187
2188 fhRePtAsym = (TH2D*)outputList->At(index++);
2189 fhRePtAsymPi0 = (TH2D*)outputList->At(index++);
2190 fhRePtAsymEta = (TH2D*)outputList->At(index++);
eee5fcf1 2191
5ae09196 2192 if(fMultiCutAna){
2193
eee5fcf1 2194 if(!fhRePtNCellAsymCuts) fhRePtNCellAsymCuts = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
2195 if(!fhRePIDBits) fhRePIDBits = new TH2D*[fNPIDBits];
2196
5ae09196 2197 for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
2198 fhRePIDBits[ipid] = (TH2D*) outputList->At(index++);
2199 }// ipid loop
2200
2201 for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
2202 for(Int_t icell=0; icell<fNCellNCuts; icell++){
2203 for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
2204 fhRePtNCellAsymCuts[((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym] = (TH2D*) outputList->At(index++);
2205 }// iasym
2206 }// icell loop
2207 }// ipt loop
af7b3903 2208
2209 if(!fhRePtMult) fhRePtMult = new TH3D*[fNAsymCuts] ;
2210 for(Int_t iasym = 0; iasym < fNAsymCuts; iasym++)
2211 fhRePtMult[iasym] = (TH3D*) outputList->At(index++);
5ae09196 2212 }// multi cut analysis
50f39b97 2213
745f04da 2214 fhEvents = (TH3D *) outputList->At(index++);
2215 fhCentrality = (TH1D *) outputList->At(index++);
2216
af7b3903 2217 fhRealOpeningAngle = (TH2D*) outputList->At(index++);
2218 fhRealCosOpeningAngle = (TH2D*) outputList->At(index++);
6175da48 2219 if(fDoOwnMix){
2220 fhMixedOpeningAngle = (TH2D*) outputList->At(index++);
2221 fhMixedCosOpeningAngle = (TH2D*) outputList->At(index++);
2222 }
af7b3903 2223
50f39b97 2224 //Histograms filled only if MC data is requested
2225 if(IsDataMC() || (GetReader()->GetDataType() == AliCaloTrackReader::kMC) ){
6175da48 2226 fhPrimPi0Pt = (TH1D*) outputList->At(index++);
2227 fhPrimPi0AccPt = (TH1D*) outputList->At(index++);
156549ae 2228 fhPrimPi0Y = (TH2D*) outputList->At(index++);
2229 fhPrimPi0AccY = (TH2D*) outputList->At(index++);
2230 fhPrimPi0Phi = (TH2D*) outputList->At(index++);
2231 fhPrimPi0AccPhi = (TH2D*) outputList->At(index++);
2232 fhPrimEtaPt = (TH1D*) outputList->At(index++);
2233 fhPrimEtaAccPt = (TH1D*) outputList->At(index++);
2234 fhPrimEtaY = (TH2D*) outputList->At(index++);
2235 fhPrimEtaAccY = (TH2D*) outputList->At(index++);
2236 fhPrimEtaPhi = (TH2D*) outputList->At(index++);
2237 fhPrimEtaAccPhi = (TH2D*) outputList->At(index++);
6175da48 2238 for(Int_t i = 0; i<13; i++){
2239 fhMCOrgMass[i] = (TH2D*) outputList->At(index++);
2240 fhMCOrgAsym[i] = (TH2D*) outputList->At(index++);
2241 fhMCOrgDeltaEta[i] = (TH2D*) outputList->At(index++);
2242 fhMCOrgDeltaPhi[i] = (TH2D*) outputList->At(index++);
2243 }
2244
2245 if(fMultiCutAnaSim){
2246 fhMCPi0MassPtTrue = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
2247 fhMCPi0MassPtRec = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
2248 fhMCPi0PtTruePtRec = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
2249 fhMCEtaMassPtTrue = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
2250 fhMCEtaMassPtRec = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
2251 fhMCEtaPtTruePtRec = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
2252 for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
2253 for(Int_t icell=0; icell<fNCellNCuts; icell++){
2254 for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
2255 Int_t in = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
2256 fhMCPi0MassPtTrue[in] = (TH2D*) outputList->At(index++);
2257 fhMCPi0PtTruePtRec[in] = (TH2D*) outputList->At(index++);
2258 fhMCEtaMassPtTrue[in] = (TH2D*) outputList->At(index++);
2259 fhMCEtaPtTruePtRec[in] = (TH2D*) outputList->At(index++);
2260 }
2261 }
2262 }
2263 }
2264 else{
2265 fhMCPi0MassPtTrue = new TH2D*[1];
2266 fhMCPi0PtTruePtRec = new TH2D*[1];
2267 fhMCEtaMassPtTrue = new TH2D*[1];
2268 fhMCEtaPtTruePtRec = new TH2D*[1];
2269
2270 fhMCPi0MassPtTrue[0] = (TH2D*) outputList->At(index++);
2271 fhMCPi0PtTruePtRec[0] = (TH2D*) outputList->At(index++);
2272 fhMCEtaMassPtTrue[0] = (TH2D*) outputList->At(index++);
2273 fhMCEtaPtTruePtRec[0] = (TH2D*) outputList->At(index++);
2274 }
50f39b97 2275 }
2276
6175da48 2277 for(Int_t imod=0; imod < fNModules; imod++){
8d230fa8 2278 fhReMod[imod] = (TH2D*) outputList->At(index++);
2279 if(fCalorimeter=="EMCAL"){
2280 if(imod < fNModules/2) fhReSameSectorEMCALMod[imod] = (TH2D*) outputList->At(index++);
2281 if(imod < fNModules-2) fhReSameSideEMCALMod[imod] = (TH2D*) outputList->At(index++);
2282 }
2283 else fhReDiffPHOSMod[imod] = (TH2D*) outputList->At(index++);
2284
6175da48 2285 if(fDoOwnMix){
8d230fa8 2286 fhMiMod[imod] = (TH2D*) outputList->At(index++);
2287 if(fCalorimeter=="EMCAL"){
2288 if(imod < fNModules/2) fhMiSameSectorEMCALMod[imod] = (TH2D*) outputList->At(index++);
2289 if(imod < fNModules-2) fhMiSameSideEMCALMod[imod] = (TH2D*) outputList->At(index++);
2290 }
2291 else fhMiDiffPHOSMod[imod] = (TH2D*) outputList->At(index++);
6175da48 2292 }
2293 }
eee5fcf1 2294
a5cc4f03 2295}
2296
2297
6639984f 2298//____________________________________________________________________________________________________________________________________________________
a5cc4f03 2299void AliAnaPi0::Terminate(TList* outputList)
6639984f 2300{
2301 //Do some calculations and plots from the final histograms.
477d6cee 2302
fbeaf916 2303 printf(" *** %s Terminate:\n", GetName()) ;
50f39b97 2304
a5cc4f03 2305 //Recover histograms from output histograms list, needed for distributed analysis.
2306 ReadHistograms(outputList);
50f39b97 2307
2e557d1c 2308 if (!fhRe1) {
50f39b97 2309 printf("AliAnaPi0::Terminate() - Error: Remote output histograms not imported in AliAnaPi0 object");
2310 return;
2e557d1c 2311 }
50f39b97 2312
a3aebfff 2313 printf("AliAnaPi0::Terminate() Mgg Real : %5.3f , RMS : %5.3f \n", fhRe1[0]->GetMean(), fhRe1[0]->GetRMS() ) ;
5ae09196 2314
2315 const Int_t buffersize = 255;
2316
2317 char nameIM[buffersize];
2318 snprintf(nameIM, buffersize,"AliAnaPi0_%s_cPt",fCalorimeter.Data());
71dd883b 2319 TCanvas * cIM = new TCanvas(nameIM, "", 400, 10, 600, 700) ;
6639984f 2320 cIM->Divide(2, 2);
50f39b97 2321
6639984f 2322 cIM->cd(1) ;
2323 //gPad->SetLogy();
af7b3903 2324 TH1D * hIMAllPt = (TH1D*) fhRe1[0]->ProjectionY(Form("IMPtAll_%s",fCalorimeter.Data()));
6639984f 2325 hIMAllPt->SetLineColor(2);
2326 hIMAllPt->SetTitle("No cut on p_{T, #gamma#gamma} ");
2327 hIMAllPt->Draw();
2328
2329 cIM->cd(2) ;
af7b3903 2330 TH1D * hIMPt5 = (TH1D*) fhRe1[0]->ProjectionY(Form("IMPt0-5_%s",fCalorimeter.Data()),0, fhRe1[0]->GetXaxis()->FindBin(5.));
2244659d 2331// hRe1Pt5->GetXaxis()->SetRangeUser(0,5);
2332// TH1D * hIMPt5 = (TH1D*) hRe1Pt5->Project3D(Form("IMPt5_%s_pz",fCalorimeter.Data()));
6639984f 2333 hIMPt5->SetLineColor(2);
2334 hIMPt5->SetTitle("0 < p_{T, #gamma#gamma} < 5 GeV/c");
2335 hIMPt5->Draw();
2336
2337 cIM->cd(3) ;
af7b3903 2338 TH1D * hIMPt10 = (TH1D*) fhRe1[0]->ProjectionY(Form("IMPt5-10_%s",fCalorimeter.Data()), fhRe1[0]->GetXaxis()->FindBin(5.),fhRe1[0]->GetXaxis()->FindBin(10.));
2244659d 2339// hRe1Pt10->GetXaxis()->SetRangeUser(5,10);
2340// TH1D * hIMPt10 = (TH1D*) hRe1Pt10->Project3D(Form("IMPt10_%s_pz",fCalorimeter.Data()));
6639984f 2341 hIMPt10->SetLineColor(2);
2342 hIMPt10->SetTitle("5 < p_{T, #gamma#gamma} < 10 GeV/c");
2343 hIMPt10->Draw();
2344
2345 cIM->cd(4) ;
af7b3903 2346 TH1D * hIMPt20 = (TH1D*) fhRe1[0]->ProjectionY(Form("IMPt10-20_%s",fCalorimeter.Data()), fhRe1[0]->GetXaxis()->FindBin(10.),fhRe1[0]->GetXaxis()->FindBin(20.));
2244659d 2347 // TH3F * hRe1Pt20 = (TH3F*)fhRe1[0]->Clone(Form("IMPt20_%s",fCalorimeter.Data()));
2348// hRe1Pt20->GetXaxis()->SetRangeUser(10,20);
2349// TH1D * hIMPt20 = (TH1D*) hRe1Pt20->Project3D(Form("IMPt20_%s_pz",fCalorimeter.Data()));
6639984f 2350 hIMPt20->SetLineColor(2);
2351 hIMPt20->SetTitle("10 < p_{T, #gamma#gamma} < 20 GeV/c");
2352 hIMPt20->Draw();
2353
5ae09196 2354 char nameIMF[buffersize];
2355 snprintf(nameIMF,buffersize,"AliAnaPi0_%s_Mgg.eps",fCalorimeter.Data());
71dd883b 2356 cIM->Print(nameIMF);
6639984f 2357
5ae09196 2358 char namePt[buffersize];
2359 snprintf(namePt,buffersize,"AliAnaPi0_%s_cPt",fCalorimeter.Data());
71dd883b 2360 TCanvas * cPt = new TCanvas(namePt, "", 400, 10, 600, 700) ;
6639984f 2361 cPt->Divide(2, 2);
2362
2363 cPt->cd(1) ;
2364 //gPad->SetLogy();
af7b3903 2365 TH1D * hPt = (TH1D*) fhRe1[0]->ProjectionX(Form("Pt0_%s",fCalorimeter.Data()),-1,-1);
6639984f 2366 hPt->SetLineColor(2);
2367 hPt->SetTitle("No cut on M_{#gamma#gamma} ");
2368 hPt->Draw();
2369
2370 cPt->cd(2) ;
af7b3903 2371 TH1D * hPtIM1 = (TH1D*)fhRe1[0]->ProjectionX(Form("Pt1_%s",fCalorimeter.Data()), fhRe1[0]->GetZaxis()->FindBin(0.05),fhRe1[0]->GetZaxis()->FindBin(0.21));
2244659d 2372// TH3F * hRe1IM1 = (TH3F*)fhRe1[0]->Clone(Form("Pt1_%s",fCalorimeter.Data()));
2373// hRe1IM1->GetZaxis()->SetRangeUser(0.05,0.21);
2374// TH1D * hPtIM1 = (TH1D*) hRe1IM1->Project3D("x");
6639984f 2375 hPtIM1->SetLineColor(2);
2376 hPtIM1->SetTitle("0.05 < M_{#gamma#gamma} < 0.21 GeV/c^{2}");
2377 hPtIM1->Draw();
2378
2379 cPt->cd(3) ;
af7b3903 2380 TH1D * hPtIM2 = (TH1D*)fhRe1[0]->ProjectionX(Form("Pt2_%s",fCalorimeter.Data()), fhRe1[0]->GetZaxis()->FindBin(0.09),fhRe1[0]->GetZaxis()->FindBin(0.17));
2244659d 2381// TH3F * hRe1IM2 = (TH3F*)fhRe1[0]->Clone(Form("Pt2_%s",fCalorimeter.Data()));
2382// hRe1IM2->GetZaxis()->SetRangeUser(0.09,0.17);
2383// TH1D * hPtIM2 = (TH1D*) hRe1IM2->Project3D("x");
6639984f 2384 hPtIM2->SetLineColor(2);
2385 hPtIM2->SetTitle("0.09 < M_{#gamma#gamma} < 0.17 GeV/c^{2}");
2386 hPtIM2->Draw();
2387
2388 cPt->cd(4) ;
af7b3903 2389 TH1D * hPtIM3 = (TH1D*)fhRe1[0]->ProjectionX(Form("Pt3_%s",fCalorimeter.Data()), fhRe1[0]->GetZaxis()->FindBin(0.11),fhRe1[0]->GetZaxis()->FindBin(0.15));
2244659d 2390// TH3F * hRe1IM3 = (TH3F*)fhRe1[0]->Clone(Form("Pt3_%s",fCalorimeter.Data()));
2391// hRe1IM3->GetZaxis()->SetRangeUser(0.11,0.15);
2392// TH1D * hPtIM3 = (TH1D*) hRe1IM1->Project3D("x");
6639984f 2393 hPtIM3->SetLineColor(2);
2394 hPtIM3->SetTitle("0.11 < M_{#gamma#gamma} < 0.15 GeV/c^{2}");
2395 hPtIM3->Draw();
2396
164a1d84 2397 char namePtF[buffersize];
5ae09196 2398 snprintf(namePtF,buffersize,"AliAnaPi0_%s_Pt.eps",fCalorimeter.Data());
71dd883b 2399 cPt->Print(namePtF);
1c5acb87 2400
5ae09196 2401 char line[buffersize] ;
2402 snprintf(line,buffersize,".!tar -zcf %s_%s.tar.gz *.eps", GetName(),fCalorimeter.Data()) ;
6639984f 2403 gROOT->ProcessLine(line);
5ae09196 2404 snprintf(line, buffersize,".!rm -fR AliAnaPi0_%s*.eps",fCalorimeter.Data());
6639984f 2405 gROOT->ProcessLine(line);
2406
71dd883b 2407 printf(" AliAnaPi0::Terminate() - !! All the eps files are in %s_%s.tar.gz !!!\n", GetName(), fCalorimeter.Data());
1c5acb87 2408
6639984f 2409}
c8fe2783 2410 //____________________________________________________________________________________________________________________________________________________
2411Int_t AliAnaPi0::GetEventIndex(AliAODPWG4Particle * part, Double_t * vert)
2412{
f8006433 2413 // retieves the event index and checks the vertex
2414 // in the mixed buffer returns -2 if vertex NOK
2415 // for normal events returns 0 if vertex OK and -1 if vertex NOK
2416
2417 Int_t evtIndex = -1 ;
2418 if(GetReader()->GetDataType()!=AliCaloTrackReader::kMC){
2419
2420 if (GetMixedEvent()){
2421
2422 evtIndex = GetMixedEvent()->EventIndexForCaloCluster(part->GetCaloLabel(0)) ;
2423 GetVertex(vert,evtIndex);
2424
5025c139 2425 if(TMath::Abs(vert[2])> GetZvertexCut())
f8006433 2426 evtIndex = -2 ; //Event can not be used (vertex, centrality,... cuts not fulfilled)
2427 } else {// Single event
2428
2429 GetVertex(vert);
2430
5025c139 2431 if(TMath::Abs(vert[2])> GetZvertexCut())
f8006433 2432 evtIndex = -1 ; //Event can not be used (vertex, centrality,... cuts not fulfilled)
2433 else
2434 evtIndex = 0 ;
c8fe2783 2435 }
0ae57829 2436 }//No MC reader
f8006433 2437 else {
2438 evtIndex = 0;
2439 vert[0] = 0. ;
2440 vert[1] = 0. ;
2441 vert[2] = 0. ;
2442 }
0ae57829 2443
f8006433 2444 return evtIndex ;
c8fe2783 2445}
f8006433 2446