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