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