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