]> 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
Be sure to load mapping when needed
[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;
c4a7d28a 1349 TVector3 prodVertex;
6175da48 1350 Int_t ancLabel = GetMCAnalysisUtils()->CheckCommonAncestor(index1, index2,
c4a7d28a 1351 GetReader(), ancPDG, ancStatus,ancMomentum, prodVertex);
6175da48 1352
08a56f5f 1353 Int_t momindex = -1;
1354 Int_t mompdg = -1;
1355 Int_t momstatus = -1;
6175da48 1356 if(GetDebug() > 1) printf("AliAnaPi0::FillMCVersusRecDataHistograms() - Common ancestor label %d, pdg %d, name %s, status %d; \n",
1357 ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus);
1358
1359 if(ancLabel > -1){
1360 if(ancPDG==22){//gamma
1361 fhMCOrgMass[0]->Fill(pt,mass);
1362 fhMCOrgAsym[0]->Fill(pt,asym);
1363 fhMCOrgDeltaEta[0]->Fill(pt,deta);
1364 fhMCOrgDeltaPhi[0]->Fill(pt,dphi);
1365 }
1366 else if(TMath::Abs(ancPDG)==11){//e
1367 fhMCOrgMass[1]->Fill(pt,mass);
1368 fhMCOrgAsym[1]->Fill(pt,asym);
1369 fhMCOrgDeltaEta[1]->Fill(pt,deta);
1370 fhMCOrgDeltaPhi[1]->Fill(pt,dphi);
1371 }
1372 else if(ancPDG==111){//Pi0
91e1ea12 1373 fhMCOrgMass[2]->Fill(pt,mass);
1374 fhMCOrgAsym[2]->Fill(pt,asym);
1375 fhMCOrgDeltaEta[2]->Fill(pt,deta);
1376 fhMCOrgDeltaPhi[2]->Fill(pt,dphi);
6175da48 1377 if(fMultiCutAnaSim){
1378 for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
1379 for(Int_t icell=0; icell<fNCellNCuts; icell++){
1380 for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
1381 Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
1382 if(pt1 > fPtCuts[ipt] && pt2 > fPtCuts[ipt] &&
1383 asym < fAsymCuts[iasym] &&
1384 ncell1 >= fCellNCuts[icell] && ncell2 >= fCellNCuts[icell]){
91e1ea12 1385 fhMCPi0MassPtRec [index]->Fill(pt,mass);
1386 fhMCPi0MassPtTrue[index]->Fill(ancMomentum.Pt(),mass);
1387 if(mass < 0.17 && mass > 0.1) fhMCPi0PtTruePtRec[index]->Fill(ancMomentum.Pt(),pt);
6175da48 1388 }//pass the different cuts
1389 }// pid bit cut loop
1390 }// icell loop
1391 }// pt cut loop
1392 }//Multi cut ana sim
1393 else {
91e1ea12 1394 fhMCPi0MassPtTrue[0]->Fill(ancMomentum.Pt(),mass);
08a56f5f 1395 if(mass < 0.17 && mass > 0.1) {
91e1ea12 1396 fhMCPi0PtTruePtRec[0]->Fill(ancMomentum.Pt(),pt);
08a56f5f 1397
1398 if(GetReader()->ReadStack()){
1399 TParticle* ancestor = GetMCStack()->Particle(ancLabel);
1400 momindex = ancestor->GetFirstMother();
41121cfe 1401 if(momindex < 0) return;
08a56f5f 1402 TParticle* mother = GetMCStack()->Particle(momindex);
1403 mompdg = TMath::Abs(mother->GetPdgCode());
1404 momstatus = mother->GetStatusCode();
1405 }
1406 else {
1407 TClonesArray * mcparticles = GetReader()->GetAODMCParticles(0);
1408 AliAODMCParticle* ancestor = (AliAODMCParticle *) mcparticles->At(ancLabel);
1409 momindex = ancestor->GetMother();
41121cfe 1410 if(momindex < 0) return;
08a56f5f 1411 AliAODMCParticle* mother = (AliAODMCParticle *) mcparticles->At(momindex);
1412 mompdg = TMath::Abs(mother->GetPdgCode());
1413 momstatus = mother->GetStatus();
1414 }
1415
1416 if (momstatus == 21) fhMCPi0PtOrigin->Fill(pt,0.5);//parton
1417 else if(mompdg < 22 ) fhMCPi0PtOrigin->Fill(pt,1.5);//quark
1418 else if(mompdg > 2100 && mompdg < 2210) fhMCPi0PtOrigin->Fill(pt,2.5);// resonances
1419 else if(mompdg == 221) fhMCPi0PtOrigin->Fill(pt,8.5);//eta
1420 else if(mompdg == 331) fhMCPi0PtOrigin->Fill(pt,9.5);//eta prime
1421 else if(mompdg == 213) fhMCPi0PtOrigin->Fill(pt,4.5);//rho
1422 else if(mompdg == 223) fhMCPi0PtOrigin->Fill(pt,5.5);//omega
1423 else if(mompdg >= 310 && mompdg <= 323) fhMCPi0PtOrigin->Fill(pt,6.5);//k0S, k+-,k*
1424 else if(mompdg == 130) fhMCPi0PtOrigin->Fill(pt,6.5);//k0L
1425 else if(momstatus == 11 || momstatus == 12 ) fhMCPi0PtOrigin->Fill(pt,3.5);//resonances
91e1ea12 1426 else fhMCPi0PtOrigin->Fill(pt,7.5);//other?
08a56f5f 1427
1428 }//pi0 mass region
1429
6175da48 1430 }
1431 }
1432 else if(ancPDG==221){//Eta
1433 fhMCOrgMass[3]->Fill(pt,mass);
1434 fhMCOrgAsym[3]->Fill(pt,asym);
1435 fhMCOrgDeltaEta[3]->Fill(pt,deta);
1436 fhMCOrgDeltaPhi[3]->Fill(pt,dphi);
1437 if(fMultiCutAnaSim){
1438 for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
1439 for(Int_t icell=0; icell<fNCellNCuts; icell++){
1440 for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
1441 Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
1442 if(pt1 > fPtCuts[ipt] && pt2 > fPtCuts[ipt] &&
1443 asym < fAsymCuts[iasym] &&
1444 ncell1 >= fCellNCuts[icell] && ncell2 >= fCellNCuts[icell]){
1445 fhMCEtaMassPtRec [index]->Fill(pt,mass);
1446 fhMCEtaMassPtTrue[index]->Fill(ancMomentum.Pt(),mass);
156549ae 1447 if(mass < 0.65 && mass > 0.45) fhMCEtaPtTruePtRec[index]->Fill(ancMomentum.Pt(),pt);
6175da48 1448 }//pass the different cuts
1449 }// pid bit cut loop
1450 }// icell loop
1451 }// pt cut loop
1452 } //Multi cut ana sim
1453 else {
1454 fhMCEtaMassPtTrue[0]->Fill(ancMomentum.Pt(),mass);
156549ae 1455 if(mass < 0.65 && mass > 0.45) fhMCEtaPtTruePtRec[0]->Fill(ancMomentum.Pt(),pt);
08a56f5f 1456
1457 if(GetReader()->ReadStack()){
1458 TParticle* ancestor = GetMCStack()->Particle(ancLabel);
1459 momindex = ancestor->GetFirstMother();
41121cfe 1460 if(momindex < 0) return;
08a56f5f 1461 TParticle* mother = GetMCStack()->Particle(momindex);
1462 mompdg = TMath::Abs(mother->GetPdgCode());
1463 momstatus = mother->GetStatusCode();
1464 }
1465 else {
1466 TClonesArray * mcparticles = GetReader()->GetAODMCParticles(0);
1467 AliAODMCParticle* ancestor = (AliAODMCParticle *) mcparticles->At(ancLabel);
1468 momindex = ancestor->GetMother();
41121cfe 1469 if(momindex < 0) return;
08a56f5f 1470 AliAODMCParticle* mother = (AliAODMCParticle *) mcparticles->At(momindex);
1471 mompdg = TMath::Abs(mother->GetPdgCode());
1472 momstatus = mother->GetStatus();
1473 }
1474
1475 if (momstatus == 21 ) fhMCEtaPtOrigin->Fill(pt,0.5);//parton
1476 else if(mompdg < 22 ) fhMCEtaPtOrigin->Fill(pt,1.5);//quark
1477 else if(mompdg > 2100 && mompdg < 2210) fhMCEtaPtOrigin->Fill(pt,2.5);//qq resonances
1478 else if(mompdg == 331) fhMCEtaPtOrigin->Fill(pt,5.5);//eta prime
1479 else if(momstatus == 11 || momstatus == 12 ) fhMCEtaPtOrigin->Fill(pt,3.5);//resonances
1480 else fhMCEtaPtOrigin->Fill(pt,4.5);//stable, conversions?
1481 //printf("Other Meson pdg %d, Mother %s, pdg %d, status %d\n",pdg, TDatabasePDG::Instance()->GetParticle(mompdg)->GetName(),mompdg, momstatus );
1482 }// eta mass region
6175da48 1483 }
1484 else if(ancPDG==-2212){//AProton
1485 fhMCOrgMass[4]->Fill(pt,mass);
1486 fhMCOrgAsym[4]->Fill(pt,asym);
1487 fhMCOrgDeltaEta[4]->Fill(pt,deta);
1488 fhMCOrgDeltaPhi[4]->Fill(pt,dphi);
1489 }
1490 else if(ancPDG==-2112){//ANeutron
1491 fhMCOrgMass[5]->Fill(pt,mass);
1492 fhMCOrgAsym[5]->Fill(pt,asym);
1493 fhMCOrgDeltaEta[5]->Fill(pt,deta);
1494 fhMCOrgDeltaPhi[5]->Fill(pt,dphi);
1495 }
1496 else if(TMath::Abs(ancPDG)==13){//muons
1497 fhMCOrgMass[6]->Fill(pt,mass);
1498 fhMCOrgAsym[6]->Fill(pt,asym);
1499 fhMCOrgDeltaEta[6]->Fill(pt,deta);
1500 fhMCOrgDeltaPhi[6]->Fill(pt,dphi);
1501 }
1502 else if (TMath::Abs(ancPDG) > 100 && ancLabel > 7) {
1503 if(ancStatus==1){//Stable particles, converted? not decayed resonances
1504 fhMCOrgMass[6]->Fill(pt,mass);
1505 fhMCOrgAsym[6]->Fill(pt,asym);
1506 fhMCOrgDeltaEta[6]->Fill(pt,deta);
1507 fhMCOrgDeltaPhi[6]->Fill(pt,dphi);
1508 }
1509 else{//resonances and other decays, more hadron conversions?
1510 fhMCOrgMass[7]->Fill(pt,mass);
1511 fhMCOrgAsym[7]->Fill(pt,asym);
1512 fhMCOrgDeltaEta[7]->Fill(pt,deta);
1513 fhMCOrgDeltaPhi[7]->Fill(pt,dphi);
1514 }
1515 }
1516 else {//Partons, colliding protons, strings, intermediate corrections
1517 if(ancStatus==11 || ancStatus==12){//String fragmentation
1518 fhMCOrgMass[8]->Fill(pt,mass);
1519 fhMCOrgAsym[8]->Fill(pt,asym);
1520 fhMCOrgDeltaEta[8]->Fill(pt,deta);
1521 fhMCOrgDeltaPhi[8]->Fill(pt,dphi);
1522 }
1523 else if (ancStatus==21){
1524 if(ancLabel < 2) {//Colliding protons
1525 fhMCOrgMass[11]->Fill(pt,mass);
1526 fhMCOrgAsym[11]->Fill(pt,asym);
1527 fhMCOrgDeltaEta[11]->Fill(pt,deta);
1528 fhMCOrgDeltaPhi[11]->Fill(pt,dphi);
1529 }//colliding protons
1530 else if(ancLabel < 6){//partonic initial states interactions
1531 fhMCOrgMass[9]->Fill(pt,mass);
1532 fhMCOrgAsym[9]->Fill(pt,asym);
1533 fhMCOrgDeltaEta[9]->Fill(pt,deta);
1534 fhMCOrgDeltaPhi[9]->Fill(pt,dphi);
1535 }
1536 else if(ancLabel < 8){//Final state partons radiations?
1537 fhMCOrgMass[10]->Fill(pt,mass);
1538 fhMCOrgAsym[10]->Fill(pt,asym);
1539 fhMCOrgDeltaEta[10]->Fill(pt,deta);
1540 fhMCOrgDeltaPhi[10]->Fill(pt,dphi);
1541 }
4512d3e8 1542 // else {
1543 // printf("AliAnaPi0::FillMCVersusRecDataHistograms() - Check ** Common ancestor label %d, pdg %d, name %s, status %d; \n",
1544 // ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus);
1545 // }
6175da48 1546 }//status 21
4512d3e8 1547 //else {
1548 // printf("AliAnaPi0::FillMCVersusRecDataHistograms() - Check *** Common ancestor label %d, pdg %d, name %s, status %d; \n",
1549 // ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus);
1550 // }
6175da48 1551 }////Partons, colliding protons, strings, intermediate corrections
1552 }//ancLabel > -1
1553 else { //ancLabel <= -1
1554 //printf("Not related at all label = %d\n",ancLabel);
1555 fhMCOrgMass[12]->Fill(pt,mass);
1556 fhMCOrgAsym[12]->Fill(pt,asym);
1557 fhMCOrgDeltaEta[12]->Fill(pt,deta);
1558 fhMCOrgDeltaPhi[12]->Fill(pt,dphi);
1559 }
1560}
1561
20218aea 1562//____________________________________________________________________________________________________________________________________________________
1563void AliAnaPi0::CountAndGetAverages(Int_t &nClus,Int_t &nCell, Float_t &eClusTot,Float_t &eCellTot, Float_t &eDenClus,Float_t &eDenCell) {
1564// Count the number of clusters and cells, deposited energy, and do some averages in case multiplicity bins dependent on such numbers
1565// are requested
1566 if(fCalorimeter=="EMCAL"){
1567 nClus = GetEMCALClusters() ->GetEntriesFast();
1568 nCell = GetEMCALCells()->GetNumberOfCells();
1569 for(Int_t icl=0; icl < nClus; icl++) {
1570 Float_t e1 = ((AliVCluster*)GetEMCALClusters()->At(icl))->E();
1571 eClusTot += e1;
1572 // if(e1 > emax) emax = e1;
1573 // ((AliVCluster*)GetEMCALClusters()->At(icl))->GetPosition(pos1);
1574 // for(Int_t icl2=icl+1; icl2 < nClus; icl2++) {
1575 // Float_t e2 = ((AliVCluster*)GetEMCALClusters()->At(icl2))->E();
1576 // ((AliVCluster*)GetEMCALClusters()->At(icl2))->GetPosition(pos2);
1577 // rtmp = TMath::Sqrt((pos1[0]-pos2[0])*(pos1[0]-pos2[0]) + (pos1[2]-pos2[2])*(pos1[2]-pos2[2]));
1578 // 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);
1579 // rxz += rtmp;
1580 // rxzw += rtmpw;
1581 // ncomb++;
1582 // fhClusterPairDist ->Fill(rtmp);
1583 // fhClusterPairDistWeight->Fill(rtmpw);
1584 // //printf("Distance: %f; weighted %f\n ",rtmp,rtmp/(e1+((AliVCluster*)GetEMCALClusters()->At(icl2))->E()));
1585 //
1586 // }// second cluster loop
1587 }// first cluster
1588
1589 for(Int_t jce=0; jce < nCell; jce++) eCellTot += GetEMCALCells()->GetAmplitude(jce);
1590 }
1591 else {
1592 nClus = GetPHOSClusters()->GetEntriesFast();
1593 nCell = GetPHOSCells() ->GetNumberOfCells();
1594 for(Int_t icl=0; icl < nClus; icl++) {
1595 Float_t e1 = ((AliVCluster*)GetPHOSClusters()->At(icl))->E();
1596 eClusTot += e1;
1597 // ((AliVCluster*)GetPHOSClusters()->At(icl))->GetPosition(pos1);
1598 // for(Int_t icl2=icl+1; icl2 < nClus; icl2++) {
1599 // Float_t e2 = ((AliVCluster*)GetPHOSClusters()->At(icl2))->E();
1600 // ((AliVCluster*)GetPHOSClusters()->At(icl2))->GetPosition(pos2);
1601 // rtmp = TMath::Sqrt((pos1[0]-pos2[0])*(pos1[0]-pos2[0]) + (pos1[2]-pos2[2])*(pos1[2]-pos2[2]));
1602 // 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);
1603 // rxz += rtmp;
1604 // rxzw += rtmpw;
1605 // ncomb++;
1606 // fhClusterPairDist ->Fill(rtmp);
1607 // fhClusterPairDistWeight->Fill(rtmpw);
1608 // }// second cluster loop
1609 }// first cluster
1610 for(Int_t jce=0; jce < nCell; jce++) eCellTot += GetPHOSCells()->GetAmplitude(jce);
1611 }
1612 if(GetDebug() > 1)
1613 printf("AliAnaPi0::MakeAnalysisFillHistograms() - # Clusters %d, sum cluster E per SM %f,# Cells %d, sum cell E per SM %f\n", nClus,eClusTot,nCell,eCellTot);
1614
1615 //Fill histograms with "energy density", ncell and nclust will be > 0 since there are at least 2 "photons"
1616 eDenClus = eClusTot/nClus;
1617 eDenCell = eCellTot/nCell;
1618 fhEDensityCluster ->Fill(eDenClus);
1619 fhEDensityCell ->Fill(eDenCell);
1620 fhEDensityCellvsCluster->Fill(eDenClus, eDenCell);
1621 //Fill the average number of cells or clusters per SM
1622 eClusTot /=fNModules;
1623 eCellTot /=fNModules;
1624 fhAverTotECluster ->Fill(eClusTot);
1625 fhAverTotECell ->Fill(eCellTot);
1626 fhAverTotECellvsCluster->Fill(eClusTot, eCellTot);
1627 //printf("Average Cluster: E %f, density %f; Average Cell E %f, density %f\n ",eClusTot,eDenClus,eCellTot,eDenCell);
1628
1629 // //Average weighted pair distance
1630 // rxz /= ncomb;
1631 // rxzw /= ncomb;
1632 //
1633 // fhAverClusterPairDist ->Fill(rxz );
1634 // fhAverClusterPairDistWeight ->Fill(rxzw);
1635 // fhAverClusterPairDistvsAverE ->Fill(rxz ,eDenClus);
1636 // fhAverClusterPairDistWeightvsAverE->Fill(rxzw,eDenClus);
1637 // fhAverClusterPairDistvsN ->Fill(rxz ,nClus);
1638 // fhAverClusterPairDistWeightvsN ->Fill(rxzw,nClus);
1639 //
1640 // //emax
1641 // fhMaxEvsClustEDen->Fill(emax,eDenClus);
1642 // fhMaxEvsClustMult->Fill(emax,nPhot);
1643
1644 //printf("Average Distance: %f; weighted %f\n ",rxz,rxzw);
1645
1646}
1647
1c5acb87 1648//____________________________________________________________________________________________________________________________________________________
6639984f 1649void AliAnaPi0::MakeAnalysisFillHistograms()
1c5acb87 1650{
477d6cee 1651 //Process one event and extract photons from AOD branch
1652 // filled with AliAnaPhoton and fill histos with invariant mass
1653
6175da48 1654 //In case of simulated data, fill acceptance histograms
1655 if(IsDataMC())FillAcceptanceHistograms();
08a56f5f 1656
1657 //if (GetReader()->GetEventNumber()%10000 == 0)
1658 // printf("--- Event %d ---\n",GetReader()->GetEventNumber());
1659
6175da48 1660 //Init some variables
1661//Int_t iRun = (GetReader()->GetInputEvent())->GetRunNumber() ;
1662 Int_t nPhot = GetInputAODBranch()->GetEntriesFast() ;
1663 Int_t nClus = 0;
1664 Int_t nCell = 0;
1665 Float_t eClusTot = 0;
1666 Float_t eCellTot = 0;
156549ae 1667 Float_t eDenClus = 0;
1668 Float_t eDenCell = 0;
1669// Int_t ncomb = 0;
1670// Float_t rtmp = 0;
1671// Float_t rtmpw = 0;
1672// Float_t rxz = 0;
1673// Float_t rxzw = 0;
1674// Float_t pos1[3];
1675// Float_t pos2[3];
1676// Float_t emax = 0;
477d6cee 1677
72542aba 1678 if(GetNCentrBin() > 1 && (fUseAverCellEBins||fUseAverClusterEBins||fUseAverClusterEDenBins))
20218aea 1679 CountAndGetAverages(nClus,nCell,eClusTot,eCellTot,eDenClus,eDenCell);
1680
1681
156549ae 1682 if(GetDebug() > 1)
1683 printf("AliAnaPi0::MakeAnalysisFillHistograms() - Photon entries %d\n", nPhot);
1684
1685 //If less than photon 2 entries in the list, skip this event
20218aea 1686 if(nPhot < 2 ) {
156549ae 1687
20218aea 1688 if(GetDebug() > 2)
1689 printf("AliAnaPi0::MakeAnalysisFillHistograms() - nPhotons %d, cent bin %d continue to next event\n",nPhot, GetEventCentrality());
1690
72542aba 1691 if(GetNCentrBin() > 1) fhCentralityNoPair->Fill(GetEventCentrality() * GetNCentrBin() / GetReader()->GetCentralityOpt());
20218aea 1692
1693 return ;
6175da48 1694 }
6175da48 1695
1696 //Init variables
1697 Int_t module1 = -1;
1698 Int_t module2 = -1;
1699 Double_t vert[] = {0.0, 0.0, 0.0} ; //vertex
1700 Int_t evtIndex1 = 0 ;
1701 Int_t currentEvtIndex = -1;
1702 Int_t curCentrBin = 0 ;
1703 Int_t curRPBin = 0 ;
1704 Int_t curZvertBin = 0 ;
1705
c4a7d28a 1706 //Get shower shape information of clusters
1707 TObjArray *clusters = 0;
1708 if (fCalorimeter="EMCAL") clusters = GetEMCALClusters();
1709 else if(fCalorimeter="PHOS" ) clusters = GetPHOSClusters() ;
1710
6175da48 1711 //---------------------------------
1712 //First loop on photons/clusters
1713 //---------------------------------
477d6cee 1714 for(Int_t i1=0; i1<nPhot-1; i1++){
1715 AliAODPWG4Particle * p1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i1)) ;
6175da48 1716 //printf("AliAnaPi0::MakeAnalysisFillHistograms() : cluster1 id %d\n",p1->GetCaloLabel(0));
1717
7e7694bb 1718 // get the event index in the mixed buffer where the photon comes from
1719 // in case of mixing with analysis frame, not own mixing
c8fe2783 1720 evtIndex1 = GetEventIndex(p1, vert) ;
5025c139 1721 //printf("charge = %d\n", track->Charge());
c8fe2783 1722 if ( evtIndex1 == -1 )
1723 return ;
1724 if ( evtIndex1 == -2 )
1725 continue ;
41121cfe 1726
1727 //printf("z vertex %f < %f\n",vert[2],GetZvertexCut());
2244659d 1728 if(TMath::Abs(vert[2]) > GetZvertexCut()) continue ; //vertex cut
41121cfe 1729
6175da48 1730
1731 //----------------------------------------------------------------------------
1732 // Get the multiplicity bin. Different cases: centrality (PbPb),
1733 // average cluster multiplicity, average cell multiplicity, track multiplicity
1734 // default is centrality bins
1735 //----------------------------------------------------------------------------
c8fe2783 1736 if (evtIndex1 != currentEvtIndex) {
6175da48 1737 if(fUseTrackMultBins){ // Track multiplicity bins
1738 //printf("track mult %d\n",GetTrackMultiplicity());
1739 curCentrBin = (GetTrackMultiplicity()-1)/5;
72542aba 1740 if(curCentrBin > GetNCentrBin()-1) curCentrBin=GetNCentrBin()-1;
6175da48 1741 //printf("track mult bin %d\n",curCentrBin);
1742 }
1743 else if(fUsePhotonMultBins){ // Photon multiplicity bins
1744 //printf("photon mult %d cluster mult %d\n",nPhot, nClus);
72542aba 1745 curCentrBin = nPhot-2;
1746 if(curCentrBin > GetNCentrBin() -1) curCentrBin=GetNCentrBin()-1;
156549ae 1747 //printf("photon mult bin %d\n",curRPBin);
6175da48 1748 }
156549ae 1749 else if(fUseAverClusterEBins){ // Cluster average energy bins
6175da48 1750 //Bins for pp, if needed can be done in a more general way
72542aba 1751 curCentrBin = (Int_t) eClusTot/10 * GetNCentrBin();
1752 if(curCentrBin > GetNCentrBin()-1) curCentrBin=GetNCentrBin()-1;
6175da48 1753 //printf("cluster E average %f, bin %d \n",eClusTot,curCentrBin);
1754 }
156549ae 1755 else if(fUseAverCellEBins){ // Cell average energy bins
6175da48 1756 //Bins for pp, if needed can be done in a more general way
72542aba 1757 curCentrBin = (Int_t) eCellTot/10*GetNCentrBin();
1758 if(curCentrBin > GetNCentrBin()-1) curCentrBin=GetNCentrBin()-1;
6175da48 1759 //printf("cell E average %f, bin %d \n",eCellTot,curCentrBin);
1760 }
156549ae 1761 else if(fUseAverClusterEDenBins){ // Energy density bins
1762 //Bins for pp, if needed can be done in a more general way
72542aba 1763 curCentrBin = (Int_t) eDenClus/10*GetNCentrBin();
1764 if(curCentrBin > GetNCentrBin()-1) curCentrBin=GetNCentrBin()-1;
156549ae 1765 //printf("cluster Eden average %f, bin %d \n",eDenClus,curCentrBin);
1766 }
1767// else if(fUseAverClusterPairRBins){ // Cluster average distance bins
1768// //Bins for pp, if needed can be done in a more general way
72542aba 1769// curCentrBin = rxz/650*GetNCentrBin();
1770// if(curCentrBin > GetNCentrBin()-1) curCentrBin=GetNCentrBin()-1;
156549ae 1771// //printf("cluster pair R average %f, bin %d \n",rxz,curCentrBin);
1772// }
1773// else if(fUseAverClusterPairRWeightBins){ // Cluster average distance bins
1774// //Bins for pp, if needed can be done in a more general way
72542aba 1775// curCentrBin = rxzw/350*GetNCentrBin();
1776// if(curCentrBin > GetNCentrBin()-1) curCentrBin=GetNCentrBin()-1;
156549ae 1777// //printf("cluster pair rW average %f, bin %d \n",rxzw,curCentrBin);
1778// }
1779// else if(fUseEMaxBins){ // Cluster average distance bins
1780// //Bins for pp, if needed can be done in a more general way
72542aba 1781// curCentrBin = emax/20*GetNCentrBin();
1782// if(curCentrBin > GetNCentrBin()-1) curCentrBin=GetNCentrBin()-1;
156549ae 1783// //printf("cluster pair rW average %f, bin %d \n",rxzw,curCentrBin);
1784// }
6175da48 1785 else { //Event centrality
20218aea 1786 // Centrality task returns at maximum 10, 20 or 100, depending on option chosen and
1787 // number of bins, the bin has to be corrected
72542aba 1788 curCentrBin = GetEventCentrality() * GetNCentrBin() / GetReader()->GetCentralityOpt();
20218aea 1789 if(GetDebug() > 0 )printf("AliAnaPi0::MakeAnalysisFillHistograms() - curCentrBin %d, centrality %d, n bins %d, max bin from centrality %d\n",
72542aba 1790 curCentrBin, GetEventCentrality(), GetNCentrBin(), GetReader()->GetCentralityOpt());
6175da48 1791 }
6eb99ccd 1792
72542aba 1793 if (curCentrBin < 0 || curCentrBin >= GetNCentrBin()){
aba0ecc2 1794 if(GetDebug() > 0)
72542aba 1795 printf("AliAnaPi0::MakeAnalysisFillHistograms() - Centrality bin <%d> not expected, n bins <%d> , return\n",curCentrBin,GetNCentrBin());
6eb99ccd 1796 return;
1797 }
1798
156549ae 1799 //Reaction plane bin
c8fe2783 1800 curRPBin = 0 ;
606cbcb1 1801 if(GetNRPBin()>1 && GetEventPlane()){
72542aba 1802 Float_t epAngle = GetEventPlane()->GetEventplane(GetEventPlaneMethod());
1803 fhEventPlaneAngle->Fill(epAngle);
1804 curRPBin = TMath::Nint(epAngle*(GetNRPBin()-1)/TMath::Pi());
1805 if(curRPBin >= GetNRPBin()) printf("RP Bin %d out of range %d\n",curRPBin,GetNRPBin());
1806 //printf("RP: %d, %f, angle %f, n bin %d\n", curRPBin,epAngle*(GetNRPBin()-1)/TMath::Pi(),epAngle,GetNRPBin());
1807 }
1808
156549ae 1809 //Get vertex z bin
5025c139 1810 curZvertBin = (Int_t)(0.5*GetNZvertBin()*(vert[2]+GetZvertexCut())/GetZvertexCut()) ;
6175da48 1811
1812 //Fill event bin info
745f04da 1813 fhEvents ->Fill(curCentrBin+0.5,curZvertBin+0.5,curRPBin+0.5) ;
72542aba 1814 if(GetNCentrBin() > 1) {
1815 fhCentrality->Fill(curCentrBin);
606cbcb1 1816 if(GetNRPBin() > 1 && GetEventPlane()) fhEventPlaneResolution->Fill(curCentrBin,TMath::Cos(2.*GetEventPlane()->GetQsubRes()));
72542aba 1817 }
c8fe2783 1818 currentEvtIndex = evtIndex1 ;
ca468d44 1819 if(GetDebug() > 1)
6175da48 1820 printf("AliAnaPi0::MakeAnalysisFillHistograms() - Centrality %d, Vertex Bin %d, RP bin %d \n",curCentrBin,curRPBin,curZvertBin);
c8fe2783 1821 }
7e7694bb 1822
f8006433 1823 //printf("AliAnaPi0::MakeAnalysisFillHistograms(): Photon 1 Evt %d Vertex : %f,%f,%f\n",evtIndex1, GetVertex(evtIndex1)[0] ,GetVertex(evtIndex1)[1],GetVertex(evtIndex1)[2]);
af7b3903 1824
6175da48 1825 //Get the momentum of this cluster
477d6cee 1826 TLorentzVector photon1(p1->Px(),p1->Py(),p1->Pz(),p1->E());
6175da48 1827
1828 //Get (Super)Module number of this cluster
59b6bd99 1829 module1 = GetModuleNumber(p1);
6175da48 1830
c4a7d28a 1831
1832 //Get original cluster time,
1833 AliVCluster *cluster1 = 0;
1834 Bool_t bFound1 = kFALSE;
9ab9e937 1835 Int_t caloLabel1 = p1->GetCaloLabel(0);
1836 Bool_t iclus1 =-1;
1837 if(clusters){
1838 for(Int_t iclus = 0; iclus < clusters->GetEntriesFast(); iclus++){
1839 AliVCluster *cluster= dynamic_cast<AliVCluster*> (clusters->At(iclus));
1840 if(cluster){
1841 if (cluster->GetID()==caloLabel1) {
1842 bFound1 = kTRUE ;
1843 cluster1 = cluster;
1844 iclus1 = iclus;
1845 }
1846 }
1847 if(bFound1) break;
1848 }
c4a7d28a 1849 }// calorimeter clusters loop
1850
6175da48 1851 //---------------------------------
1852 //Second loop on photons/clusters
1853 //---------------------------------
477d6cee 1854 for(Int_t i2=i1+1; i2<nPhot; i2++){
1855 AliAODPWG4Particle * p2 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i2)) ;
c4a7d28a 1856
6175da48 1857 //In case of mixing frame, check we are not in the same event as the first cluster
c8fe2783 1858 Int_t evtIndex2 = GetEventIndex(p2, vert) ;
1859 if ( evtIndex2 == -1 )
1860 return ;
1861 if ( evtIndex2 == -2 )
1862 continue ;
1863 if (GetMixedEvent() && (evtIndex1 == evtIndex2))
7e7694bb 1864 continue ;
c4a7d28a 1865
9ab9e937 1866 //------------------------------------------
c4a7d28a 1867 //Check if time of clusters is similar
1868 AliVCluster *cluster2 = 0;
1869 Bool_t bFound2 = kFALSE;
1870 Int_t caloLabel2 = p2->GetCaloLabel(0);
9ab9e937 1871 if(clusters){
1872 for(Int_t iclus = iclus1+1; iclus < clusters->GetEntriesFast(); iclus++){
1873 AliVCluster *cluster= dynamic_cast<AliVCluster*> (clusters->At(iclus));
1874 if(cluster){
1875 if(cluster->GetID()==caloLabel2) {
1876 bFound2 = kTRUE ;
1877 cluster2 = cluster;
1878 }
1879 }
1880 if(bFound2) break;
1881 }// calorimeter clusters loop
1882 }
c4a7d28a 1883
1884 Float_t tof1 = -1;
1885 if(cluster1 && bFound1){
1886 tof1 = cluster1->GetTOF()*1e9;
1887 }
9ab9e937 1888// else printf("cluster1 not available: calo label %d / %d, cluster ID %d\n",
1889// p1->GetCaloLabel(0),(GetReader()->GetInputEvent())->GetNumberOfCaloClusters()-1,cluster1->GetID());
c4a7d28a 1890
1891 Float_t tof2 = -1;
1892 if(cluster2 && bFound2){
1893 tof2 = cluster2->GetTOF()*1e9;
1894 }
9ab9e937 1895// else printf("cluster2 not available: calo label %d / %d, cluster ID %d\n",
1896// p2->GetCaloLabel(0),(GetReader()->GetInputEvent())->GetNumberOfCaloClusters()-1,cluster2->GetID());
6175da48 1897
9ab9e937 1898 if(clusters){
1899 Double_t t12diff = tof1-tof2;
1900 if(TMath::Abs(t12diff) > GetPairTimeCut()) continue;
1901 }
1902 //------------------------------------------
1903
f8006433 1904 //printf("AliAnaPi0::MakeAnalysisFillHistograms(): Photon 2 Evt %d Vertex : %f,%f,%f\n",evtIndex2, GetVertex(evtIndex2)[0] ,GetVertex(evtIndex2)[1],GetVertex(evtIndex2)[2]);
9ab9e937 1905
6175da48 1906 //Get the momentum of this cluster
477d6cee 1907 TLorentzVector photon2(p2->Px(),p2->Py(),p2->Pz(),p2->E());
59b6bd99 1908 //Get module number
6175da48 1909 module2 = GetModuleNumber(p2);
1910
1911 //---------------------------------
1912 // Get pair kinematics
1913 //---------------------------------
1914 Double_t m = (photon1 + photon2).M() ;
1915 Double_t pt = (photon1 + photon2).Pt();
1916 Double_t deta = photon1.Eta() - photon2.Eta();
1917 Double_t dphi = photon1.Phi() - photon2.Phi();
1918 Double_t a = TMath::Abs(p1->E()-p2->E())/(p1->E()+p2->E()) ;
1919
477d6cee 1920 if(GetDebug() > 2)
6175da48 1921 printf(" E: photon1 %f, photon2 %f; Pair: pT %f, mass %f, a %f\n", p1->E(), p2->E(), (photon1 + photon2).E(),m,a);
1922
1923 //--------------------------------
1924 // Opening angle selection
1925 //--------------------------------
50f39b97 1926 //Check if opening angle is too large or too small compared to what is expected
1927 Double_t angle = photon1.Angle(photon2.Vect());
6175da48 1928 if(fUseAngleEDepCut && !GetNeutralMesonSelection()->IsAngleInWindow((photon1+photon2).E(),angle+0.05)) {
1929 if(GetDebug() > 2)
1930 printf("AliAnaPi0::MakeAnalysisFillHistograms() -Real pair angle %f not in E %f window\n",angle, (photon1+photon2).E());
c8fe2783 1931 continue;
6175da48 1932 }
af7b3903 1933
6175da48 1934 if(fUseAngleCut && (angle < fAngleCut || angle > fAngleMaxCut)) {
1935 if(GetDebug() > 2)
1936 printf("AliAnaPi0::MakeAnalysisFillHistograms() - Real pair cut %f < angle %f < cut %f\n",fAngleCut, angle, fAngleMaxCut);
1937 continue;
1938 }
1939
1940 //-------------------------------------------------------------------------------------------------
af7b3903 1941 //Fill module dependent histograms, put a cut on assymmetry on the first available cut in the array
6175da48 1942 //-------------------------------------------------------------------------------------------------
20218aea 1943 if(a < fAsymCuts[0] && fFillSMCombinations){
af7b3903 1944 if(module1==module2 && module1 >=0 && module1<fNModules)
91e1ea12 1945 fhReMod[module1]->Fill(pt,m) ;
8d230fa8 1946
af7b3903 1947 if(fCalorimeter=="EMCAL"){
8d230fa8 1948
1949 // Same sector
1950 Int_t j=0;
1951 for(Int_t i = 0; i < fNModules/2; i++){
1952 j=2*i;
91e1ea12 1953 if((module1==j && module2==j+1) || (module1==j+1 && module2==j)) fhReSameSectorEMCALMod[i]->Fill(pt,m) ;
8d230fa8 1954 }
1955
1956 // Same side
1957 for(Int_t i = 0; i < fNModules-2; i++){
91e1ea12 1958 if((module1==i && module2==i+2) || (module1==i+2 && module2==i)) fhReSameSideEMCALMod[i]->Fill(pt,m);
8d230fa8 1959 }
1960 }//EMCAL
1961 else {//PHOS
91e1ea12 1962 if((module1==0 && module2==1) || (module1==1 && module2==0)) fhReDiffPHOSMod[0]->Fill(pt,m) ;
1963 if((module1==0 && module2==2) || (module1==2 && module2==0)) fhReDiffPHOSMod[1]->Fill(pt,m) ;
1964 if((module1==1 && module2==2) || (module1==2 && module2==1)) fhReDiffPHOSMod[2]->Fill(pt,m) ;
8d230fa8 1965 }//PHOS
821c8090 1966 }
7e7694bb 1967
af7b3903 1968 //In case we want only pairs in same (super) module, check their origin.
1969 Bool_t ok = kTRUE;
1970 if(fSameSM && module1!=module2) ok=kFALSE;
1971 if(ok){
6175da48 1972
1973 //Check if one of the clusters comes from a conversion
d39cba7e 1974 if(fCheckConversion){
1975 if (p1->IsTagged() && p2->IsTagged()) fhReConv2->Fill(pt,m);
1976 else if(p1->IsTagged() || p2->IsTagged()) fhReConv ->Fill(pt,m);
1977 }
6175da48 1978
af7b3903 1979 //Fill histograms for different bad channel distance, centrality, assymmetry cut and pid bit
5ae09196 1980 for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
af7b3903 1981 if((p1->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton)) && (p2->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton))){
1982 for(Int_t iasym=0; iasym < fNAsymCuts; iasym++){
1983 if(a < fAsymCuts[iasym]){
1984 Int_t index = ((curCentrBin*fNPIDBits)+ipid)*fNAsymCuts + iasym;
6175da48 1985 //printf("index %d :(cen %d * nPID %d + ipid %d)*nasym %d + iasym %d\n",index,curCentrBin,fNPIDBits,ipid,fNAsymCuts,iasym);
91e1ea12 1986 fhRe1 [index]->Fill(pt,m);
1987 if(fMakeInvPtPlots)fhReInvPt1[index]->Fill(pt,m,1./pt) ;
6175da48 1988 if(fFillBadDistHisto){
1989 if(p1->DistToBad()>0 && p2->DistToBad()>0){
91e1ea12 1990 fhRe2 [index]->Fill(pt,m) ;
1991 if(fMakeInvPtPlots)fhReInvPt2[index]->Fill(pt,m,1./pt) ;
6175da48 1992 if(p1->DistToBad()>1 && p2->DistToBad()>1){
91e1ea12 1993 fhRe3 [index]->Fill(pt,m) ;
1994 if(fMakeInvPtPlots)fhReInvPt3[index]->Fill(pt,m,1./pt) ;
6175da48 1995 }// bad 3
1996 }// bad2
1997 }// Fill bad dist histos
1998 }//assymetry cut
1999 }// asymmetry cut loop
af7b3903 2000 }// bad 1
2001 }// pid bit loop
5ae09196 2002
af7b3903 2003 //Fill histograms with opening angle
91e1ea12 2004 fhRealOpeningAngle ->Fill(pt,angle);
2005 fhRealCosOpeningAngle->Fill(pt,TMath::Cos(angle));
af7b3903 2006
2007 //Fill histograms with pair assymmetry
91e1ea12 2008 fhRePtAsym->Fill(pt,a);
2009 if(m > 0.10 && m < 0.17) fhRePtAsymPi0->Fill(pt,a);
2010 if(m > 0.45 && m < 0.65) fhRePtAsymEta->Fill(pt,a);
af7b3903 2011
6175da48 2012 //-------------------------------------------------------
2013 //Get the number of cells needed for multi cut analysis.
2014 //-------------------------------------------------------
2015 Int_t ncell1 = 0;
2016 Int_t ncell2 = 0;
2017 if(fMultiCutAna || (IsDataMC() && fMultiCutAnaSim)){
2018
af7b3903 2019 AliVEvent * event = GetReader()->GetInputEvent();
2020 if(event){
2021 for(Int_t iclus = 0; iclus < event->GetNumberOfCaloClusters(); iclus++){
2022 AliVCluster *cluster = event->GetCaloCluster(iclus);
5ae09196 2023
af7b3903 2024 Bool_t is = kFALSE;
2025 if (fCalorimeter == "EMCAL" && GetReader()->IsEMCALCluster(cluster)) is = kTRUE;
2026 else if(fCalorimeter == "PHOS" && GetReader()->IsPHOSCluster (cluster)) is = kTRUE;
5ae09196 2027
af7b3903 2028 if(is){
2029 if (p1->GetCaloLabel(0) == cluster->GetID()) ncell1 = cluster->GetNCells();
2030 else if (p2->GetCaloLabel(0) == cluster->GetID()) ncell2 = cluster->GetNCells();
2031 } // PHOS or EMCAL cluster as requested in analysis
2032
2033 if(ncell2 > 0 && ncell1 > 0) break; // No need to continue the iteration
2034
2035 }
2036 //printf("e 1: %2.2f, e 2: %2.2f, ncells: n1 %d, n2 %d\n", p1->E(), p2->E(),ncell1,ncell2);
2037 }
6175da48 2038 }
2039
2040 //---------
2041 // MC data
2042 //---------
2043 //Do some MC checks on the origin of the pair, is there any common ancestor and if there is one, who?
2044 if(IsDataMC()) FillMCVersusRecDataHistograms(p1->GetLabel(), p2->GetLabel(),p1->Pt(), p2->Pt(),ncell1, ncell2, m, pt, a,deta, dphi);
2045
2046 //-----------------------
2047 //Multi cuts analysis
2048 //-----------------------
2049 if(fMultiCutAna){
2050 //Histograms for different PID bits selection
2051 for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
2052
2053 if(p1->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton) &&
91e1ea12 2054 p2->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton)) fhRePIDBits[ipid]->Fill(pt,m) ;
6175da48 2055
2056 //printf("ipt %d, ipid%d, name %s\n",ipt, ipid, fhRePtPIDCuts[ipt*fNPIDBitsBits+ipid]->GetName());
2057 } // pid bit cut loop
2058
2059 //Several pt,ncell and asymmetry cuts
af7b3903 2060 for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
2061 for(Int_t icell=0; icell<fNCellNCuts; icell++){
2062 for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
2063 Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
ba1eeb1f 2064 if(p1->E() > fPtCuts[ipt] && p2->E() > fPtCuts[ipt] &&
af7b3903 2065 a < fAsymCuts[iasym] &&
6175da48 2066 ncell1 >= fCellNCuts[icell] && ncell2 >= fCellNCuts[icell]){
91e1ea12 2067 fhRePtNCellAsymCuts[index]->Fill(pt,m) ;
6175da48 2068 //printf("ipt %d, icell%d, iasym %d, name %s\n",ipt, icell, iasym, fhRePtNCellAsymCuts[((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym]->GetName());
9302260a 2069 if(fFillSMCombinations && module1==module2){
91e1ea12 2070 if (module1==0) fhRePtNCellAsymCutsSM0[index]->Fill(pt,m) ;
2071 else if(module1==1) fhRePtNCellAsymCutsSM1[index]->Fill(pt,m) ;
2072 else if(module1==2) fhRePtNCellAsymCutsSM2[index]->Fill(pt,m) ;
2073 else if(module1==3) fhRePtNCellAsymCutsSM3[index]->Fill(pt,m) ;
20218aea 2074 //else printf("AliAnaPi0::FillHistograms() - WRONG SM NUMBER\n");
6175da48 2075 }
2076 }
af7b3903 2077 }// pid bit cut loop
2078 }// icell loop
2079 }// pt cut loop
91e1ea12 2080 if(GetHistoTrackMultiplicityBins()){
2081 for(Int_t iasym = 0; iasym < fNAsymCuts; iasym++){
2082 if(a < fAsymCuts[iasym])fhRePtMult[iasym]->Fill(pt,GetTrackMultiplicity(),m) ;
2083 }
2084 }
af7b3903 2085 }// multiple cuts analysis
2086 }// ok if same sm
7e7694bb 2087 }// second same event particle
2088 }// first cluster
6175da48 2089
2090 //-------------------------------------------------------------
2091 // Mixing
2092 //-------------------------------------------------------------
7e7694bb 2093 if(fDoOwnMix){
156549ae 2094 //printf("Cen bin %d, RP bin %d, e aver %f, mult %d\n",curCentrBin,curRPBin, eClusTot, nPhot);
6175da48 2095 //Recover events in with same characteristics as the current event
5025c139 2096 TList * evMixList=fEventsList[curCentrBin*GetNZvertBin()*GetNRPBin()+curZvertBin*GetNRPBin()+curRPBin] ;
7e7694bb 2097 Int_t nMixed = evMixList->GetSize() ;
2098 for(Int_t ii=0; ii<nMixed; ii++){
2099 TClonesArray* ev2= (TClonesArray*) (evMixList->At(ii));
2100 Int_t nPhot2=ev2->GetEntriesFast() ;
2101 Double_t m = -999;
6175da48 2102 if(GetDebug() > 1)
2103 printf("AliAnaPi0::MakeAnalysisFillHistograms() - Mixed event %d photon entries %d, centrality bin %d\n", ii, nPhot2, curCentrBin);
7e7694bb 2104
6175da48 2105 //---------------------------------
2106 //First loop on photons/clusters
2107 //---------------------------------
7e7694bb 2108 for(Int_t i1=0; i1<nPhot; i1++){
2109 AliAODPWG4Particle * p1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i1)) ;
6175da48 2110 if(fSameSM && GetModuleNumber(p1)!=module1) continue;
2111
2112 //Get kinematics of cluster and (super) module of this cluster
7e7694bb 2113 TLorentzVector photon1(p1->Px(),p1->Py(),p1->Pz(),p1->E());
af7b3903 2114 module1 = GetModuleNumber(p1);
91e1ea12 2115
6175da48 2116 //---------------------------------
2117 //First loop on photons/clusters
2118 //---------------------------------
7e7694bb 2119 for(Int_t i2=0; i2<nPhot2; i2++){
2120 AliAODPWG4Particle * p2 = (AliAODPWG4Particle*) (ev2->At(i2)) ;
91e1ea12 2121
6175da48 2122 //Get kinematics of second cluster and calculate those of the pair
7e7694bb 2123 TLorentzVector photon2(p2->Px(),p2->Py(),p2->Pz(),p2->E());
6175da48 2124 m = (photon1+photon2).M() ;
7e7694bb 2125 Double_t pt = (photon1 + photon2).Pt();
2126 Double_t a = TMath::Abs(p1->E()-p2->E())/(p1->E()+p2->E()) ;
2127
2128 //Check if opening angle is too large or too small compared to what is expected
2129 Double_t angle = photon1.Angle(photon2.Vect());
6175da48 2130 if(fUseAngleEDepCut && !GetNeutralMesonSelection()->IsAngleInWindow((photon1+photon2).E(),angle+0.05)){
2131 if(GetDebug() > 2)
2132 printf("AliAnaPi0::MakeAnalysisFillHistograms() -Mix pair angle %f not in E %f window\n",angle, (photon1+photon2).E());
2133 continue;
2134 }
2135 if(fUseAngleCut && (angle < fAngleCut || angle > fAngleMaxCut)) {
2136 if(GetDebug() > 2)
2137 printf("AliAnaPi0::MakeAnalysisFillHistograms() -Mix pair angle %f < cut %f\n",angle,fAngleCut);
2138 continue;
2139
2140 }
7e7694bb 2141
2142 if(GetDebug() > 2)
2143 printf("AliAnaPi0::MakeAnalysisFillHistograms() - Mixed Event: pT: photon1 %2.2f, photon2 %2.2f; Pair: pT %2.2f, mass %2.3f, a %f2.3\n",
af7b3903 2144 p1->Pt(), p2->Pt(), pt,m,a);
6175da48 2145
af7b3903 2146 //In case we want only pairs in same (super) module, check their origin.
2147 module2 = GetModuleNumber(p2);
6175da48 2148
2149 //-------------------------------------------------------------------------------------------------
2150 //Fill module dependent histograms, put a cut on assymmetry on the first available cut in the array
2151 //-------------------------------------------------------------------------------------------------
20218aea 2152 if(a < fAsymCuts[0] && fFillSMCombinations){
6175da48 2153 if(module1==module2 && module1 >=0 && module1<fNModules)
91e1ea12 2154 fhMiMod[module1]->Fill(pt,m) ;
6175da48 2155
8d230fa8 2156 if(fCalorimeter=="EMCAL"){
2157
2158 // Same sector
2159 Int_t j=0;
2160 for(Int_t i = 0; i < fNModules/2; i++){
2161 j=2*i;
91e1ea12 2162 if((module1==j && module2==j+1) || (module1==j+1 && module2==j)) fhMiSameSectorEMCALMod[i]->Fill(pt,m) ;
8d230fa8 2163 }
2164
2165 // Same side
2166 for(Int_t i = 0; i < fNModules-2; i++){
91e1ea12 2167 if((module1==i && module2==i+2) || (module1==i+2 && module2==i)) fhMiSameSideEMCALMod[i]->Fill(pt,m);
8d230fa8 2168 }
2169 }//EMCAL
2170 else {//PHOS
91e1ea12 2171 if((module1==0 && module2==1) || (module1==1 && module2==0)) fhMiDiffPHOSMod[0]->Fill(pt,m) ;
2172 if((module1==0 && module2==2) || (module1==2 && module2==0)) fhMiDiffPHOSMod[1]->Fill(pt,m) ;
2173 if((module1==1 && module2==2) || (module1==2 && module2==1)) fhMiDiffPHOSMod[2]->Fill(pt,m) ;
8d230fa8 2174 }//PHOS
2175
2176
6175da48 2177 }
2178
af7b3903 2179 Bool_t ok = kTRUE;
2180 if(fSameSM && module1!=module2) ok=kFALSE;
2181 if(ok){
6175da48 2182
2183 //Check if one of the clusters comes from a conversion
d39cba7e 2184 if(fCheckConversion){
2185 if (p1->IsTagged() && p2->IsTagged()) fhMiConv2->Fill(pt,m);
2186 else if(p1->IsTagged() || p2->IsTagged()) fhMiConv ->Fill(pt,m);
2187 }
6175da48 2188 //Fill histograms for different bad channel distance, centrality, assymmetry cut and pid bit
af7b3903 2189 for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
2190 if((p1->IsPIDOK(ipid,AliCaloPID::kPhoton)) && (p2->IsPIDOK(ipid,AliCaloPID::kPhoton))){
2191 for(Int_t iasym=0; iasym < fNAsymCuts; iasym++){
2192 if(a < fAsymCuts[iasym]){
2193 Int_t index = ((curCentrBin*fNPIDBits)+ipid)*fNAsymCuts + iasym;
91e1ea12 2194 fhMi1 [index]->Fill(pt,m) ;
2195 if(fMakeInvPtPlots)fhMiInvPt1[index]->Fill(pt,m,1./pt) ;
6175da48 2196 if(fFillBadDistHisto){
2197 if(p1->DistToBad()>0 && p2->DistToBad()>0){
2198 fhMi2 [index]->Fill(pt,m) ;
91e1ea12 2199 if(fMakeInvPtPlots)fhMiInvPt2[index]->Fill(pt,m,1./pt) ;
6175da48 2200 if(p1->DistToBad()>1 && p2->DistToBad()>1){
2201 fhMi3 [index]->Fill(pt,m) ;
91e1ea12 2202 if(fMakeInvPtPlots)fhMiInvPt3[index]->Fill(pt,m,1./pt) ;
6175da48 2203 }
af7b3903 2204 }
6175da48 2205 }// Fill bad dist histo
af7b3903 2206 }//Asymmetry cut
2207 }// Asymmetry loop
2208 }//PID cut
2209 }//loop for histograms
6175da48 2210
2211 //-----------------------
2212 //Multi cuts analysis
2213 //-----------------------
2214 if(fMultiCutAna){
2215 //Several pt,ncell and asymmetry cuts
2216 for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
2217 for(Int_t icell=0; icell<fNCellNCuts; icell++){
2218 for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
2219 Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
2220 if(p1->Pt() > fPtCuts[ipt] && p2->Pt() > fPtCuts[ipt] &&
2221 a < fAsymCuts[iasym] &&
2222 p1->GetBtag() >= fCellNCuts[icell] && p2->GetBtag() >= fCellNCuts[icell]){
91e1ea12 2223 fhMiPtNCellAsymCuts[index]->Fill(pt,m) ;
6175da48 2224 //printf("ipt %d, icell%d, iasym %d, name %s\n",ipt, icell, iasym, fhRePtNCellAsymCuts[((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym]->GetName());
2225 }
2226 }// pid bit cut loop
2227 }// icell loop
2228 }// pt cut loop
2229 } // Multi cut ana
2230
2231 //Fill histograms with opening angle
91e1ea12 2232 fhMixedOpeningAngle ->Fill(pt,angle);
2233 fhMixedCosOpeningAngle->Fill(pt,TMath::Cos(angle));
af7b3903 2234 }//ok
7e7694bb 2235 }// second cluster loop
2236 }//first cluster loop
2237 }//loop on mixed events
2238
6175da48 2239 //--------------------------------------------------------
2240 //Add the current event to the list of events for mixing
2241 //--------------------------------------------------------
7e7694bb 2242 TClonesArray *currentEvent = new TClonesArray(*GetInputAODBranch());
af7b3903 2243 //Add current event to buffer and Remove redundant events
7e7694bb 2244 if(currentEvent->GetEntriesFast()>0){
2245 evMixList->AddFirst(currentEvent) ;
2246 currentEvent=0 ; //Now list of particles belongs to buffer and it will be deleted with buffer
72542aba 2247 if(evMixList->GetSize() >= GetNMaxEvMix())
7e7694bb 2248 {
2249 TClonesArray * tmp = (TClonesArray*) (evMixList->Last()) ;
2250 evMixList->RemoveLast() ;
2251 delete tmp ;
2252 }
2253 }
2254 else{ //empty event
2255 delete currentEvent ;
2256 currentEvent=0 ;
477d6cee 2257 }
7e7694bb 2258 }// DoOwnMix
c8fe2783 2259
1c5acb87 2260}
2261
a5cc4f03 2262//________________________________________________________________________
2263void AliAnaPi0::ReadHistograms(TList* outputList)
2264{
50f39b97 2265 // Needed when Terminate is executed in distributed environment
2266 // Refill analysis histograms of this class with corresponding histograms in output list.
2267
2268 // Histograms of this analsys are kept in the same list as other analysis, recover the position of
2269 // the first one and then add the next.
2270 Int_t index = outputList->IndexOf(outputList->FindObject(GetAddedHistogramsStringToName()+"hRe_cen0_pid0_dist1"));
2271
72542aba 2272 if(!fhRe1) fhRe1 = new TH2D*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
2273 if(!fhRe2) fhRe2 = new TH2D*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
2274 if(!fhRe3) fhRe3 = new TH2D*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
2275 if(!fhMi1) fhMi1 = new TH2D*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
2276 if(!fhMi2) fhMi2 = new TH2D*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
2277 if(!fhMi3) fhMi3 = new TH2D*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
2278 if(!fhReInvPt1) fhReInvPt1 = new TH2D*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
2279 if(!fhReInvPt2) fhReInvPt2 = new TH2D*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
2280 if(!fhReInvPt3) fhReInvPt3 = new TH2D*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
2281 if(!fhMiInvPt1) fhMiInvPt1 = new TH2D*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
2282 if(!fhMiInvPt2) fhMiInvPt2 = new TH2D*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
2283 if(!fhMiInvPt3) fhMiInvPt3 = new TH2D*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
20218aea 2284
2285 if(fFillSMCombinations){
2286 if(!fhReMod) fhReMod = new TH2D*[fNModules] ;
2287 if(!fhReDiffPHOSMod) fhReDiffPHOSMod = new TH2D*[fNModules] ;
2288 if(!fhReSameSectorEMCALMod)fhReSameSectorEMCALMod = new TH2D*[fNModules/2] ;
2289 if(!fhReSameSideEMCALMod) fhReSameSideEMCALMod = new TH2D*[fNModules-2] ;
2290 if(!fhMiMod) fhMiMod = new TH2D*[fNModules] ;
2291 if(!fhMiDiffPHOSMod) fhMiDiffPHOSMod = new TH2D*[fNModules] ;
2292 if(!fhMiSameSectorEMCALMod)fhMiSameSectorEMCALMod = new TH2D*[fNModules/2] ;
2293 if(!fhMiSameSideEMCALMod) fhMiSameSideEMCALMod = new TH2D*[fNModules-2] ;
2294 }
2295
2296 if(fCheckConversion){
2297 fhReConv = (TH2D*) outputList->At(index++);
2298 fhMiConv = (TH2D*) outputList->At(index++);
2299 fhReConv2 = (TH2D*) outputList->At(index++);
2300 fhMiConv2 = (TH2D*) outputList->At(index++);
2301 }
2302
72542aba 2303 for(Int_t ic=0; ic<GetNCentrBin(); ic++){
af7b3903 2304 for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
2305 for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
2306 Int_t ihisto = ((ic*fNPIDBits)+ipid)*fNAsymCuts + iasym;
2307
2308 fhRe1[ihisto] = (TH2D*) outputList->At(index++);
2309 fhRe2[ihisto] = (TH2D*) outputList->At(index++);
2310 fhRe3[ihisto] = (TH2D*) outputList->At(index++);
6175da48 2311
2312 fhReInvPt1[ihisto] = (TH2D*) outputList->At(index++);
2313 fhReInvPt2[ihisto] = (TH2D*) outputList->At(index++);
2314 fhReInvPt3[ihisto] = (TH2D*) outputList->At(index++);
5ae09196 2315
af7b3903 2316 if(fDoOwnMix){
2317 fhMi1[ihisto] = (TH2D*) outputList->At(index++);
2318 fhMi2[ihisto] = (TH2D*) outputList->At(index++);
2319 fhMi3[ihisto] = (TH2D*) outputList->At(index++);
6175da48 2320
2321 fhMiInvPt1[ihisto] = (TH2D*) outputList->At(index++);
2322 fhMiInvPt2[ihisto] = (TH2D*) outputList->At(index++);
2323 fhMiInvPt3[ihisto] = (TH2D*) outputList->At(index++);
af7b3903 2324 }//Own mix
2325 }//asymmetry loop
2326 }// pid loop
2327 }// centrality loop
2328
2329 fhRePtAsym = (TH2D*)outputList->At(index++);
2330 fhRePtAsymPi0 = (TH2D*)outputList->At(index++);
2331 fhRePtAsymEta = (TH2D*)outputList->At(index++);
eee5fcf1 2332
5ae09196 2333 if(fMultiCutAna){
2334
eee5fcf1 2335 if(!fhRePtNCellAsymCuts) fhRePtNCellAsymCuts = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
2336 if(!fhRePIDBits) fhRePIDBits = new TH2D*[fNPIDBits];
2337
5ae09196 2338 for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
2339 fhRePIDBits[ipid] = (TH2D*) outputList->At(index++);
2340 }// ipid loop
2341
2342 for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
2343 for(Int_t icell=0; icell<fNCellNCuts; icell++){
2344 for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
2345 fhRePtNCellAsymCuts[((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym] = (TH2D*) outputList->At(index++);
2346 }// iasym
2347 }// icell loop
2348 }// ipt loop
af7b3903 2349
2350 if(!fhRePtMult) fhRePtMult = new TH3D*[fNAsymCuts] ;
2351 for(Int_t iasym = 0; iasym < fNAsymCuts; iasym++)
2352 fhRePtMult[iasym] = (TH3D*) outputList->At(index++);
5ae09196 2353 }// multi cut analysis
50f39b97 2354
745f04da 2355 fhEvents = (TH3D *) outputList->At(index++);
72542aba 2356 if(GetNCentrBin()>1)fhCentrality = (TH1D *) outputList->At(index++);
2357 if(GetNCentrBin()>1)fhCentralityNoPair = (TH1D *) outputList->At(index++);
2358 if(GetNRPBin()>1) fhEventPlaneAngle = (TH1D *) outputList->At(index++);
2359 if(GetNRPBin()>1 && GetNCentrBin()>1)fhEventPlaneResolution = (TH2D *) outputList->At(index++);
745f04da 2360
af7b3903 2361 fhRealOpeningAngle = (TH2D*) outputList->At(index++);
2362 fhRealCosOpeningAngle = (TH2D*) outputList->At(index++);
6175da48 2363 if(fDoOwnMix){
2364 fhMixedOpeningAngle = (TH2D*) outputList->At(index++);
2365 fhMixedCosOpeningAngle = (TH2D*) outputList->At(index++);
2366 }
af7b3903 2367
50f39b97 2368 //Histograms filled only if MC data is requested
2369 if(IsDataMC() || (GetReader()->GetDataType() == AliCaloTrackReader::kMC) ){
6175da48 2370 fhPrimPi0Pt = (TH1D*) outputList->At(index++);
2371 fhPrimPi0AccPt = (TH1D*) outputList->At(index++);
156549ae 2372 fhPrimPi0Y = (TH2D*) outputList->At(index++);
2373 fhPrimPi0AccY = (TH2D*) outputList->At(index++);
2374 fhPrimPi0Phi = (TH2D*) outputList->At(index++);
2375 fhPrimPi0AccPhi = (TH2D*) outputList->At(index++);
2376 fhPrimEtaPt = (TH1D*) outputList->At(index++);
2377 fhPrimEtaAccPt = (TH1D*) outputList->At(index++);
2378 fhPrimEtaY = (TH2D*) outputList->At(index++);
2379 fhPrimEtaAccY = (TH2D*) outputList->At(index++);
2380 fhPrimEtaPhi = (TH2D*) outputList->At(index++);
2381 fhPrimEtaAccPhi = (TH2D*) outputList->At(index++);
6175da48 2382 for(Int_t i = 0; i<13; i++){
2383 fhMCOrgMass[i] = (TH2D*) outputList->At(index++);
2384 fhMCOrgAsym[i] = (TH2D*) outputList->At(index++);
2385 fhMCOrgDeltaEta[i] = (TH2D*) outputList->At(index++);
2386 fhMCOrgDeltaPhi[i] = (TH2D*) outputList->At(index++);
2387 }
2388
2389 if(fMultiCutAnaSim){
2390 fhMCPi0MassPtTrue = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
2391 fhMCPi0MassPtRec = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
2392 fhMCPi0PtTruePtRec = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
2393 fhMCEtaMassPtTrue = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
2394 fhMCEtaMassPtRec = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
2395 fhMCEtaPtTruePtRec = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
2396 for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
2397 for(Int_t icell=0; icell<fNCellNCuts; icell++){
2398 for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
2399 Int_t in = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
2400 fhMCPi0MassPtTrue[in] = (TH2D*) outputList->At(index++);
2401 fhMCPi0PtTruePtRec[in] = (TH2D*) outputList->At(index++);
2402 fhMCEtaMassPtTrue[in] = (TH2D*) outputList->At(index++);
2403 fhMCEtaPtTruePtRec[in] = (TH2D*) outputList->At(index++);
2404 }
2405 }
2406 }
2407 }
2408 else{
2409 fhMCPi0MassPtTrue = new TH2D*[1];
2410 fhMCPi0PtTruePtRec = new TH2D*[1];
2411 fhMCEtaMassPtTrue = new TH2D*[1];
2412 fhMCEtaPtTruePtRec = new TH2D*[1];
2413
2414 fhMCPi0MassPtTrue[0] = (TH2D*) outputList->At(index++);
2415 fhMCPi0PtTruePtRec[0] = (TH2D*) outputList->At(index++);
2416 fhMCEtaMassPtTrue[0] = (TH2D*) outputList->At(index++);
2417 fhMCEtaPtTruePtRec[0] = (TH2D*) outputList->At(index++);
2418 }
50f39b97 2419 }
2420
6175da48 2421 for(Int_t imod=0; imod < fNModules; imod++){
8d230fa8 2422 fhReMod[imod] = (TH2D*) outputList->At(index++);
2423 if(fCalorimeter=="EMCAL"){
2424 if(imod < fNModules/2) fhReSameSectorEMCALMod[imod] = (TH2D*) outputList->At(index++);
2425 if(imod < fNModules-2) fhReSameSideEMCALMod[imod] = (TH2D*) outputList->At(index++);
2426 }
2427 else fhReDiffPHOSMod[imod] = (TH2D*) outputList->At(index++);
2428
6175da48 2429 if(fDoOwnMix){
8d230fa8 2430 fhMiMod[imod] = (TH2D*) outputList->At(index++);
2431 if(fCalorimeter=="EMCAL"){
2432 if(imod < fNModules/2) fhMiSameSectorEMCALMod[imod] = (TH2D*) outputList->At(index++);
2433 if(imod < fNModules-2) fhMiSameSideEMCALMod[imod] = (TH2D*) outputList->At(index++);
2434 }
2435 else fhMiDiffPHOSMod[imod] = (TH2D*) outputList->At(index++);
6175da48 2436 }
2437 }
eee5fcf1 2438
a5cc4f03 2439}
2440
2441
6639984f 2442//____________________________________________________________________________________________________________________________________________________
a5cc4f03 2443void AliAnaPi0::Terminate(TList* outputList)
6639984f 2444{
2445 //Do some calculations and plots from the final histograms.
477d6cee 2446
fbeaf916 2447 printf(" *** %s Terminate:\n", GetName()) ;
50f39b97 2448
a5cc4f03 2449 //Recover histograms from output histograms list, needed for distributed analysis.
2450 ReadHistograms(outputList);
50f39b97 2451
2e557d1c 2452 if (!fhRe1) {
50f39b97 2453 printf("AliAnaPi0::Terminate() - Error: Remote output histograms not imported in AliAnaPi0 object");
2454 return;
2e557d1c 2455 }
50f39b97 2456
a3aebfff 2457 printf("AliAnaPi0::Terminate() Mgg Real : %5.3f , RMS : %5.3f \n", fhRe1[0]->GetMean(), fhRe1[0]->GetRMS() ) ;
5ae09196 2458
2459 const Int_t buffersize = 255;
2460
2461 char nameIM[buffersize];
2462 snprintf(nameIM, buffersize,"AliAnaPi0_%s_cPt",fCalorimeter.Data());
71dd883b 2463 TCanvas * cIM = new TCanvas(nameIM, "", 400, 10, 600, 700) ;
6639984f 2464 cIM->Divide(2, 2);
50f39b97 2465
6639984f 2466 cIM->cd(1) ;
2467 //gPad->SetLogy();
af7b3903 2468 TH1D * hIMAllPt = (TH1D*) fhRe1[0]->ProjectionY(Form("IMPtAll_%s",fCalorimeter.Data()));
6639984f 2469 hIMAllPt->SetLineColor(2);
2470 hIMAllPt->SetTitle("No cut on p_{T, #gamma#gamma} ");
2471 hIMAllPt->Draw();
2472
2473 cIM->cd(2) ;
af7b3903 2474 TH1D * hIMPt5 = (TH1D*) fhRe1[0]->ProjectionY(Form("IMPt0-5_%s",fCalorimeter.Data()),0, fhRe1[0]->GetXaxis()->FindBin(5.));
2244659d 2475// hRe1Pt5->GetXaxis()->SetRangeUser(0,5);
2476// TH1D * hIMPt5 = (TH1D*) hRe1Pt5->Project3D(Form("IMPt5_%s_pz",fCalorimeter.Data()));
6639984f 2477 hIMPt5->SetLineColor(2);
2478 hIMPt5->SetTitle("0 < p_{T, #gamma#gamma} < 5 GeV/c");
2479 hIMPt5->Draw();
2480
2481 cIM->cd(3) ;
af7b3903 2482 TH1D * hIMPt10 = (TH1D*) fhRe1[0]->ProjectionY(Form("IMPt5-10_%s",fCalorimeter.Data()), fhRe1[0]->GetXaxis()->FindBin(5.),fhRe1[0]->GetXaxis()->FindBin(10.));
2244659d 2483// hRe1Pt10->GetXaxis()->SetRangeUser(5,10);
2484// TH1D * hIMPt10 = (TH1D*) hRe1Pt10->Project3D(Form("IMPt10_%s_pz",fCalorimeter.Data()));
6639984f 2485 hIMPt10->SetLineColor(2);
2486 hIMPt10->SetTitle("5 < p_{T, #gamma#gamma} < 10 GeV/c");
2487 hIMPt10->Draw();
2488
2489 cIM->cd(4) ;
af7b3903 2490 TH1D * hIMPt20 = (TH1D*) fhRe1[0]->ProjectionY(Form("IMPt10-20_%s",fCalorimeter.Data()), fhRe1[0]->GetXaxis()->FindBin(10.),fhRe1[0]->GetXaxis()->FindBin(20.));
2244659d 2491 // TH3F * hRe1Pt20 = (TH3F*)fhRe1[0]->Clone(Form("IMPt20_%s",fCalorimeter.Data()));
2492// hRe1Pt20->GetXaxis()->SetRangeUser(10,20);
2493// TH1D * hIMPt20 = (TH1D*) hRe1Pt20->Project3D(Form("IMPt20_%s_pz",fCalorimeter.Data()));
6639984f 2494 hIMPt20->SetLineColor(2);
2495 hIMPt20->SetTitle("10 < p_{T, #gamma#gamma} < 20 GeV/c");
2496 hIMPt20->Draw();
2497
5ae09196 2498 char nameIMF[buffersize];
2499 snprintf(nameIMF,buffersize,"AliAnaPi0_%s_Mgg.eps",fCalorimeter.Data());
71dd883b 2500 cIM->Print(nameIMF);
6639984f 2501
5ae09196 2502 char namePt[buffersize];
2503 snprintf(namePt,buffersize,"AliAnaPi0_%s_cPt",fCalorimeter.Data());
71dd883b 2504 TCanvas * cPt = new TCanvas(namePt, "", 400, 10, 600, 700) ;
6639984f 2505 cPt->Divide(2, 2);
2506
2507 cPt->cd(1) ;
2508 //gPad->SetLogy();
af7b3903 2509 TH1D * hPt = (TH1D*) fhRe1[0]->ProjectionX(Form("Pt0_%s",fCalorimeter.Data()),-1,-1);
6639984f 2510 hPt->SetLineColor(2);
2511 hPt->SetTitle("No cut on M_{#gamma#gamma} ");
2512 hPt->Draw();
2513
2514 cPt->cd(2) ;
af7b3903 2515 TH1D * hPtIM1 = (TH1D*)fhRe1[0]->ProjectionX(Form("Pt1_%s",fCalorimeter.Data()), fhRe1[0]->GetZaxis()->FindBin(0.05),fhRe1[0]->GetZaxis()->FindBin(0.21));
2244659d 2516// TH3F * hRe1IM1 = (TH3F*)fhRe1[0]->Clone(Form("Pt1_%s",fCalorimeter.Data()));
2517// hRe1IM1->GetZaxis()->SetRangeUser(0.05,0.21);
2518// TH1D * hPtIM1 = (TH1D*) hRe1IM1->Project3D("x");
6639984f 2519 hPtIM1->SetLineColor(2);
2520 hPtIM1->SetTitle("0.05 < M_{#gamma#gamma} < 0.21 GeV/c^{2}");
2521 hPtIM1->Draw();
2522
2523 cPt->cd(3) ;
af7b3903 2524 TH1D * hPtIM2 = (TH1D*)fhRe1[0]->ProjectionX(Form("Pt2_%s",fCalorimeter.Data()), fhRe1[0]->GetZaxis()->FindBin(0.09),fhRe1[0]->GetZaxis()->FindBin(0.17));
2244659d 2525// TH3F * hRe1IM2 = (TH3F*)fhRe1[0]->Clone(Form("Pt2_%s",fCalorimeter.Data()));
2526// hRe1IM2->GetZaxis()->SetRangeUser(0.09,0.17);
2527// TH1D * hPtIM2 = (TH1D*) hRe1IM2->Project3D("x");
6639984f 2528 hPtIM2->SetLineColor(2);
2529 hPtIM2->SetTitle("0.09 < M_{#gamma#gamma} < 0.17 GeV/c^{2}");
2530 hPtIM2->Draw();
2531
2532 cPt->cd(4) ;
af7b3903 2533 TH1D * hPtIM3 = (TH1D*)fhRe1[0]->ProjectionX(Form("Pt3_%s",fCalorimeter.Data()), fhRe1[0]->GetZaxis()->FindBin(0.11),fhRe1[0]->GetZaxis()->FindBin(0.15));
2244659d 2534// TH3F * hRe1IM3 = (TH3F*)fhRe1[0]->Clone(Form("Pt3_%s",fCalorimeter.Data()));
2535// hRe1IM3->GetZaxis()->SetRangeUser(0.11,0.15);
2536// TH1D * hPtIM3 = (TH1D*) hRe1IM1->Project3D("x");
6639984f 2537 hPtIM3->SetLineColor(2);
2538 hPtIM3->SetTitle("0.11 < M_{#gamma#gamma} < 0.15 GeV/c^{2}");
2539 hPtIM3->Draw();
2540
164a1d84 2541 char namePtF[buffersize];
5ae09196 2542 snprintf(namePtF,buffersize,"AliAnaPi0_%s_Pt.eps",fCalorimeter.Data());
71dd883b 2543 cPt->Print(namePtF);
1c5acb87 2544
5ae09196 2545 char line[buffersize] ;
2546 snprintf(line,buffersize,".!tar -zcf %s_%s.tar.gz *.eps", GetName(),fCalorimeter.Data()) ;
6639984f 2547 gROOT->ProcessLine(line);
5ae09196 2548 snprintf(line, buffersize,".!rm -fR AliAnaPi0_%s*.eps",fCalorimeter.Data());
6639984f 2549 gROOT->ProcessLine(line);
2550
71dd883b 2551 printf(" AliAnaPi0::Terminate() - !! All the eps files are in %s_%s.tar.gz !!!\n", GetName(), fCalorimeter.Data());
1c5acb87 2552
6639984f 2553}
c8fe2783 2554 //____________________________________________________________________________________________________________________________________________________
2555Int_t AliAnaPi0::GetEventIndex(AliAODPWG4Particle * part, Double_t * vert)
2556{
f8006433 2557 // retieves the event index and checks the vertex
2558 // in the mixed buffer returns -2 if vertex NOK
2559 // for normal events returns 0 if vertex OK and -1 if vertex NOK
2560
2561 Int_t evtIndex = -1 ;
2562 if(GetReader()->GetDataType()!=AliCaloTrackReader::kMC){
2563
2564 if (GetMixedEvent()){
2565
2566 evtIndex = GetMixedEvent()->EventIndexForCaloCluster(part->GetCaloLabel(0)) ;
2567 GetVertex(vert,evtIndex);
2568
5025c139 2569 if(TMath::Abs(vert[2])> GetZvertexCut())
f8006433 2570 evtIndex = -2 ; //Event can not be used (vertex, centrality,... cuts not fulfilled)
2571 } else {// Single event
2572
2573 GetVertex(vert);
2574
5025c139 2575 if(TMath::Abs(vert[2])> GetZvertexCut())
f8006433 2576 evtIndex = -1 ; //Event can not be used (vertex, centrality,... cuts not fulfilled)
2577 else
2578 evtIndex = 0 ;
c8fe2783 2579 }
0ae57829 2580 }//No MC reader
f8006433 2581 else {
2582 evtIndex = 0;
2583 vert[0] = 0. ;
2584 vert[1] = 0. ;
2585 vert[2] = 0. ;
2586 }
0ae57829 2587
f8006433 2588 return evtIndex ;
c8fe2783 2589}
f8006433 2590