AliPhysics  vAN-20150810 (e477c3b)
 All Classes Namespaces Files Functions Variables Enumerations Enumerator Macros
AliAnaPi0.cxx
Go to the documentation of this file.
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 
16 // --- ROOT system ---
17 #include "TH3.h"
18 #include "TH2F.h"
19 //#include "Riostream.h"
20 #include "TCanvas.h"
21 #include "TPad.h"
22 #include "TROOT.h"
23 #include "TClonesArray.h"
24 #include "TObjString.h"
25 #include "TDatabasePDG.h"
26 
27 //---- AliRoot system ----
28 #include "AliAnaPi0.h"
29 #include "AliCaloTrackReader.h"
30 #include "AliCaloPID.h"
31 #include "AliStack.h"
32 #include "AliFiducialCut.h"
33 #include "TParticle.h"
34 #include "AliVEvent.h"
35 #include "AliESDCaloCluster.h"
36 #include "AliESDEvent.h"
37 #include "AliAODEvent.h"
38 #include "AliNeutralMesonSelection.h"
39 #include "AliMixedEvent.h"
40 #include "AliAODMCParticle.h"
41 
42 // --- Detectors ---
43 #include "AliPHOSGeoUtils.h"
44 #include "AliEMCALGeometry.h"
45 
47 ClassImp(AliAnaPi0) ;
49 
50 //______________________________________________________
52 //______________________________________________________
53 AliAnaPi0::AliAnaPi0() : AliAnaCaloTrackCorrBaseClass(),
54 fEventsList(0x0),
55 fNModules(22),
56 fUseAngleCut(kFALSE), fUseAngleEDepCut(kFALSE), fAngleCut(0), fAngleMaxCut(0.),
57 fMultiCutAna(kFALSE), fMultiCutAnaSim(kFALSE),
58 fNPtCuts(0), fNAsymCuts(0), fNCellNCuts(0), fNPIDBits(0),
59 fMakeInvPtPlots(kFALSE), fSameSM(kFALSE),
60 fFillSMCombinations(kFALSE), fCheckConversion(kFALSE),
61 fFillBadDistHisto(kFALSE), fFillSSCombinations(kFALSE),
62 fFillAngleHisto(kFALSE), fFillAsymmetryHisto(kFALSE), fFillOriginHisto(0),
63 fFillArmenterosThetaStar(0), fFillOnlyMCAcceptanceHisto(0),
64 fCheckAccInSector(kFALSE),
65 fPhotonMom1(), fPhotonMom1Boost(), fPhotonMom2(), fPi0Mom(),
66 fProdVertex(),
67 
68 // Histograms
69 fhAverTotECluster(0), fhAverTotECell(0), fhAverTotECellvsCluster(0),
70 fhEDensityCluster(0), fhEDensityCell(0), fhEDensityCellvsCluster(0),
71 fhReMod(0x0), fhReSameSideEMCALMod(0x0), fhReSameSectorEMCALMod(0x0), fhReDiffPHOSMod(0x0),
72 fhMiMod(0x0), fhMiSameSideEMCALMod(0x0), fhMiSameSectorEMCALMod(0x0), fhMiDiffPHOSMod(0x0),
73 fhReConv(0x0), fhMiConv(0x0), fhReConv2(0x0), fhMiConv2(0x0),
74 fhRe1(0x0), fhMi1(0x0), fhRe2(0x0), fhMi2(0x0),
75 fhRe3(0x0), fhMi3(0x0), fhReInvPt1(0x0), fhMiInvPt1(0x0),
76 fhReInvPt2(0x0), fhMiInvPt2(0x0), fhReInvPt3(0x0), fhMiInvPt3(0x0),
77 fhRePtNCellAsymCuts(0x0), fhMiPtNCellAsymCuts(0x0), fhRePtNCellAsymCutsSM(),
78 fhRePIDBits(0x0), fhRePtMult(0x0), fhReSS(),
79 fhRePtAsym(0x0), fhRePtAsymPi0(0x0), fhRePtAsymEta(0x0),
80 fhEventBin(0), fhEventMixBin(0),
81 fhCentrality(0x0), fhCentralityNoPair(0x0),
82 fhEventPlaneResolution(0x0),
83 fhRealOpeningAngle(0x0), fhRealCosOpeningAngle(0x0), fhMixedOpeningAngle(0x0), fhMixedCosOpeningAngle(0x0),
84 
85 // MC histograms
86 fhPrimPi0E(0x0), fhPrimPi0Pt(0x0),
87 fhPrimPi0AccE(0x0), fhPrimPi0AccPt(0x0),
88 fhPrimPi0Y(0x0), fhPrimPi0AccY(0x0),
89 fhPrimPi0Yeta(0x0), fhPrimPi0YetaYcut(0x0), fhPrimPi0AccYeta(0x0),
90 fhPrimPi0Phi(0x0), fhPrimPi0AccPhi(0x0),
91 fhPrimPi0OpeningAngle(0x0), fhPrimPi0OpeningAngleAsym(0x0),fhPrimPi0CosOpeningAngle(0x0),
92 fhPrimPi0PtCentrality(0), fhPrimPi0PtEventPlane(0),
93 fhPrimPi0AccPtCentrality(0), fhPrimPi0AccPtEventPlane(0),
94 fhPrimEtaE(0x0), fhPrimEtaPt(0x0),
95 fhPrimEtaAccE(0x0), fhPrimEtaAccPt(0x0),
96 fhPrimEtaY(0x0), fhPrimEtaAccY(0x0),
97 fhPrimEtaYeta(0x0), fhPrimEtaYetaYcut(0x0), fhPrimEtaAccYeta(0x0),
98 fhPrimEtaPhi(0x0), fhPrimEtaAccPhi(0x0),
99 fhPrimEtaOpeningAngle(0x0), fhPrimEtaOpeningAngleAsym(0x0),fhPrimEtaCosOpeningAngle(0x0),
100 fhPrimEtaPtCentrality(0), fhPrimEtaPtEventPlane(0),
101 fhPrimEtaAccPtCentrality(0), fhPrimEtaAccPtEventPlane(0),
102 fhPrimPi0PtOrigin(0x0), fhPrimEtaPtOrigin(0x0),
103 fhPrimNotResonancePi0PtOrigin(0x0), fhPrimPi0PtStatus(0x0),
104 fhMCOrgMass(), fhMCOrgAsym(), fhMCOrgDeltaEta(), fhMCOrgDeltaPhi(),
105 fhMCPi0MassPtRec(), fhMCPi0MassPtTrue(), fhMCPi0PtTruePtRec(),
106 fhMCEtaMassPtRec(), fhMCEtaMassPtTrue(), fhMCEtaPtTruePtRec(),
107 fhMCPi0PtOrigin(0x0), fhMCEtaPtOrigin(0x0),
108 fhMCNotResonancePi0PtOrigin(0),fhMCPi0PtStatus(0x0),
109 fhMCPi0ProdVertex(0), fhMCEtaProdVertex(0),
110 fhPrimPi0ProdVertex(0), fhPrimEtaProdVertex(0),
111 fhReMCFromConversion(0), fhReMCFromNotConversion(0), fhReMCFromMixConversion(0),
112 fhCosThStarPrimPi0(0), fhCosThStarPrimEta(0),
113 fhEPairDiffTime(0)
114 {
115  InitParameters();
116 
117  for(Int_t i = 0; i < 4; i++)
118  {
119  fhArmPrimEta[i] = 0;
120  fhArmPrimPi0[i] = 0;
121  }
122 }
123 
124 //_____________________
126 //_____________________
127 AliAnaPi0::~AliAnaPi0()
128 {
129  // Remove event containers
130 
131  if(DoOwnMix() && fEventsList)
132  {
133  for(Int_t ic=0; ic<GetNCentrBin(); ic++)
134  {
135  for(Int_t iz=0; iz<GetNZvertBin(); iz++)
136  {
137  for(Int_t irp=0; irp<GetNRPBin(); irp++)
138  {
139  Int_t bin = GetEventMixBin(ic,iz,irp);
140  fEventsList[bin]->Delete() ;
141  delete fEventsList[bin] ;
142  }
143  }
144  }
145  delete[] fEventsList;
146  }
147 }
148 
149 //______________________________
152 //______________________________
153 void AliAnaPi0::InitParameters()
154 {
155  SetInputAODName("PWG4Particle");
156 
157  AddToHistogramsName("AnaPi0_");
158 
159  fUseAngleEDepCut = kFALSE;
160 
161  fUseAngleCut = kTRUE;
162  fAngleCut = 0.;
163  fAngleMaxCut = DegToRad(100.); // 100 degrees cut, avoid EMCal/DCal combinations
164 
165  fMultiCutAna = kFALSE;
166 
167  fNPtCuts = 3;
168  fPtCuts[0] = 0.; fPtCuts[1] = 0.3; fPtCuts[2] = 0.5;
169  for(Int_t i = fNPtCuts; i < 10; i++)fPtCuts[i] = 0.;
170 
171  fNAsymCuts = 2;
172  fAsymCuts[0] = 1.; fAsymCuts[1] = 0.7; //fAsymCuts[2] = 0.6; // fAsymCuts[3] = 0.1;
173  for(Int_t i = fNAsymCuts; i < 10; i++)fAsymCuts[i] = 0.;
174 
175  fNCellNCuts = 3;
176  fCellNCuts[0] = 0; fCellNCuts[1] = 1; fCellNCuts[2] = 2;
177  for(Int_t i = fNCellNCuts; i < 10; i++)fCellNCuts[i] = 0;
178 
179  fNPIDBits = 2;
180  fPIDBits[0] = 0; fPIDBits[1] = 2; // fPIDBits[2] = 4; fPIDBits[3] = 6;// check, no cut, dispersion, neutral, dispersion&&neutral
181  for(Int_t i = fNPIDBits; i < 10; i++)fPIDBits[i] = 0;
182 }
183 
184 //_______________________________________
186 //_______________________________________
187 TObjString * AliAnaPi0::GetAnalysisCuts()
188 {
189  TString parList ; //this will be list of parameters used for this analysis.
190  const Int_t buffersize = 255;
191  char onePar[buffersize] ;
192  snprintf(onePar,buffersize,"--- AliAnaPi0 ---:") ;
193  parList+=onePar ;
194  snprintf(onePar,buffersize,"Number of bins in Centrality: %d;",GetNCentrBin()) ;
195  parList+=onePar ;
196  snprintf(onePar,buffersize,"Number of bins in Z vert. pos: %d;",GetNZvertBin()) ;
197  parList+=onePar ;
198  snprintf(onePar,buffersize,"Number of bins in Reac. Plain: %d;",GetNRPBin()) ;
199  parList+=onePar ;
200  snprintf(onePar,buffersize,"Depth of event buffer: %d;",GetNMaxEvMix()) ;
201  parList+=onePar ;
202  snprintf(onePar,buffersize,"Select pairs with their angle: %d, edep %d, min angle %2.3f, max angle %2.3f;",fUseAngleCut, fUseAngleEDepCut,fAngleCut,fAngleMaxCut) ;
203  parList+=onePar ;
204  snprintf(onePar,buffersize," Asymmetry cuts: n = %d, asymmetry < ",fNAsymCuts) ;
205  for(Int_t i = 0; i < fNAsymCuts; i++) snprintf(onePar,buffersize,"%s %2.2f;",onePar,fAsymCuts[i]);
206  parList+=onePar ;
207  snprintf(onePar,buffersize," PID selection bits: n = %d, PID bit =",fNPIDBits) ;
208  for(Int_t i = 0; i < fNPIDBits; i++) snprintf(onePar,buffersize,"%s %d;",onePar,fPIDBits[i]);
209  parList+=onePar ;
210  snprintf(onePar,buffersize,"Cuts:") ;
211  parList+=onePar ;
212  snprintf(onePar,buffersize,"Z vertex position: -%f < z < %f;",GetZvertexCut(),GetZvertexCut()) ;
213  parList+=onePar ;
214  snprintf(onePar,buffersize,"Calorimeter: %s;",GetCalorimeterString().Data()) ;
215  parList+=onePar ;
216  snprintf(onePar,buffersize,"Number of modules: %d:",fNModules) ;
217  parList+=onePar ;
218  if(fMultiCutAna)
219  {
220  snprintf(onePar, buffersize," pT cuts: n = %d, pt > ",fNPtCuts) ;
221  for(Int_t i = 0; i < fNPtCuts; i++) snprintf(onePar,buffersize,"%s %2.2f;",onePar,fPtCuts[i]);
222  parList+=onePar ;
223  snprintf(onePar,buffersize, " N cell in cluster cuts: n = %d, nCell > ",fNCellNCuts) ;
224  for(Int_t i = 0; i < fNCellNCuts; i++) snprintf(onePar,buffersize,"%s %d;",onePar,fCellNCuts[i]);
225  parList+=onePar ;
226  }
227 
228  return new TObjString(parList) ;
229 }
230 
231 //_________________________________________
234 //_________________________________________
235 TList * AliAnaPi0::GetCreateOutputObjects()
236 {
237  TList * outputContainer = new TList() ;
238  outputContainer->SetName(GetName());
239 
240  // Set sizes and ranges
241  const Int_t buffersize = 255;
242  char key[buffersize] ;
243  char title[buffersize] ;
244 
245  Int_t nptbins = GetHistogramRanges()->GetHistoPtBins();
246  Int_t nphibins = GetHistogramRanges()->GetHistoPhiBins();
247  Int_t netabins = GetHistogramRanges()->GetHistoEtaBins();
248  Float_t ptmax = GetHistogramRanges()->GetHistoPtMax();
249  Float_t phimax = GetHistogramRanges()->GetHistoPhiMax();
250  Float_t etamax = GetHistogramRanges()->GetHistoEtaMax();
251  Float_t ptmin = GetHistogramRanges()->GetHistoPtMin();
252  Float_t phimin = GetHistogramRanges()->GetHistoPhiMin();
253  Float_t etamin = GetHistogramRanges()->GetHistoEtaMin();
254 
255  Int_t nmassbins = GetHistogramRanges()->GetHistoMassBins();
256  Int_t nasymbins = GetHistogramRanges()->GetHistoAsymmetryBins();
257  Float_t massmax = GetHistogramRanges()->GetHistoMassMax();
258  Float_t asymmax = GetHistogramRanges()->GetHistoAsymmetryMax();
259  Float_t massmin = GetHistogramRanges()->GetHistoMassMin();
260  Float_t asymmin = GetHistogramRanges()->GetHistoAsymmetryMin();
261  Int_t ntrmbins = GetHistogramRanges()->GetHistoTrackMultiplicityBins();
262  Int_t ntrmmax = GetHistogramRanges()->GetHistoTrackMultiplicityMax();
263  Int_t ntrmmin = GetHistogramRanges()->GetHistoTrackMultiplicityMin();
264  Int_t tdbins = GetHistogramRanges()->GetHistoDiffTimeBins() ;
265  Float_t tdmax = GetHistogramRanges()->GetHistoDiffTimeMax();
266  Float_t tdmin = GetHistogramRanges()->GetHistoDiffTimeMin();
267 
268  // Start with pure MC kinematics histograms
269  // In case other tasks just need this info like AliAnaPi0EbE
270  if(IsDataMC())
271  {
272  // Pi0
273 
274  fhPrimPi0E = new TH1F("hPrimPi0E","Primary #pi^{0} E, |#it{Y}|<1",nptbins,ptmin,ptmax) ;
275  fhPrimPi0AccE = new TH1F("hPrimPi0AccE","Primary #pi^{0} #it{E} with both photons in acceptance",nptbins,ptmin,ptmax) ;
276  fhPrimPi0E ->SetXTitle("#it{E} (GeV)");
277  fhPrimPi0AccE->SetXTitle("#it{E} (GeV)");
278  outputContainer->Add(fhPrimPi0E) ;
279  outputContainer->Add(fhPrimPi0AccE) ;
280 
281  fhPrimPi0Pt = new TH1F("hPrimPi0Pt","Primary #pi^{0} #it{p}_{T} , |#it{Y}|<1",nptbins,ptmin,ptmax) ;
282  fhPrimPi0AccPt = new TH1F("hPrimPi0AccPt","Primary #pi^{0} #it{p}_{T} with both photons in acceptance",nptbins,ptmin,ptmax) ;
283  fhPrimPi0Pt ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
284  fhPrimPi0AccPt->SetXTitle("#it{p}_{T} (GeV/#it{c})");
285  outputContainer->Add(fhPrimPi0Pt) ;
286  outputContainer->Add(fhPrimPi0AccPt) ;
287 
288  Int_t netabinsopen = TMath::Nint(netabins*4/(etamax-etamin));
289  fhPrimPi0Y = new TH2F("hPrimPi0Rapidity","Rapidity of primary #pi^{0}",nptbins,ptmin,ptmax,netabinsopen,-2, 2) ;
290  fhPrimPi0Y ->SetYTitle("#it{Rapidity}");
291  fhPrimPi0Y ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
292  outputContainer->Add(fhPrimPi0Y) ;
293 
294  fhPrimPi0Yeta = new TH2F("hPrimPi0PseudoRapidity","PseudoRapidity of primary #pi^{0}",nptbins,ptmin,ptmax,netabinsopen,-2, 2) ;
295  fhPrimPi0Yeta ->SetYTitle("#eta");
296  fhPrimPi0Yeta ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
297  outputContainer->Add(fhPrimPi0Yeta) ;
298 
299  fhPrimPi0YetaYcut = new TH2F("hPrimPi0PseudoRapidityYcut","PseudoRapidity of primary #pi^{0}, |#it{Y}|<1",nptbins,ptmin,ptmax,netabinsopen,-2, 2) ;
300  fhPrimPi0YetaYcut ->SetYTitle("#eta");
301  fhPrimPi0YetaYcut ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
302  outputContainer->Add(fhPrimPi0YetaYcut) ;
303 
304  fhPrimPi0AccY = new TH2F("hPrimPi0AccRapidity","Rapidity of primary #pi^{0} with accepted daughters",nptbins,ptmin,ptmax,netabins,etamin,etamax) ;
305  fhPrimPi0AccY->SetYTitle("Rapidity");
306  fhPrimPi0AccY->SetXTitle("#it{p}_{T} (GeV/#it{c})");
307  outputContainer->Add(fhPrimPi0AccY) ;
308 
309  fhPrimPi0AccYeta = new TH2F("hPrimPi0AccPseudoRapidity","PseudoRapidity of primary #pi^{0} with accepted daughters",nptbins,ptmin,ptmax,netabins,etamin,etamax) ;
310  fhPrimPi0AccYeta ->SetYTitle("#eta");
311  fhPrimPi0AccYeta ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
312  outputContainer->Add(fhPrimPi0AccYeta) ;
313 
314  Int_t nphibinsopen = TMath::Nint(nphibins*TMath::TwoPi()/(phimax-phimin));
315  fhPrimPi0Phi = new TH2F("hPrimPi0Phi","#phi of primary #pi^{0}, |#it{Y}|<1",nptbins,ptmin,ptmax,nphibinsopen,0,360) ;
316  fhPrimPi0Phi->SetYTitle("#phi (deg)");
317  fhPrimPi0Phi->SetXTitle("#it{p}_{T} (GeV/#it{c})");
318  outputContainer->Add(fhPrimPi0Phi) ;
319 
320  fhPrimPi0AccPhi = new TH2F("hPrimPi0AccPhi","#phi of primary #pi^{0} with accepted daughters",nptbins,ptmin,ptmax,
321  nphibins,phimin*TMath::RadToDeg(),phimax*TMath::RadToDeg()) ;
322  fhPrimPi0AccPhi->SetYTitle("#phi (deg)");
323  fhPrimPi0AccPhi->SetXTitle("#it{p}_{T} (GeV/#it{c})");
324  outputContainer->Add(fhPrimPi0AccPhi) ;
325 
326  // Eta
327 
328  fhPrimEtaE = new TH1F("hPrimEtaE","Primary eta E",nptbins,ptmin,ptmax) ;
329  fhPrimEtaAccE = new TH1F("hPrimEtaAccE","Primary #eta #it{E} with both photons in acceptance",nptbins,ptmin,ptmax) ;
330  fhPrimEtaE ->SetXTitle("#it{E} (GeV)");
331  fhPrimEtaAccE->SetXTitle("#it{E} (GeV)");
332  outputContainer->Add(fhPrimEtaE) ;
333  outputContainer->Add(fhPrimEtaAccE) ;
334 
335  fhPrimEtaPt = new TH1F("hPrimEtaPt","Primary #eta #it{p}_{T}",nptbins,ptmin,ptmax) ;
336  fhPrimEtaAccPt = new TH1F("hPrimEtaAccPt","Primary eta #it{p}_{T} with both photons in acceptance",nptbins,ptmin,ptmax) ;
337  fhPrimEtaPt ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
338  fhPrimEtaAccPt->SetXTitle("#it{p}_{T} (GeV/#it{c})");
339  outputContainer->Add(fhPrimEtaPt) ;
340  outputContainer->Add(fhPrimEtaAccPt) ;
341 
342  fhPrimEtaY = new TH2F("hPrimEtaRapidity","Rapidity of primary #eta",nptbins,ptmin,ptmax,netabinsopen,-2, 2) ;
343  fhPrimEtaY->SetYTitle("#it{Rapidity}");
344  fhPrimEtaY->SetXTitle("#it{p}_{T} (GeV/#it{c})");
345  outputContainer->Add(fhPrimEtaY) ;
346 
347  fhPrimEtaYeta = new TH2F("hPrimEtaPseudoRapidityEta","PseudoRapidity of primary #eta",nptbins,ptmin,ptmax,netabinsopen,-2, 2) ;
348  fhPrimEtaYeta->SetYTitle("#it{Rapidity}");
349  fhPrimEtaYeta->SetXTitle("#it{p}_{T} (GeV/#it{c})");
350  outputContainer->Add(fhPrimEtaYeta) ;
351 
352  fhPrimEtaYetaYcut = new TH2F("hPrimEtaPseudoRapidityEtaYcut","PseudoRapidity of primary #eta, |#it{Y}|<1",nptbins,ptmin,ptmax,netabinsopen,-2, 2) ;
353  fhPrimEtaYetaYcut->SetYTitle("#it{Pseudorapidity}");
354  fhPrimEtaYetaYcut->SetXTitle("#it{p}_{T} (GeV/#it{c})");
355  outputContainer->Add(fhPrimEtaYetaYcut) ;
356 
357  fhPrimEtaAccY = new TH2F("hPrimEtaAccRapidity","Rapidity of primary #eta",nptbins,ptmin,ptmax, netabins,etamin,etamax) ;
358  fhPrimEtaAccY->SetYTitle("#it{Rapidity}");
359  fhPrimEtaAccY->SetXTitle("#it{p}_{T} (GeV/#it{c})");
360  outputContainer->Add(fhPrimEtaAccY) ;
361 
362  fhPrimEtaAccYeta = new TH2F("hPrimEtaAccPseudoRapidity","PseudoRapidity of primary #eta",nptbins,ptmin,ptmax, netabins,etamin,etamax) ;
363  fhPrimEtaAccYeta->SetYTitle("#it{Pseudorapidity}");
364  fhPrimEtaAccYeta->SetXTitle("#it{p}_{T} (GeV/#it{c})");
365  outputContainer->Add(fhPrimEtaAccYeta) ;
366 
367  fhPrimEtaPhi = new TH2F("hPrimEtaPhi","Azimuthal of primary #eta",nptbins,ptmin,ptmax, nphibinsopen,0,360) ;
368  fhPrimEtaPhi->SetYTitle("#phi (deg)");
369  fhPrimEtaPhi->SetXTitle("#it{p}_{T} (GeV/#it{c})");
370  outputContainer->Add(fhPrimEtaPhi) ;
371 
372  fhPrimEtaAccPhi = new TH2F("hPrimEtaAccPhi","Azimuthal of primary #eta with accepted daughters",nptbins,ptmin,ptmax, nphibins,phimin*TMath::RadToDeg(),phimax*TMath::RadToDeg()) ;
373  fhPrimEtaAccPhi->SetYTitle("#phi (deg)");
374  fhPrimEtaAccPhi->SetXTitle("#it{p}_{T} (GeV/#it{c})");
375  outputContainer->Add(fhPrimEtaAccPhi) ;
376 
377  // Create histograms only for PbPb or high multiplicity analysis analysis
378  if( IsHighMultiplicityAnalysisOn() )
379  {
380  fhPrimPi0PtCentrality = new TH2F("hPrimPi0PtCentrality","Primary #pi^{0} #it{p}_{T} vs reco centrality, |#it{Y}|<1",
381  nptbins,ptmin,ptmax, 100, 0, 100) ;
382  fhPrimPi0AccPtCentrality = new TH2F("hPrimPi0AccPtCentrality","Primary #pi^{0} with both photons in acceptance #it{p}_{T} vs reco centrality",
383  nptbins,ptmin,ptmax, 100, 0, 100) ;
384  fhPrimPi0PtCentrality ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
385  fhPrimPi0AccPtCentrality->SetXTitle("#it{p}_{T} (GeV/#it{c})");
386  fhPrimPi0PtCentrality ->SetYTitle("Centrality");
387  fhPrimPi0AccPtCentrality->SetYTitle("Centrality");
388  outputContainer->Add(fhPrimPi0PtCentrality) ;
389  outputContainer->Add(fhPrimPi0AccPtCentrality) ;
390 
391  fhPrimEtaPtCentrality = new TH2F("hPrimEtaPtCentrality","Primary #eta #it{p}_{T} vs reco centrality, |#it{Y}|<1",nptbins,ptmin,ptmax, 100, 0, 100) ;
392  fhPrimEtaAccPtCentrality = new TH2F("hPrimEtaAccPtCentrality","Primary #eta with both photons in acceptance #it{p}_{T} vs reco centrality",nptbins,ptmin,ptmax, 100, 0, 100) ;
393  fhPrimEtaPtCentrality ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
394  fhPrimEtaAccPtCentrality->SetXTitle("#it{p}_{T} (GeV/#it{c})");
395  fhPrimEtaPtCentrality ->SetYTitle("Centrality");
396  fhPrimEtaAccPtCentrality->SetYTitle("Centrality");
397  outputContainer->Add(fhPrimEtaPtCentrality) ;
398  outputContainer->Add(fhPrimEtaAccPtCentrality) ;
399 
400  fhPrimPi0PtEventPlane = new TH2F("hPrimPi0PtEventPlane","Primary #pi^{0} #it{p}_{T} vs reco event plane angle, |#it{Y}|<1",
401  nptbins,ptmin,ptmax, 100, 0, TMath::Pi()) ;
402  fhPrimPi0AccPtEventPlane = new TH2F("hPrimPi0AccPtEventPlane","Primary #pi^{0} with both photons in acceptance #it{p}_{T} vs reco event plane angle",
403  nptbins,ptmin,ptmax, 100, 0, TMath::Pi()) ;
404  fhPrimPi0PtEventPlane ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
405  fhPrimPi0AccPtEventPlane->SetXTitle("#it{p}_{T} (GeV/#it{c})");
406  fhPrimPi0PtEventPlane ->SetYTitle("Event Plane Angle (rad)");
407  fhPrimPi0AccPtEventPlane->SetYTitle("Event Plane Angle (rad)");
408  outputContainer->Add(fhPrimPi0PtEventPlane) ;
409  outputContainer->Add(fhPrimPi0AccPtEventPlane) ;
410 
411  fhPrimEtaPtEventPlane = new TH2F("hPrimEtaPtEventPlane","Primary #eta #it{p}_{T} vs reco event plane angle, |#it{Y}|<1",nptbins,ptmin,ptmax, 100, 0, TMath::Pi()) ;
412  fhPrimEtaAccPtEventPlane = new TH2F("hPrimEtaAccPtEventPlane","Primary #eta with both #gamma_{decay} in acceptance #it{p}_{T} vs reco event plane angle",nptbins,ptmin,ptmax, 100, 0, TMath::Pi()) ;
413  fhPrimEtaPtEventPlane ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
414  fhPrimEtaAccPtEventPlane->SetXTitle("#it{p}_{T} (GeV/#it{c})");
415  fhPrimEtaPtEventPlane ->SetYTitle("Event Plane Angle (rad)");
416  fhPrimEtaAccPtEventPlane->SetYTitle("Event Plane Angle (rad)");
417  outputContainer->Add(fhPrimEtaPtEventPlane) ;
418  outputContainer->Add(fhPrimEtaAccPtEventPlane) ;
419  }
420 
421  if(fFillAngleHisto)
422  {
423  fhPrimPi0OpeningAngle = new TH2F
424  ("hPrimPi0OpeningAngle","Angle between all primary #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,100,0,0.5);
425  fhPrimPi0OpeningAngle->SetYTitle("#theta(rad)");
426  fhPrimPi0OpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
427  outputContainer->Add(fhPrimPi0OpeningAngle) ;
428 
429  fhPrimPi0OpeningAngleAsym = new TH2F
430  ("hPrimPi0OpeningAngleAsym","Angle between all primary #gamma pair vs #it{Asymmetry}, #it{p}_{T}>5 GeV/#it{c}",100,0,1,100,0,0.5);
431  fhPrimPi0OpeningAngleAsym->SetXTitle("|A|=| (E_{1}-E_{2}) / (E_{1}+E_{2}) |");
432  fhPrimPi0OpeningAngleAsym->SetYTitle("#theta(rad)");
433  outputContainer->Add(fhPrimPi0OpeningAngleAsym) ;
434 
435  fhPrimPi0CosOpeningAngle = new TH2F
436  ("hPrimPi0CosOpeningAngle","Cosinus of angle between all primary #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,100,-1,1);
437  fhPrimPi0CosOpeningAngle->SetYTitle("cos (#theta) ");
438  fhPrimPi0CosOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
439  outputContainer->Add(fhPrimPi0CosOpeningAngle) ;
440 
441  fhPrimEtaOpeningAngle = new TH2F
442  ("hPrimEtaOpeningAngle","Angle between all primary #gamma pair vs E_{#eta}",nptbins,ptmin,ptmax,100,0,0.5);
443  fhPrimEtaOpeningAngle->SetYTitle("#theta(rad)");
444  fhPrimEtaOpeningAngle->SetXTitle("E_{#eta} (GeV)");
445  outputContainer->Add(fhPrimEtaOpeningAngle) ;
446 
447  fhPrimEtaOpeningAngleAsym = new TH2F
448  ("hPrimEtaOpeningAngleAsym","Angle between all primary #gamma pair vs #it{Asymmetry}, #it{p}_{T}>5 GeV/#it{c}",100,0,1,100,0,0.5);
449  fhPrimEtaOpeningAngleAsym->SetXTitle("|#it{A}|=| (#it{E}_{1}-#it{E}_{2}) / (#it{E}_{1}+#it{E}_{2}) |");
450  fhPrimEtaOpeningAngleAsym->SetYTitle("#theta(rad)");
451  outputContainer->Add(fhPrimEtaOpeningAngleAsym) ;
452 
453 
454  fhPrimEtaCosOpeningAngle = new TH2F
455  ("hPrimEtaCosOpeningAngle","Cosinus of angle between all primary #gamma pair vs E_{#eta}",nptbins,ptmin,ptmax,100,-1,1);
456  fhPrimEtaCosOpeningAngle->SetYTitle("cos (#theta) ");
457  fhPrimEtaCosOpeningAngle->SetXTitle("#it{E}_{ #eta} (GeV)");
458  outputContainer->Add(fhPrimEtaCosOpeningAngle) ;
459  }
460 
461  // Primary origin
462  if(fFillOriginHisto)
463  {
464  // Pi0
465  fhPrimPi0PtOrigin = new TH2F("hPrimPi0PtOrigin","Primary #pi^{0} #it{p}_{T} vs origin",nptbins,ptmin,ptmax,11,0,11) ;
466  fhPrimPi0PtOrigin->SetXTitle("#it{p}_{T} (GeV/#it{c})");
467  fhPrimPi0PtOrigin->SetYTitle("Origin");
468  fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(1 ,"Status 21");
469  fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(2 ,"Quark");
470  fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(3 ,"qq Resonances ");
471  fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(4 ,"Resonances");
472  fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(5 ,"#rho");
473  fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(6 ,"#omega");
474  fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(7 ,"K");
475  fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(8 ,"Other");
476  fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(9 ,"#eta");
477  fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(10 ,"#eta prime");
478  outputContainer->Add(fhPrimPi0PtOrigin) ;
479 
480  fhPrimNotResonancePi0PtOrigin = new TH2F("hPrimNotResonancePi0PtOrigin","Primary #pi^{0} #it{p}_{T} vs origin",nptbins,ptmin,ptmax,11,0,11) ;
481  fhPrimNotResonancePi0PtOrigin->SetXTitle("#it{p}_{T} (GeV/#it{c})");
482  fhPrimNotResonancePi0PtOrigin->SetYTitle("Origin");
483  fhPrimNotResonancePi0PtOrigin->GetYaxis()->SetBinLabel(1 ,"Status 21");
484  fhPrimNotResonancePi0PtOrigin->GetYaxis()->SetBinLabel(2 ,"Quark");
485  fhPrimNotResonancePi0PtOrigin->GetYaxis()->SetBinLabel(3 ,"qq Resonances");
486  fhPrimNotResonancePi0PtOrigin->GetYaxis()->SetBinLabel(4 ,"Resonances");
487  fhPrimNotResonancePi0PtOrigin->GetYaxis()->SetBinLabel(5 ,"#rho");
488  fhPrimNotResonancePi0PtOrigin->GetYaxis()->SetBinLabel(6 ,"#omega");
489  fhPrimNotResonancePi0PtOrigin->GetYaxis()->SetBinLabel(7 ,"K");
490  fhPrimNotResonancePi0PtOrigin->GetYaxis()->SetBinLabel(8 ,"Other");
491  fhPrimNotResonancePi0PtOrigin->GetYaxis()->SetBinLabel(9 ,"#eta");
492  fhPrimNotResonancePi0PtOrigin->GetYaxis()->SetBinLabel(10 ,"#eta prime");
493  outputContainer->Add(fhPrimNotResonancePi0PtOrigin) ;
494 
495  fhPrimPi0PtStatus = new TH2F("hPrimPi0PtStatus","Primary #pi^{0} #it{p}_{T} vs status",nptbins,ptmin,ptmax,101,-50,50) ;
496  fhPrimPi0PtStatus->SetXTitle("#it{p}_{T} (GeV/#it{c})");
497  fhPrimPi0PtStatus->SetYTitle("Status");
498  outputContainer->Add(fhPrimPi0PtStatus) ;
499 
500  // Eta
501  fhPrimEtaPtOrigin = new TH2F("hPrimEtaPtOrigin","Primary #pi^{0} #it{p}_{T} vs origin",nptbins,ptmin,ptmax,7,0,7) ;
502  fhPrimEtaPtOrigin->SetXTitle("#it{p}_{T} (GeV/#it{c})");
503  fhPrimEtaPtOrigin->SetYTitle("Origin");
504  fhPrimEtaPtOrigin->GetYaxis()->SetBinLabel(1 ,"Status 21");
505  fhPrimEtaPtOrigin->GetYaxis()->SetBinLabel(2 ,"Quark");
506  fhPrimEtaPtOrigin->GetYaxis()->SetBinLabel(3 ,"qq Resonances");
507  fhPrimEtaPtOrigin->GetYaxis()->SetBinLabel(4 ,"Resonances");
508  fhPrimEtaPtOrigin->GetYaxis()->SetBinLabel(5 ,"Other");
509  fhPrimEtaPtOrigin->GetYaxis()->SetBinLabel(6 ,"#eta prime ");
510  outputContainer->Add(fhPrimEtaPtOrigin) ;
511 
512  // Production vertex
513  fhPrimPi0ProdVertex = new TH2F("hPrimPi0ProdVertex","generated #pi^{0} #it{p}_{T} vs production vertex",
514  200,0.,20.,5000,0,500) ;
515  fhPrimPi0ProdVertex->SetXTitle("#it{p}_{T} (GeV/#it{c})");
516  fhPrimPi0ProdVertex->SetYTitle("#it{R} (cm)");
517  outputContainer->Add(fhPrimPi0ProdVertex) ;
518 
519  fhPrimEtaProdVertex = new TH2F("hPrimEtaProdVertex","generated #eta #it{p}_{T} vs production vertex",
520  200,0.,20.,5000,0,500) ;
521  fhPrimEtaProdVertex->SetXTitle("#it{p}_{T} (GeV/#it{c})");
522  fhPrimEtaProdVertex->SetYTitle("#it{R} (cm)");
523  outputContainer->Add(fhPrimEtaProdVertex) ;
524  }
525 
526  if(fFillArmenterosThetaStar)
527  {
528  TString ebin[] = {"8 < E < 12 GeV","12 < E < 16 GeV", "16 < E < 20 GeV", "E > 20 GeV" };
529  Int_t narmbins = 400;
530  Float_t armmin = 0;
531  Float_t armmax = 0.4;
532 
533  for(Int_t i = 0; i < 4; i++)
534  {
535  fhArmPrimPi0[i] = new TH2F(Form("hArmenterosPrimPi0EBin%d",i),
536  Form("Armenteros of primary #pi^{0}, %s",ebin[i].Data()),
537  200, -1, 1, narmbins,armmin,armmax);
538  fhArmPrimPi0[i]->SetYTitle("#it{p}_{T}^{Arm}");
539  fhArmPrimPi0[i]->SetXTitle("#alpha^{Arm}");
540  outputContainer->Add(fhArmPrimPi0[i]) ;
541 
542  fhArmPrimEta[i] = new TH2F(Form("hArmenterosPrimEtaEBin%d",i),
543  Form("Armenteros of primary #eta, %s",ebin[i].Data()),
544  200, -1, 1, narmbins,armmin,armmax);
545  fhArmPrimEta[i]->SetYTitle("#it{p}_{T}^{Arm}");
546  fhArmPrimEta[i]->SetXTitle("#alpha^{Arm}");
547  outputContainer->Add(fhArmPrimEta[i]) ;
548  }
549 
550  // Same as asymmetry ...
551  fhCosThStarPrimPi0 = new TH2F
552  ("hCosThStarPrimPi0","cos(#theta *) for primary #pi^{0}",nptbins,ptmin,ptmax,200,-1,1);
553  fhCosThStarPrimPi0->SetYTitle("cos(#theta *)");
554  fhCosThStarPrimPi0->SetXTitle("E_{ #pi^{0}} (GeV)");
555  outputContainer->Add(fhCosThStarPrimPi0) ;
556 
557  fhCosThStarPrimEta = new TH2F
558  ("hCosThStarPrimEta","cos(#theta *) for primary #eta",nptbins,ptmin,ptmax,200,-1,1);
559  fhCosThStarPrimEta->SetYTitle("cos(#theta *)");
560  fhCosThStarPrimEta->SetXTitle("E_{ #eta} (GeV)");
561  outputContainer->Add(fhCosThStarPrimEta) ;
562  }
563 
564  if(fFillOnlyMCAcceptanceHisto) return outputContainer;
565  }
566 
567  //
568  // Create mixed event containers
569  //
570  fEventsList = new TList*[GetNCentrBin()*GetNZvertBin()*GetNRPBin()] ;
571 
572  for(Int_t ic=0; ic<GetNCentrBin(); ic++)
573  {
574  for(Int_t iz=0; iz<GetNZvertBin(); iz++)
575  {
576  for(Int_t irp=0; irp<GetNRPBin(); irp++)
577  {
578  Int_t bin = GetEventMixBin(ic,iz,irp);
579  fEventsList[bin] = new TList() ;
580  fEventsList[bin]->SetOwner(kFALSE);
581  }
582  }
583  }
584 
585  // Init the number of modules, set in the class AliCalorimeterUtils
586  fNModules = GetCaloUtils()->GetNumberOfSuperModulesUsed();
587  if(GetCalorimeter()==kPHOS && fNModules > 4) fNModules = 4;
588 
589  fhReMod = new TH2F*[fNModules] ;
590  fhMiMod = new TH2F*[fNModules] ;
591 
592  if(GetCalorimeter() == kPHOS)
593  {
594  fhReDiffPHOSMod = new TH2F*[fNModules] ;
595  fhMiDiffPHOSMod = new TH2F*[fNModules] ;
596  }
597  else
598  {
599  fhReSameSectorEMCALMod = new TH2F*[fNModules/2] ;
600  fhReSameSideEMCALMod = new TH2F*[fNModules-2] ;
601  fhMiSameSectorEMCALMod = new TH2F*[fNModules/2] ;
602  fhMiSameSideEMCALMod = new TH2F*[fNModules-2] ;
603  }
604 
605 
606  fhRe1 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
607  fhMi1 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
608  if(fFillBadDistHisto)
609  {
610  fhRe2 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
611  fhRe3 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
612  fhMi2 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
613  fhMi3 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
614  }
615  if(fMakeInvPtPlots)
616  {
617  fhReInvPt1 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
618  fhMiInvPt1 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
619  if(fFillBadDistHisto){
620  fhReInvPt2 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
621  fhReInvPt3 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
622  fhMiInvPt2 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
623  fhMiInvPt3 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
624  }
625  }
626 
627  if(fCheckConversion)
628  {
629  fhReConv = new TH2F("hReConv","Real Pair with one recombined conversion ",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
630  fhReConv->SetXTitle("#it{p}_{T} (GeV/#it{c})");
631  fhReConv->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
632  outputContainer->Add(fhReConv) ;
633 
634  fhReConv2 = new TH2F("hReConv2","Real Pair with 2 recombined conversion ",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
635  fhReConv2->SetXTitle("#it{p}_{T} (GeV/#it{c})");
636  fhReConv2->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
637  outputContainer->Add(fhReConv2) ;
638 
639  if(DoOwnMix())
640  {
641  fhMiConv = new TH2F("hMiConv","Mixed Pair with one recombined conversion ",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
642  fhMiConv->SetXTitle("#it{p}_{T} (GeV/#it{c})");
643  fhMiConv->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
644  outputContainer->Add(fhMiConv) ;
645 
646  fhMiConv2 = new TH2F("hMiConv2","Mixed Pair with 2 recombined conversion ",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
647  fhMiConv2->SetXTitle("#it{p}_{T} (GeV/#it{c})");
648  fhMiConv2->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
649  outputContainer->Add(fhMiConv2) ;
650  }
651  }
652 
653  for(Int_t ic=0; ic<GetNCentrBin(); ic++)
654  {
655  for(Int_t ipid=0; ipid<fNPIDBits; ipid++)
656  {
657  for(Int_t iasym=0; iasym<fNAsymCuts; iasym++)
658  {
659  Int_t index = ((ic*fNPIDBits)+ipid)*fNAsymCuts + iasym;
660  //printf("cen %d, pid %d, asy %d, Index %d\n",ic,ipid,iasym,index);
661  //Distance to bad module 1
662  snprintf(key, buffersize,"hRe_cen%d_pidbit%d_asy%d_dist1",ic,ipid,iasym) ;
663  snprintf(title, buffersize,"Real #it{M}_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 1",
664  ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
665  fhRe1[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
666  fhRe1[index]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
667  fhRe1[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
668  //printf("name: %s\n ",fhRe1[index]->GetName());
669  outputContainer->Add(fhRe1[index]) ;
670 
671  if(fFillBadDistHisto)
672  {
673  //Distance to bad module 2
674  snprintf(key, buffersize,"hRe_cen%d_pidbit%d_asy%d_dist2",ic,ipid,iasym) ;
675  snprintf(title, buffersize,"Real #it{M}_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 2",
676  ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
677  fhRe2[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
678  fhRe2[index]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
679  fhRe2[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
680  outputContainer->Add(fhRe2[index]) ;
681 
682  //Distance to bad module 3
683  snprintf(key, buffersize,"hRe_cen%d_pidbit%d_asy%d_dist3",ic,ipid,iasym) ;
684  snprintf(title, buffersize,"Real #it{M}_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 3",
685  ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
686  fhRe3[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
687  fhRe3[index]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
688  fhRe3[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
689  outputContainer->Add(fhRe3[index]) ;
690  }
691 
692  //Inverse pT
693  if(fMakeInvPtPlots)
694  {
695  //Distance to bad module 1
696  snprintf(key, buffersize,"hReInvPt_cen%d_pidbit%d_asy%d_dist1",ic,ipid,iasym) ;
697  snprintf(title, buffersize,"Real #it{M}_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 1",
698  ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
699  fhReInvPt1[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
700  fhReInvPt1[index]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
701  fhReInvPt1[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
702  outputContainer->Add(fhReInvPt1[index]) ;
703 
704  if(fFillBadDistHisto){
705  //Distance to bad module 2
706  snprintf(key, buffersize,"hReInvPt_cen%d_pidbit%d_asy%d_dist2",ic,ipid,iasym) ;
707  snprintf(title, buffersize,"Real #it{M}_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 2",
708  ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
709  fhReInvPt2[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
710  fhReInvPt2[index]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
711  fhReInvPt2[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
712  outputContainer->Add(fhReInvPt2[index]) ;
713 
714  //Distance to bad module 3
715  snprintf(key, buffersize,"hReInvPt_cen%d_pidbit%d_asy%d_dist3",ic,ipid,iasym) ;
716  snprintf(title, buffersize,"Real #it{M}_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 3",
717  ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
718  fhReInvPt3[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
719  fhReInvPt3[index]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
720  fhReInvPt3[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
721  outputContainer->Add(fhReInvPt3[index]) ;
722  }
723  }
724 
725  if(DoOwnMix())
726  {
727  //Distance to bad module 1
728  snprintf(key, buffersize,"hMi_cen%d_pidbit%d_asy%d_dist1",ic,ipid,iasym) ;
729  snprintf(title, buffersize,"Mixed #it{M}_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 1",
730  ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
731  fhMi1[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
732  fhMi1[index]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
733  fhMi1[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
734  outputContainer->Add(fhMi1[index]) ;
735  if(fFillBadDistHisto){
736  //Distance to bad module 2
737  snprintf(key, buffersize,"hMi_cen%d_pidbit%d_asy%d_dist2",ic,ipid,iasym) ;
738  snprintf(title, buffersize,"Mixed #it{M}_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 2",
739  ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
740  fhMi2[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
741  fhMi2[index]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
742  fhMi2[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
743  outputContainer->Add(fhMi2[index]) ;
744 
745  //Distance to bad module 3
746  snprintf(key, buffersize,"hMi_cen%d_pidbit%d_asy%d_dist3",ic,ipid,iasym) ;
747  snprintf(title, buffersize,"Mixed #it{M}_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 3",
748  ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
749  fhMi3[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
750  fhMi3[index]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
751  fhMi3[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
752  outputContainer->Add(fhMi3[index]) ;
753  }
754 
755  // Inverse pT
756  if(fMakeInvPtPlots)
757  {
758  // Distance to bad module 1
759  snprintf(key, buffersize,"hMiInvPt_cen%d_pidbit%d_asy%d_dist1",ic,ipid,iasym) ;
760  snprintf(title, buffersize,"Mixed #it{M}_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 1",
761  ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
762  fhMiInvPt1[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
763  fhMiInvPt1[index]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
764  fhMiInvPt1[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
765  outputContainer->Add(fhMiInvPt1[index]) ;
766  if(fFillBadDistHisto){
767  // Distance to bad module 2
768  snprintf(key, buffersize,"hMiInvPt_cen%d_pidbit%d_asy%d_dist2",ic,ipid,iasym) ;
769  snprintf(title, buffersize,"Mixed #it{M}_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 2",
770  ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
771  fhMiInvPt2[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
772  fhMiInvPt2[index]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
773  fhMiInvPt2[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
774  outputContainer->Add(fhMiInvPt2[index]) ;
775 
776  // Distance to bad module 3
777  snprintf(key, buffersize,"hMiInvPt_cen%d_pidbit%d_asy%d_dist3",ic,ipid,iasym) ;
778  snprintf(title, buffersize,"Mixed #it{M}_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f,dist bad 3",
779  ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
780  fhMiInvPt3[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
781  fhMiInvPt3[index]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
782  fhMiInvPt3[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
783  outputContainer->Add(fhMiInvPt3[index]) ;
784  }
785  }
786  }
787  }
788  }
789  }
790 
791  fhEPairDiffTime = new TH2F("hEPairDiffTime","cluster pair time difference vs #it{p}_{T}",nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
792  fhEPairDiffTime->SetXTitle("#it{p}_{T,pair} (GeV/#it{c})");
793  fhEPairDiffTime->SetYTitle("#Delta t (ns)");
794  outputContainer->Add(fhEPairDiffTime);
795 
796  if(fFillAsymmetryHisto)
797  {
798  fhRePtAsym = new TH2F("hRePtAsym","#it{Asymmetry} vs #it{p}_{T} , for pairs",nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax) ;
799  fhRePtAsym->SetXTitle("#it{p}_{T} (GeV/#it{c})");
800  fhRePtAsym->SetYTitle("#it{Asymmetry}");
801  outputContainer->Add(fhRePtAsym);
802 
803  fhRePtAsymPi0 = new TH2F("hRePtAsymPi0","#it{Asymmetry} vs #it{p}_{T} , for pairs close to #pi^{0} mass",nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax) ;
804  fhRePtAsymPi0->SetXTitle("#it{p}_{T} (GeV/#it{c})");
805  fhRePtAsymPi0->SetYTitle("Asymmetry");
806  outputContainer->Add(fhRePtAsymPi0);
807 
808  fhRePtAsymEta = new TH2F("hRePtAsymEta","#it{Asymmetry} vs #it{p}_{T} , for pairs close to #eta mass",nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax) ;
809  fhRePtAsymEta->SetXTitle("#it{p}_{T} (GeV/#it{c})");
810  fhRePtAsymEta->SetYTitle("#it{Asymmetry}");
811  outputContainer->Add(fhRePtAsymEta);
812  }
813 
814  if(fMultiCutAna)
815  {
816  fhRePIDBits = new TH2F*[fNPIDBits];
817  for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
818  snprintf(key, buffersize,"hRe_pidbit%d",ipid) ;
819  snprintf(title, buffersize,"Real #it{M}_{#gamma#gamma} distr. for PIDBit=%d",fPIDBits[ipid]) ;
820  fhRePIDBits[ipid] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
821  fhRePIDBits[ipid]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
822  fhRePIDBits[ipid]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
823  outputContainer->Add(fhRePIDBits[ipid]) ;
824  }// pid bit loop
825 
826  fhRePtNCellAsymCuts = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
827  fhMiPtNCellAsymCuts = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
828 
829  if(fFillSMCombinations)
830  for(Int_t iSM = 0; iSM < fNModules; iSM++) fhRePtNCellAsymCutsSM[iSM] = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
831 
832 
833  for(Int_t ipt=0; ipt<fNPtCuts; ipt++)
834  {
835  for(Int_t icell=0; icell<fNCellNCuts; icell++)
836  {
837  for(Int_t iasym=0; iasym<fNAsymCuts; iasym++)
838  {
839  snprintf(key, buffersize,"hRe_pt%d_cell%d_asym%d",ipt,icell,iasym) ;
840  snprintf(title, buffersize,"Real #it{M}_{#gamma#gamma} distr. for pt >%2.2f, ncell>%d and asym >%1.2f ",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]) ;
841  Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
842  //printf("ipt %d, icell %d, iassym %d, index %d\n",ipt, icell, iasym, index);
843  fhRePtNCellAsymCuts[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
844  fhRePtNCellAsymCuts[index]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
845  fhRePtNCellAsymCuts[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
846  outputContainer->Add(fhRePtNCellAsymCuts[index]) ;
847 
848  snprintf(key, buffersize,"hMi_pt%d_cell%d_asym%d",ipt,icell,iasym) ;
849  snprintf(title, buffersize,"Mixed #it{M}_{#gamma#gamma} distr. for pt >%2.2f, ncell>%d and asym >%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]) ;
850  fhMiPtNCellAsymCuts[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
851  fhMiPtNCellAsymCuts[index]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
852  fhMiPtNCellAsymCuts[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
853  outputContainer->Add(fhMiPtNCellAsymCuts[index]) ;
854 
855  if(fFillSMCombinations)
856  {
857  for(Int_t iSM = 0; iSM < fNModules; iSM++)
858  {
859  snprintf(key, buffersize,"hRe_pt%d_cell%d_asym%d_SM%d",ipt,icell,iasym,iSM) ;
860  snprintf(title, buffersize,"Real #it{M}_{#gamma#gamma} distr. for pt >%2.2f, ncell>%d and asym >%1.2f, SM %d ",
861  fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym],iSM) ;
862  fhRePtNCellAsymCutsSM[iSM][index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
863  fhRePtNCellAsymCutsSM[iSM][index]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
864  fhRePtNCellAsymCutsSM[iSM][index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
865  outputContainer->Add(fhRePtNCellAsymCutsSM[iSM][index]) ;
866  }
867  }
868  }
869  }
870  }
871 
872  if(ntrmbins!=0)
873  {
874  fhRePtMult = new TH3F*[fNAsymCuts] ;
875  for(Int_t iasym = 0; iasym<fNAsymCuts; iasym++)
876  {
877  fhRePtMult[iasym] = new TH3F(Form("hRePtMult_asym%d",iasym),Form("(#it{p}_{T},C,M)_{#gamma#gamma}, A<%1.2f",fAsymCuts[iasym]),
878  nptbins,ptmin,ptmax,ntrmbins,ntrmmin,ntrmmax,nmassbins,massmin,massmax);
879  fhRePtMult[iasym]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
880  fhRePtMult[iasym]->SetYTitle("Track multiplicity");
881  fhRePtMult[iasym]->SetZTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
882  outputContainer->Add(fhRePtMult[iasym]) ;
883  }
884  }
885  } // multi cuts analysis
886 
887  if(fFillSSCombinations)
888  {
889  fhReSS[0] = new TH2F("hRe_SS_Tight"," 0.01 < #lambda_{0}^{2} < 0.4",
890  nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
891  fhReSS[0]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
892  fhReSS[0]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
893  outputContainer->Add(fhReSS[0]) ;
894 
895 
896  fhReSS[1] = new TH2F("hRe_SS_Loose"," #lambda_{0}^{2} > 0.4",
897  nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
898  fhReSS[1]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
899  fhReSS[1]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
900  outputContainer->Add(fhReSS[1]) ;
901 
902 
903  fhReSS[2] = new TH2F("hRe_SS_Both"," cluster_{1} #lambda_{0}^{2} > 0.4; cluster_{2} 0.01 < #lambda_{0}^{2} < 0.4",
904  nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
905  fhReSS[2]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
906  fhReSS[2]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
907  outputContainer->Add(fhReSS[2]) ;
908  }
909 
910  if(DoOwnMix())
911  {
912  fhEventBin=new TH1I("hEventBin","Number of real pairs per bin(cen,vz,rp)",
913  GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1,0,
914  GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1) ;
915  fhEventBin->SetXTitle("bin");
916  outputContainer->Add(fhEventBin) ;
917 
918 
919  fhEventMixBin=new TH1I("hEventMixBin","Number of mixed pairs per bin(cen,vz,rp)",
920  GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1,0,
921  GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1) ;
922  fhEventMixBin->SetXTitle("bin");
923  outputContainer->Add(fhEventMixBin) ;
924  }
925 
926  if( IsHighMultiplicityAnalysisOn() )
927  {
928  fhCentrality=new TH1F("hCentralityBin","Number of events in centrality bin",GetNCentrBin(),0.,1.*GetNCentrBin()) ;
929  fhCentrality->SetXTitle("Centrality bin");
930  outputContainer->Add(fhCentrality) ;
931 
932  fhCentralityNoPair=new TH1F("hCentralityBinNoPair","Number of events in centrality bin, with no cluster pairs",GetNCentrBin(),0.,1.*GetNCentrBin()) ;
933  fhCentralityNoPair->SetXTitle("Centrality bin");
934  outputContainer->Add(fhCentralityNoPair) ;
935 
936  fhEventPlaneResolution=new TH2F("hEventPlaneResolution","Event plane resolution",GetNCentrBin(),0,GetNCentrBin(),100,0.,TMath::TwoPi()) ;
937  fhEventPlaneResolution->SetYTitle("Resolution");
938  fhEventPlaneResolution->SetXTitle("Centrality Bin");
939  outputContainer->Add(fhEventPlaneResolution) ;
940  }
941 
942  if(fFillAngleHisto)
943  {
944  fhRealOpeningAngle = new TH2F
945  ("hRealOpeningAngle","Angle between all #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,300,0,TMath::Pi());
946  fhRealOpeningAngle->SetYTitle("#theta(rad)");
947  fhRealOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
948  outputContainer->Add(fhRealOpeningAngle) ;
949 
950  fhRealCosOpeningAngle = new TH2F
951  ("hRealCosOpeningAngle","Cosinus of angle between all #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,100,0,1);
952  fhRealCosOpeningAngle->SetYTitle("cos (#theta) ");
953  fhRealCosOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
954  outputContainer->Add(fhRealCosOpeningAngle) ;
955 
956  if(DoOwnMix())
957  {
958  fhMixedOpeningAngle = new TH2F
959  ("hMixedOpeningAngle","Angle between all #gamma pair vs E_{#pi^{0}}, Mixed pairs",nptbins,ptmin,ptmax,300,0,TMath::Pi());
960  fhMixedOpeningAngle->SetYTitle("#theta(rad)");
961  fhMixedOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
962  outputContainer->Add(fhMixedOpeningAngle) ;
963 
964  fhMixedCosOpeningAngle = new TH2F
965  ("hMixedCosOpeningAngle","Cosinus of angle between all #gamma pair vs E_{#pi^{0}}, Mixed pairs",nptbins,ptmin,ptmax,100,0,1);
966  fhMixedCosOpeningAngle->SetYTitle("cos (#theta) ");
967  fhMixedCosOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
968  outputContainer->Add(fhMixedCosOpeningAngle) ;
969  }
970  }
971 
972  // Histograms filled only if MC data is requested
973  if(IsDataMC())
974  {
975  fhReMCFromConversion = new TH2F("hReMCFromConversion","Invariant mass of 2 clusters originated in conversions",
976  nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
977  fhReMCFromConversion->SetXTitle("#it{p}_{T} (GeV/#it{c})");
978  fhReMCFromConversion->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
979  outputContainer->Add(fhReMCFromConversion) ;
980 
981  fhReMCFromNotConversion = new TH2F("hReMCNotFromConversion","Invariant mass of 2 clusters not originated in conversions",
982  nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
983  fhReMCFromNotConversion->SetXTitle("#it{p}_{T} (GeV/#it{c})");
984  fhReMCFromNotConversion->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
985  outputContainer->Add(fhReMCFromNotConversion) ;
986 
987  fhReMCFromMixConversion = new TH2F("hReMCFromMixConversion","Invariant mass of 2 clusters one from conversion and the other not",
988  nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
989  fhReMCFromMixConversion->SetXTitle("#it{p}_{T} (GeV/#it{c})");
990  fhReMCFromMixConversion->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
991  outputContainer->Add(fhReMCFromMixConversion) ;
992 
993  if(fFillOriginHisto)
994  {
995 
996  fhMCPi0PtOrigin = new TH2F("hMCPi0PtOrigin","Reconstructed pair from generated #pi^{0} #it{p}_{T} vs origin",nptbins,ptmin,ptmax,11,0,11) ;
997  fhMCPi0PtOrigin->SetXTitle("#it{p}_{T} (GeV/#it{c})");
998  fhMCPi0PtOrigin->SetYTitle("Origin");
999  fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(1 ,"Status 21");
1000  fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(2 ,"Quark");
1001  fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(3 ,"qq Resonances");
1002  fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(4 ,"Resonances");
1003  fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(5 ,"#rho");
1004  fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(6 ,"#omega");
1005  fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(7 ,"K");
1006  fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(8 ,"Other");
1007  fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(9 ,"#eta");
1008  fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(10 ,"#eta prime");
1009  outputContainer->Add(fhMCPi0PtOrigin) ;
1010 
1011  fhMCNotResonancePi0PtOrigin = new TH2F("hMCNotResonancePi0PtOrigin","Reconstructed pair from generated #pi^{0} #it{p}_{T} vs origin",nptbins,ptmin,ptmax,11,0,11) ;
1012  fhMCNotResonancePi0PtOrigin->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1013  fhMCNotResonancePi0PtOrigin->SetYTitle("Origin");
1014  fhMCNotResonancePi0PtOrigin->GetYaxis()->SetBinLabel(1 ,"Status 21");
1015  fhMCNotResonancePi0PtOrigin->GetYaxis()->SetBinLabel(2 ,"Quark");
1016  fhMCNotResonancePi0PtOrigin->GetYaxis()->SetBinLabel(3 ,"qq Resonances");
1017  fhMCNotResonancePi0PtOrigin->GetYaxis()->SetBinLabel(4 ,"Resonances");
1018  fhMCNotResonancePi0PtOrigin->GetYaxis()->SetBinLabel(5 ,"#rho");
1019  fhMCNotResonancePi0PtOrigin->GetYaxis()->SetBinLabel(6 ,"#omega");
1020  fhMCNotResonancePi0PtOrigin->GetYaxis()->SetBinLabel(7 ,"K");
1021  fhMCNotResonancePi0PtOrigin->GetYaxis()->SetBinLabel(8 ,"Other");
1022  fhMCNotResonancePi0PtOrigin->GetYaxis()->SetBinLabel(9 ,"#eta");
1023  fhMCNotResonancePi0PtOrigin->GetYaxis()->SetBinLabel(10 ,"#eta prime");
1024  outputContainer->Add(fhMCNotResonancePi0PtOrigin) ;
1025 
1026  fhMCPi0PtStatus = new TH2F("hMCPi0PtStatus","Reconstructed pair from generated #pi^{0} #it{p}_{T} vs status",nptbins,ptmin,ptmax,101,-50,50) ;
1027  fhMCPi0PtStatus->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1028  fhMCPi0PtStatus->SetYTitle("Status");
1029  outputContainer->Add(fhMCPi0PtStatus) ;
1030 
1031  // Eta
1032 
1033  fhMCEtaPtOrigin = new TH2F("hMCEtaPtOrigin","Reconstructed pair from generated #pi^{0} #it{p}_{T} vs origin",nptbins,ptmin,ptmax,7,0,7) ;
1034  fhMCEtaPtOrigin->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1035  fhMCEtaPtOrigin->SetYTitle("Origin");
1036  fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(1 ,"Status 21");
1037  fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(2 ,"Quark");
1038  fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(3 ,"qq Resonances");
1039  fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(4 ,"Resonances");
1040  fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(5 ,"Other");
1041  fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(6 ,"#eta prime");
1042  outputContainer->Add(fhMCEtaPtOrigin) ;
1043 
1044  fhMCPi0ProdVertex = new TH2F("hMCPi0ProdVertex","Selected reco pair from generated #pi^{0} #it{p}_{T} vs production vertex",
1045  200,0.,20.,5000,0,500) ;
1046  fhMCPi0ProdVertex->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1047  fhMCPi0ProdVertex->SetYTitle("#it{R} (cm)");
1048  outputContainer->Add(fhMCPi0ProdVertex) ;
1049 
1050  fhMCEtaProdVertex = new TH2F("hMCEtaProdVertex","Selected reco pair from generated #eta #it{p}_{T} vs production vertex",
1051  200,0.,20.,5000,0,500) ;
1052  fhMCEtaProdVertex->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1053  fhMCEtaProdVertex->SetYTitle("#it{R} (cm)");
1054  outputContainer->Add(fhMCEtaProdVertex) ;
1055 
1056  for(Int_t i = 0; i<13; i++)
1057  {
1058  fhMCOrgMass[i] = new TH2F(Form("hMCOrgMass_%d",i),Form("#it{M} vs #it{p}_{T}, origin %d",i),nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1059  fhMCOrgMass[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1060  fhMCOrgMass[i]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1061  outputContainer->Add(fhMCOrgMass[i]) ;
1062 
1063  fhMCOrgAsym[i]= new TH2F(Form("hMCOrgAsym_%d",i),Form("#it{Asymmetry} vs #it{p}_{T}, origin %d",i),nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax) ;
1064  fhMCOrgAsym[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1065  fhMCOrgAsym[i]->SetYTitle("A");
1066  outputContainer->Add(fhMCOrgAsym[i]) ;
1067 
1068  fhMCOrgDeltaEta[i] = new TH2F(Form("hMCOrgDeltaEta_%d",i),Form("#Delta #eta of pair vs #it{p}_{T}, origin %d",i),nptbins,ptmin,ptmax,netabins,-1.4,1.4) ;
1069  fhMCOrgDeltaEta[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1070  fhMCOrgDeltaEta[i]->SetYTitle("#Delta #eta");
1071  outputContainer->Add(fhMCOrgDeltaEta[i]) ;
1072 
1073  fhMCOrgDeltaPhi[i]= new TH2F(Form("hMCOrgDeltaPhi_%d",i),Form("#Delta #phi of pair vs #it{p}_{T}, origin %d",i),nptbins,ptmin,ptmax,nphibins,-0.7,0.7) ;
1074  fhMCOrgDeltaPhi[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1075  fhMCOrgDeltaPhi[i]->SetYTitle("#Delta #phi (rad)");
1076  outputContainer->Add(fhMCOrgDeltaPhi[i]) ;
1077  }
1078 
1079  if(fMultiCutAnaSim)
1080  {
1081  fhMCPi0MassPtTrue = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
1082  fhMCPi0MassPtRec = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
1083  fhMCPi0PtTruePtRec = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
1084  fhMCEtaMassPtRec = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
1085  fhMCEtaMassPtTrue = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
1086  fhMCEtaPtTruePtRec = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
1087 
1088  for(Int_t ipt=0; ipt<fNPtCuts; ipt++)
1089  {
1090  for(Int_t icell=0; icell<fNCellNCuts; icell++)
1091  {
1092  for(Int_t iasym=0; iasym<fNAsymCuts; iasym++)
1093  {
1094  Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
1095 
1096  fhMCPi0MassPtRec[index] = new TH2F(Form("hMCPi0MassPtRec_pt%d_cell%d_asym%d",ipt,icell,iasym),
1097  Form("Reconstructed #it{M} vs reconstructed #it{p}_{T} of true #pi^{0} cluster pairs for #it{p}_{T} >%2.2f, #it{N}^{cluster}_{cell}>%d and |#it{A}|>%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]),
1098  nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1099  fhMCPi0MassPtRec[index]->SetXTitle("#it{p}_{T, reconstructed} (GeV/#it{c})");
1100  fhMCPi0MassPtRec[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1101  outputContainer->Add(fhMCPi0MassPtRec[index]) ;
1102 
1103  fhMCPi0MassPtTrue[index] = new TH2F(Form("hMCPi0MassPtTrue_pt%d_cell%d_asym%d",ipt,icell,iasym),
1104  Form("Reconstructed #it{M} vs generated #it{p}_{T} of true #pi^{0} cluster pairs for #it{p}_{T} >%2.2f, #it{N}^{cluster}_{cell}>%d and |#it{A}|>%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]),
1105  nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1106  fhMCPi0MassPtTrue[index]->SetXTitle("#it{p}_{T, generated} (GeV/#it{c})");
1107  fhMCPi0MassPtTrue[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1108  outputContainer->Add(fhMCPi0MassPtTrue[index]) ;
1109 
1110  fhMCPi0PtTruePtRec[index] = new TH2F(Form("hMCPi0PtTruePtRec_pt%d_cell%d_asym%d",ipt,icell,iasym),
1111  Form("Generated vs reconstructed #it{p}_{T} of true #pi^{0} cluster pairs, 0.01 < rec. mass < 0.17 MeV/#it{c}^{2} for #it{p}_{T} >%2.2f, #it{N}^{cluster}_{cell}>%d and |#it{A}|>%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]),
1112  nptbins,ptmin,ptmax,nptbins,ptmin,ptmax) ;
1113  fhMCPi0PtTruePtRec[index]->SetXTitle("#it{p}_{T, generated} (GeV/#it{c})");
1114  fhMCPi0PtTruePtRec[index]->SetYTitle("#it{p}_{T, reconstructed} (GeV/#it{c})");
1115  outputContainer->Add(fhMCPi0PtTruePtRec[index]) ;
1116 
1117  fhMCEtaMassPtRec[index] = new TH2F(Form("hMCEtaMassPtRec_pt%d_cell%d_asym%d",ipt,icell,iasym),
1118  Form("Reconstructed #it{M} vs reconstructed #it{p}_{T} of true #eta cluster pairs for #it{p}_{T} >%2.2f, #it{N}^{cluster}_{cell}>%d and |#it{A}|>%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]),
1119  nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1120  fhMCEtaMassPtRec[index]->SetXTitle("#it{p}_{T, generated} (GeV/#it{c})");
1121  fhMCEtaMassPtRec[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1122  outputContainer->Add(fhMCEtaMassPtRec[index]) ;
1123 
1124  fhMCEtaMassPtTrue[index] = new TH2F(Form("hMCEtaMassPtTrue_pt%d_cell%d_asym%d",ipt,icell,iasym),
1125  Form("Reconstructed #it{M} vs generated #it{p}_{T} of true #eta cluster pairs for #it{p}_{T} >%2.2f, #it{N}^{cluster}_{cell}>%d and |#it{A}|>%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]),
1126  nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1127  fhMCEtaMassPtTrue[index]->SetXTitle("#it{p}_{T, generated} (GeV/#it{c})");
1128  fhMCEtaMassPtTrue[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1129  outputContainer->Add(fhMCEtaMassPtTrue[index]) ;
1130 
1131  fhMCEtaPtTruePtRec[index] = new TH2F(Form("hMCEtaPtTruePtRec_pt%d_cell%d_asym%d",ipt,icell,iasym),
1132  Form("Generated vs reconstructed #it{p}_{T} of true #eta cluster pairs, 0.01 < rec. mass < 0.17 MeV/#it{c}^{2} for #it{p}_{T} >%2.2f, ncell>%d and asym >%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]),
1133  nptbins,ptmin,ptmax,nptbins,ptmin,ptmax) ;
1134  fhMCEtaPtTruePtRec[index]->SetXTitle("#it{p}_{T, generated} (GeV/#it{c})");
1135  fhMCEtaPtTruePtRec[index]->SetYTitle("#it{p}_{T, reconstructed} (GeV/#it{c})");
1136  outputContainer->Add(fhMCEtaPtTruePtRec[index]) ;
1137  }
1138  }
1139  }
1140  }//multi cut ana
1141  else
1142  {
1143  fhMCPi0MassPtTrue = new TH2F*[1];
1144  fhMCPi0PtTruePtRec = new TH2F*[1];
1145  fhMCEtaMassPtTrue = new TH2F*[1];
1146  fhMCEtaPtTruePtRec = new TH2F*[1];
1147 
1148  fhMCPi0MassPtTrue[0] = new TH2F("hMCPi0MassPtTrue","Reconstructed Mass vs generated p_T of true #pi^{0} cluster pairs",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1149  fhMCPi0MassPtTrue[0]->SetXTitle("#it{p}_{T, generated} (GeV/#it{c})");
1150  fhMCPi0MassPtTrue[0]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1151  outputContainer->Add(fhMCPi0MassPtTrue[0]) ;
1152 
1153  fhMCPi0PtTruePtRec[0]= new TH2F("hMCPi0PtTruePtRec","Generated vs reconstructed p_T of true #pi^{0} cluster pairs, 0.01 < rec. mass < 0.17 MeV/#it{c}^{2}",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax) ;
1154  fhMCPi0PtTruePtRec[0]->SetXTitle("#it{p}_{T, generated} (GeV/#it{c})");
1155  fhMCPi0PtTruePtRec[0]->SetYTitle("#it{p}_{T, reconstructed} (GeV/#it{c})");
1156  outputContainer->Add(fhMCPi0PtTruePtRec[0]) ;
1157 
1158  fhMCEtaMassPtTrue[0] = new TH2F("hMCEtaMassPtTrue","Reconstructed Mass vs generated p_T of true #eta cluster pairs",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1159  fhMCEtaMassPtTrue[0]->SetXTitle("#it{p}_{T, generated} (GeV/#it{c})");
1160  fhMCEtaMassPtTrue[0]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1161  outputContainer->Add(fhMCEtaMassPtTrue[0]) ;
1162 
1163  fhMCEtaPtTruePtRec[0]= new TH2F("hMCEtaPtTruePtRec","Generated vs reconstructed p_T of true #eta cluster pairs, 0.01 < rec. mass < 0.17 MeV/#it{c}^{2}",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax) ;
1164  fhMCEtaPtTruePtRec[0]->SetXTitle("#it{p}_{T, generated} (GeV/#it{c})");
1165  fhMCEtaPtTruePtRec[0]->SetYTitle("#it{p}_{T, reconstructed} (GeV/#it{c})");
1166  outputContainer->Add(fhMCEtaPtTruePtRec[0]) ;
1167  }
1168  }
1169  }
1170 
1171  if(fFillSMCombinations)
1172  {
1173  TString pairnamePHOS[] = {"(0-1)","(0-2)","(1-2)","(0-3)","(0-4)","(1-3)","(1-4)","(2-3)","(2-4)","(3-4)"};
1174  for(Int_t imod=0; imod<fNModules; imod++)
1175  {
1176  //Module dependent invariant mass
1177  snprintf(key, buffersize,"hReMod_%d",imod) ;
1178  snprintf(title, buffersize,"Real #it{M}_{#gamma#gamma} distr. for Module %d",imod) ;
1179  fhReMod[imod] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1180  fhReMod[imod]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1181  fhReMod[imod]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1182  outputContainer->Add(fhReMod[imod]) ;
1183  if(GetCalorimeter()==kPHOS)
1184  {
1185  snprintf(key, buffersize,"hReDiffPHOSMod_%d",imod) ;
1186  snprintf(title, buffersize,"Real pairs PHOS, clusters in different Modules: %s",(pairnamePHOS[imod]).Data()) ;
1187  fhReDiffPHOSMod[imod] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1188  fhReDiffPHOSMod[imod]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1189  fhReDiffPHOSMod[imod]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1190  outputContainer->Add(fhReDiffPHOSMod[imod]) ;
1191  }
1192  else
1193  {//EMCAL
1194  if(imod<fNModules/2)
1195  {
1196  snprintf(key, buffersize,"hReSameSectorEMCAL_%d",imod) ;
1197  snprintf(title, buffersize,"Real pairs EMCAL, clusters in same sector, SM(%d,%d)",imod*2,imod*2+1) ;
1198  fhReSameSectorEMCALMod[imod] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1199  fhReSameSectorEMCALMod[imod]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1200  fhReSameSectorEMCALMod[imod]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1201  outputContainer->Add(fhReSameSectorEMCALMod[imod]) ;
1202  }
1203  if(imod<fNModules-2)
1204  {
1205  snprintf(key, buffersize,"hReSameSideEMCAL_%d",imod) ;
1206  snprintf(title, buffersize,"Real pairs EMCAL, clusters in same side SM(%d,%d)",imod, imod+2) ;
1207  fhReSameSideEMCALMod[imod] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1208  fhReSameSideEMCALMod[imod]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1209  fhReSameSideEMCALMod[imod]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1210  outputContainer->Add(fhReSameSideEMCALMod[imod]) ;
1211  }
1212  }//EMCAL
1213 
1214  if(DoOwnMix())
1215  {
1216  snprintf(key, buffersize,"hMiMod_%d",imod) ;
1217  snprintf(title, buffersize,"Mixed #it{M}_{#gamma#gamma} distr. for Module %d",imod) ;
1218  fhMiMod[imod] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1219  fhMiMod[imod]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1220  fhMiMod[imod]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1221  outputContainer->Add(fhMiMod[imod]) ;
1222 
1223  if(GetCalorimeter()==kPHOS)
1224  {
1225  snprintf(key, buffersize,"hMiDiffPHOSMod_%d",imod) ;
1226  snprintf(title, buffersize,"Mixed pairs PHOS, clusters in different Modules: %s",(pairnamePHOS[imod]).Data()) ;
1227  fhMiDiffPHOSMod[imod] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1228  fhMiDiffPHOSMod[imod]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1229  fhMiDiffPHOSMod[imod]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1230  outputContainer->Add(fhMiDiffPHOSMod[imod]) ;
1231  }//PHOS
1232  else
1233  {//EMCAL
1234  if(imod<fNModules/2)
1235  {
1236  snprintf(key, buffersize,"hMiSameSectorEMCALMod_%d",imod) ;
1237  snprintf(title, buffersize,"Mixed pairs EMCAL, clusters in same sector, SM(%d,%d)",imod*2,imod*2+1) ;
1238  fhMiSameSectorEMCALMod[imod] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1239  fhMiSameSectorEMCALMod[imod]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1240  fhMiSameSectorEMCALMod[imod]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1241  outputContainer->Add(fhMiSameSectorEMCALMod[imod]) ;
1242  }
1243  if(imod<fNModules-2){
1244 
1245  snprintf(key, buffersize,"hMiSameSideEMCALMod_%d",imod) ;
1246  snprintf(title, buffersize,"Mixed pairs EMCAL, clusters in same side SM(%d,%d)",imod, imod+2) ;
1247  fhMiSameSideEMCALMod[imod] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1248  fhMiSameSideEMCALMod[imod]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1249  fhMiSameSideEMCALMod[imod]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1250  outputContainer->Add(fhMiSameSideEMCALMod[imod]) ;
1251  }
1252  } // EMCAL
1253  } // own mix
1254  } // loop combinations
1255  } // SM combinations
1256 
1257 // for(Int_t i = 0; i < outputContainer->GetEntries() ; i++){
1258 //
1259 // printf("Histogram %d, name: %s\n ",i, outputContainer->At(i)->GetName());
1260 //
1261 // }
1262 
1263  return outputContainer;
1264 }
1265 
1266 //___________________________________________________
1268 //___________________________________________________
1269 void AliAnaPi0::Print(const Option_t * /*opt*/) const
1270 {
1271  printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
1272  AliAnaCaloTrackCorrBaseClass::Print(" ");
1273 
1274  printf("Number of bins in Centrality: %d \n",GetNCentrBin()) ;
1275  printf("Number of bins in Z vert. pos: %d \n",GetNZvertBin()) ;
1276  printf("Number of bins in Reac. Plain: %d \n",GetNRPBin()) ;
1277  printf("Depth of event buffer: %d \n",GetNMaxEvMix()) ;
1278  printf("Pair in same Module: %d \n",fSameSM) ;
1279  printf("Cuts: \n") ;
1280  // printf("Z vertex position: -%2.3f < z < %2.3f \n",GetZvertexCut(),GetZvertexCut()) ; //It crashes here, why?
1281  printf("Number of modules: %d \n",fNModules) ;
1282  printf("Select pairs with their angle: %d, edep %d, min angle %2.3f, max angle %2.3f \n",fUseAngleCut, fUseAngleEDepCut, fAngleCut, fAngleMaxCut) ;
1283  printf("Asymmetry cuts: n = %d, \n",fNAsymCuts) ;
1284  printf("\tasymmetry < ");
1285  for(Int_t i = 0; i < fNAsymCuts; i++) printf("%2.2f ",fAsymCuts[i]);
1286  printf("\n");
1287 
1288  printf("PID selection bits: n = %d, \n",fNPIDBits) ;
1289  printf("\tPID bit = ");
1290  for(Int_t i = 0; i < fNPIDBits; i++) printf("%d ",fPIDBits[i]);
1291  printf("\n");
1292 
1293  if(fMultiCutAna)
1294  {
1295  printf("pT cuts: n = %d, \n",fNPtCuts) ;
1296  printf("\tpT > ");
1297  for(Int_t i = 0; i < fNPtCuts; i++) printf("%2.2f ",fPtCuts[i]);
1298  printf("GeV/c\n");
1299 
1300  printf("N cell in cluster cuts: n = %d, \n",fNCellNCuts) ;
1301  printf("\tnCell > ");
1302  for(Int_t i = 0; i < fNCellNCuts; i++) printf("%d ",fCellNCuts[i]);
1303  printf("\n");
1304 
1305  }
1306  printf("------------------------------------------------------\n") ;
1307 }
1308 
1309 //________________________________________
1311 //________________________________________
1312 void AliAnaPi0::FillAcceptanceHistograms()
1313 {
1314  Double_t mesonY = -100 ;
1315  Double_t mesonE = -1 ;
1316  Double_t mesonPt = -1 ;
1317  Double_t mesonPhi = 100 ;
1318  Double_t mesonYeta= -1 ;
1319 
1320  Int_t pdg = 0 ;
1321  Int_t nprim = 0 ;
1322  Int_t nDaught = 0 ;
1323  Int_t iphot1 = 0 ;
1324  Int_t iphot2 = 0 ;
1325 
1326  Float_t cen = GetEventCentrality();
1327  Float_t ep = GetEventPlaneAngle();
1328 
1329  TParticle * primStack = 0;
1330  AliAODMCParticle * primAOD = 0;
1331 
1332  // Get the ESD MC particles container
1333  AliStack * stack = 0;
1334  if( GetReader()->ReadStack() )
1335  {
1336  stack = GetMCStack();
1337  if(!stack ) return;
1338  nprim = stack->GetNtrack();
1339  }
1340 
1341  // Get the AOD MC particles container
1342  TClonesArray * mcparticles = 0;
1343  if( GetReader()->ReadAODMCParticles() )
1344  {
1345  mcparticles = GetReader()->GetAODMCParticles();
1346  if( !mcparticles ) return;
1347  nprim = mcparticles->GetEntriesFast();
1348  }
1349 
1350  for(Int_t i=0 ; i < nprim; i++)
1351  {
1352  if(GetReader()->AcceptOnlyHIJINGLabels() && !GetReader()->IsHIJINGLabel(i)) continue ;
1353 
1354  if(GetReader()->ReadStack())
1355  {
1356  primStack = stack->Particle(i) ;
1357  if(!primStack)
1358  {
1359  AliWarning("ESD primaries pointer not available!!");
1360  continue;
1361  }
1362 
1363  // If too small skip
1364  if( primStack->Energy() < 0.4 ) continue;
1365 
1366  pdg = primStack->GetPdgCode();
1367  nDaught = primStack->GetNDaughters();
1368  iphot1 = primStack->GetDaughter(0) ;
1369  iphot2 = primStack->GetDaughter(1) ;
1370  if(primStack->Energy() == TMath::Abs(primStack->Pz())) continue ; //Protection against floating point exception
1371 
1372  //printf("i %d, %s %d %s %d \n",i, stack->Particle(i)->GetName(), stack->Particle(i)->GetPdgCode(),
1373  // prim->GetName(), prim->GetPdgCode());
1374 
1375  //Photon kinematics
1376  primStack->Momentum(fPi0Mom);
1377 
1378  mesonY = 0.5*TMath::Log((primStack->Energy()+primStack->Pz())/(primStack->Energy()-primStack->Pz())) ;
1379  }
1380  else
1381  {
1382  primAOD = (AliAODMCParticle *) mcparticles->At(i);
1383  if(!primAOD)
1384  {
1385  AliWarning("AOD primaries pointer not available!!");
1386  continue;
1387  }
1388 
1389  // If too small skip
1390  if( primAOD->E() < 0.4 ) continue;
1391 
1392  pdg = primAOD->GetPdgCode();
1393  nDaught = primAOD->GetNDaughters();
1394  iphot1 = primAOD->GetFirstDaughter() ;
1395  iphot2 = primAOD->GetLastDaughter() ;
1396 
1397  if(primAOD->E() == TMath::Abs(primAOD->Pz())) continue ; //Protection against floating point exception
1398 
1399  //Photon kinematics
1400  fPi0Mom.SetPxPyPzE(primAOD->Px(),primAOD->Py(),primAOD->Pz(),primAOD->E());
1401 
1402  mesonY = 0.5*TMath::Log((primAOD->E()+primAOD->Pz())/(primAOD->E()-primAOD->Pz())) ;
1403  }
1404 
1405  // Select only pi0 or eta
1406  if( pdg != 111 && pdg != 221) continue ;
1407 
1408  mesonPt = fPi0Mom.Pt () ;
1409  mesonE = fPi0Mom.E () ;
1410  mesonYeta= fPi0Mom.Eta() ;
1411  mesonPhi = fPi0Mom.Phi() ;
1412  if( mesonPhi < 0 ) mesonPhi+=TMath::TwoPi();
1413  mesonPhi *= TMath::RadToDeg();
1414 
1415  if(pdg == 111)
1416  {
1417  if(TMath::Abs(mesonY) < 1.0)
1418  {
1419  fhPrimPi0E ->Fill(mesonE , GetEventWeight()) ;
1420  fhPrimPi0Pt ->Fill(mesonPt, GetEventWeight()) ;
1421  fhPrimPi0Phi->Fill(mesonPt, mesonPhi, GetEventWeight()) ;
1422 
1423  fhPrimPi0YetaYcut->Fill(mesonPt, mesonYeta, GetEventWeight()) ;
1424 
1425  if( IsHighMultiplicityAnalysisOn() )
1426  {
1427  fhPrimPi0PtCentrality->Fill(mesonPt, cen, GetEventWeight()) ;
1428  fhPrimPi0PtEventPlane->Fill(mesonPt, ep , GetEventWeight()) ;
1429  }
1430  }
1431 
1432  fhPrimPi0Y ->Fill(mesonPt, mesonY , GetEventWeight()) ;
1433  fhPrimPi0Yeta->Fill(mesonPt, mesonYeta, GetEventWeight()) ;
1434  }
1435  else if(pdg == 221)
1436  {
1437  if(TMath::Abs(mesonY) < 1.0)
1438  {
1439  fhPrimEtaE ->Fill(mesonE , GetEventWeight()) ;
1440  fhPrimEtaPt ->Fill(mesonPt, GetEventWeight()) ;
1441  fhPrimEtaPhi->Fill(mesonPt, mesonPhi, GetEventWeight()) ;
1442 
1443  fhPrimEtaYetaYcut->Fill(mesonPt, mesonYeta, GetEventWeight()) ;
1444 
1445  if( IsHighMultiplicityAnalysisOn() )
1446  {
1447  fhPrimEtaPtCentrality->Fill(mesonPt, cen, GetEventWeight()) ;
1448  fhPrimEtaPtEventPlane->Fill(mesonPt, ep , GetEventWeight()) ;
1449  }
1450  }
1451 
1452  fhPrimEtaY ->Fill(mesonPt, mesonY , GetEventWeight()) ;
1453  fhPrimEtaYeta->Fill(mesonPt, mesonYeta, GetEventWeight()) ;
1454  }
1455 
1456  //Origin of meson
1457  if(fFillOriginHisto && TMath::Abs(mesonY) < 0.7)
1458  {
1459  Int_t momindex = -1;
1460  Int_t mompdg = -1;
1461  Int_t momstatus = -1;
1462  Int_t status = -1;
1463  Float_t momR = 0;
1464  if(GetReader()->ReadStack()) momindex = primStack->GetFirstMother();
1465  if(GetReader()->ReadAODMCParticles()) momindex = primAOD ->GetMother();
1466 
1467  if(momindex >= 0 && momindex < nprim)
1468  {
1469  if(GetReader()->ReadStack())
1470  {
1471  status = primStack->GetStatusCode();
1472  TParticle* mother = stack->Particle(momindex);
1473  mompdg = TMath::Abs(mother->GetPdgCode());
1474  momstatus = mother->GetStatusCode();
1475  momR = mother->R();
1476  }
1477 
1478  if(GetReader()->ReadAODMCParticles())
1479  {
1480  status = primAOD->GetStatus();
1481  AliAODMCParticle* mother = (AliAODMCParticle*) mcparticles->At(momindex);
1482  mompdg = TMath::Abs(mother->GetPdgCode());
1483  momstatus = mother->GetStatus();
1484  momR = TMath::Sqrt(mother->Xv()*mother->Xv()+mother->Yv()*mother->Yv());
1485  }
1486 
1487  if(pdg == 111)
1488  {
1489  fhPrimPi0ProdVertex->Fill(mesonPt, momR , GetEventWeight());
1490  fhPrimPi0PtStatus ->Fill(mesonPt, status, GetEventWeight());
1491 
1492 
1493  if (momstatus == 21) fhPrimPi0PtOrigin->Fill(mesonPt, 0.5, GetEventWeight());//parton
1494  else if(mompdg < 22 ) fhPrimPi0PtOrigin->Fill(mesonPt, 1.5, GetEventWeight());//quark
1495  else if(mompdg > 2100 && mompdg < 2210)
1496  fhPrimPi0PtOrigin->Fill(mesonPt, 2.5, GetEventWeight());// resonances
1497  else if(mompdg == 221) fhPrimPi0PtOrigin->Fill(mesonPt, 8.5, GetEventWeight());//eta
1498  else if(mompdg == 331) fhPrimPi0PtOrigin->Fill(mesonPt, 9.5, GetEventWeight());//eta prime
1499  else if(mompdg == 213) fhPrimPi0PtOrigin->Fill(mesonPt, 4.5, GetEventWeight());//rho
1500  else if(mompdg == 223) fhPrimPi0PtOrigin->Fill(mesonPt, 5.5, GetEventWeight());//omega
1501  else if(mompdg >= 310 && mompdg <= 323)
1502  fhPrimPi0PtOrigin->Fill(mesonPt, 6.5, GetEventWeight());//k0S, k+-,k*
1503  else if(mompdg == 130) fhPrimPi0PtOrigin->Fill(mesonPt, 6.5, GetEventWeight());//k0L
1504  else if(momstatus == 11 || momstatus == 12 )
1505  fhPrimPi0PtOrigin->Fill(mesonPt, 3.5, GetEventWeight());//resonances
1506  else fhPrimPi0PtOrigin->Fill(mesonPt, 7.5, GetEventWeight());//other?
1507 
1508  //printf("Prim Pi0: index %d, pt %2.2f Prod vertex %3.3f, origin pdg %d, origin status %d, origin UI %d\n",
1509  // momindex, mesonPt,mother->R(),mompdg,momstatus,mother->GetUniqueID());
1510 
1511  if(status!=11)
1512  {
1513  if (momstatus == 21) fhPrimNotResonancePi0PtOrigin->Fill(mesonPt, 0.5, GetEventWeight());//parton
1514  else if(mompdg < 22 ) fhPrimNotResonancePi0PtOrigin->Fill(mesonPt, 1.5, GetEventWeight());//quark
1515  else if(mompdg > 2100 && mompdg < 2210)
1516  fhPrimNotResonancePi0PtOrigin->Fill(mesonPt, 2.5, GetEventWeight());// resonances
1517  else if(mompdg == 221) fhPrimNotResonancePi0PtOrigin->Fill(mesonPt, 8.5, GetEventWeight());//eta
1518  else if(mompdg == 331) fhPrimNotResonancePi0PtOrigin->Fill(mesonPt, 9.5, GetEventWeight());//eta prime
1519  else if(mompdg == 213) fhPrimNotResonancePi0PtOrigin->Fill(mesonPt, 4.5, GetEventWeight());//rho
1520  else if(mompdg == 223) fhPrimNotResonancePi0PtOrigin->Fill(mesonPt, 5.5, GetEventWeight());//omega
1521  else if(mompdg >= 310 && mompdg <= 323)
1522  fhPrimNotResonancePi0PtOrigin->Fill(mesonPt, 6.5, GetEventWeight());//k0S, k+-,k*
1523  else if(mompdg == 130) fhPrimNotResonancePi0PtOrigin->Fill(mesonPt, 6.5, GetEventWeight());//k0L
1524  else if(momstatus == 11 || momstatus == 12 )
1525  fhPrimNotResonancePi0PtOrigin->Fill(mesonPt, 3.5, GetEventWeight());//resonances
1526  else fhPrimNotResonancePi0PtOrigin->Fill(mesonPt, 7.5, GetEventWeight());//other?
1527  }
1528 
1529  }//pi0
1530  else
1531  {
1532  if (momstatus == 21 ) fhPrimEtaPtOrigin->Fill(mesonPt, 0.5, GetEventWeight());//parton
1533  else if(mompdg < 22 ) fhPrimEtaPtOrigin->Fill(mesonPt, 1.5, GetEventWeight());//quark
1534  else if(mompdg > 2100 && mompdg < 2210)
1535  fhPrimEtaPtOrigin->Fill(mesonPt, 2.5, GetEventWeight());//qq resonances
1536  else if(mompdg == 331) fhPrimEtaPtOrigin->Fill(mesonPt, 5.5, GetEventWeight());//eta prime
1537  else if(momstatus == 11 || momstatus == 12 )
1538  fhPrimEtaPtOrigin->Fill(mesonPt, 3.5, GetEventWeight());//resonances
1539  else fhPrimEtaPtOrigin->Fill(mesonPt, 4.5, GetEventWeight());//stable, conversions?
1540  //printf("Other Meson pdg %d, Mother %s, pdg %d, status %d\n",pdg, TDatabasePDG::Instance()->GetParticle(mompdg)->GetName(),mompdg, momstatus );
1541 
1542  fhPrimEtaProdVertex->Fill(mesonPt, momR, GetEventWeight());
1543  }
1544  } // pi0 has mother
1545  }
1546 
1547  //Check if both photons hit Calorimeter
1548  if(nDaught != 2 ) continue; //Only interested in 2 gamma decay
1549 
1550  if(iphot1 < 0 || iphot1 >= nprim || iphot2 < 0 || iphot2 >= nprim) continue ;
1551 
1552  Int_t pdg1 = 0;
1553  Int_t pdg2 = 0;
1554  Bool_t inacceptance1 = kTRUE;
1555  Bool_t inacceptance2 = kTRUE;
1556 
1557  if(GetReader()->ReadStack())
1558  {
1559  TParticle * phot1 = stack->Particle(iphot1) ;
1560  TParticle * phot2 = stack->Particle(iphot2) ;
1561 
1562  if(!phot1 || !phot2) continue ;
1563 
1564  pdg1 = phot1->GetPdgCode();
1565  pdg2 = phot2->GetPdgCode();
1566 
1567  phot1->Momentum(fPhotonMom1);
1568  phot2->Momentum(fPhotonMom2);
1569 
1570  // Check if photons hit the Calorimeter acceptance
1571  if(IsRealCaloAcceptanceOn())
1572  {
1573  if( !GetCaloUtils()->IsMCParticleInCalorimeterAcceptance( GetCalorimeter(), phot1 )) inacceptance1 = kFALSE ;
1574  if( !GetCaloUtils()->IsMCParticleInCalorimeterAcceptance( GetCalorimeter(), phot2 )) inacceptance2 = kFALSE ;
1575  }
1576  }
1577 
1578  if(GetReader()->ReadAODMCParticles())
1579  {
1580  AliAODMCParticle * phot1 = (AliAODMCParticle *) mcparticles->At(iphot1) ;
1581  AliAODMCParticle * phot2 = (AliAODMCParticle *) mcparticles->At(iphot2) ;
1582 
1583  if(!phot1 || !phot2) continue ;
1584 
1585  pdg1 = phot1->GetPdgCode();
1586  pdg2 = phot2->GetPdgCode();
1587 
1588  fPhotonMom1.SetPxPyPzE(phot1->Px(),phot1->Py(),phot1->Pz(),phot1->E());
1589  fPhotonMom2.SetPxPyPzE(phot2->Px(),phot2->Py(),phot2->Pz(),phot2->E());
1590 
1591  // Check if photons hit the Calorimeter acceptance
1592  if(IsRealCaloAcceptanceOn())
1593  {
1594  if( !GetCaloUtils()->IsMCParticleInCalorimeterAcceptance( GetCalorimeter(), phot1 )) inacceptance1 = kFALSE ;
1595  if( !GetCaloUtils()->IsMCParticleInCalorimeterAcceptance( GetCalorimeter(), phot2 )) inacceptance2 = kFALSE ;
1596  }
1597  }
1598 
1599  if( pdg1 != 22 || pdg2 !=22) continue ;
1600 
1601  // Check if photons hit desired acceptance in the fidutial borders fixed in the analysis
1602  if(IsFiducialCutOn())
1603  {
1604  if( inacceptance1 && !GetFiducialCut()->IsInFiducialCut(fPhotonMom1.Eta(), fPhotonMom1.Phi(), GetCalorimeter()) ) inacceptance1 = kFALSE ;
1605  if( inacceptance2 && !GetFiducialCut()->IsInFiducialCut(fPhotonMom2.Eta(), fPhotonMom2.Phi(), GetCalorimeter()) ) inacceptance2 = kFALSE ;
1606  }
1607 
1608  if(fFillArmenterosThetaStar) FillArmenterosThetaStar(pdg);
1609 
1610  if(GetCalorimeter()==kEMCAL && inacceptance1 && inacceptance2 && fCheckAccInSector)
1611  {
1612  Int_t absID1=0;
1613  Int_t absID2=0;
1614 
1615  Float_t photonPhi1 = fPhotonMom1.Phi();
1616  Float_t photonPhi2 = fPhotonMom2.Phi();
1617 
1618  if(photonPhi1 < 0) photonPhi1+=TMath::TwoPi();
1619  if(photonPhi2 < 0) photonPhi2+=TMath::TwoPi();
1620 
1621  GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(fPhotonMom1.Eta(),photonPhi1,absID1);
1622  GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(fPhotonMom2.Eta(),photonPhi2,absID2);
1623 
1624  Int_t sm1 = GetEMCALGeometry()->GetSuperModuleNumber(absID1);
1625  Int_t sm2 = GetEMCALGeometry()->GetSuperModuleNumber(absID2);
1626 
1627  Int_t j=0;
1628  Bool_t sameSector = kFALSE;
1629  for(Int_t isector = 0; isector < fNModules/2; isector++)
1630  {
1631  j=2*isector;
1632  if((sm1==j && sm2==j+1) || (sm1==j+1 && sm2==j)) sameSector = kTRUE;
1633  }
1634 
1635  if(sm1!=sm2 && !sameSector)
1636  {
1637  inacceptance1 = kFALSE;
1638  inacceptance2 = kFALSE;
1639  }
1640  //if(sm1!=sm2)printf("sm1 %d, sm2 %d, same sector %d, in acceptance %d\n",sm1,sm2,sameSector,inacceptance);
1641  // if(GetEMCALGeometry()->Impact(phot1) && GetEMCALGeometry()->Impact(phot2))
1642  // inacceptance = kTRUE;
1643  }
1644 
1645  AliDebug(2,Form("Accepted in %s?: m (%2.2f,%2.2f,%2.2f), p1 (%2.2f,%2.2f,%2.2f), p2 (%2.2f,%2.2f,%2.2f) : in1 %d, in2 %d",
1646  GetCalorimeterString().Data(),
1647  mesonPt,mesonYeta,mesonPhi,
1648  fPhotonMom1.Pt(),fPhotonMom1.Eta(),fPhotonMom1.Phi()*TMath::RadToDeg(),
1649  fPhotonMom2.Pt(),fPhotonMom2.Eta(),fPhotonMom2.Phi()*TMath::RadToDeg(),
1650  inacceptance1, inacceptance2));
1651 
1652  if(inacceptance1 && inacceptance2)
1653  {
1654  Float_t asym = TMath::Abs((fPhotonMom1.E()-fPhotonMom2.E()) / (fPhotonMom1.E()+fPhotonMom2.E()));
1655  Double_t angle = fPhotonMom1.Angle(fPhotonMom2.Vect());
1656 
1657  AliDebug(2,Form("\t ACCEPTED pdg %d: pt %2.2f, phi %2.2f, eta %2.2f",pdg,mesonPt,mesonPhi,mesonYeta));
1658 
1659  if(pdg==111)
1660  {
1661  fhPrimPi0AccE ->Fill(mesonE , GetEventWeight()) ;
1662  fhPrimPi0AccPt ->Fill(mesonPt, GetEventWeight()) ;
1663  fhPrimPi0AccPhi ->Fill(mesonPt, mesonPhi , GetEventWeight()) ;
1664  fhPrimPi0AccY ->Fill(mesonPt, mesonY , GetEventWeight()) ;
1665  fhPrimPi0AccYeta->Fill(mesonPt, mesonYeta, GetEventWeight()) ;
1666 
1667  if( IsHighMultiplicityAnalysisOn() )
1668  {
1669  fhPrimPi0AccPtCentrality->Fill(mesonPt, cen, GetEventWeight()) ;
1670  fhPrimPi0AccPtEventPlane->Fill(mesonPt, ep , GetEventWeight()) ;
1671  }
1672 
1673  if(fFillAngleHisto)
1674  {
1675  fhPrimPi0OpeningAngle ->Fill(mesonPt, angle, GetEventWeight());
1676  if(mesonPt > 5)fhPrimPi0OpeningAngleAsym->Fill(asym, angle, GetEventWeight());
1677  fhPrimPi0CosOpeningAngle ->Fill(mesonPt, TMath::Cos(angle), GetEventWeight());
1678  }
1679  }
1680  else if(pdg==221)
1681  {
1682  fhPrimEtaAccE ->Fill(mesonE , GetEventWeight()) ;
1683  fhPrimEtaAccPt ->Fill(mesonPt, GetEventWeight()) ;
1684  fhPrimEtaAccPhi ->Fill(mesonPt, mesonPhi , GetEventWeight()) ;
1685  fhPrimEtaAccY ->Fill(mesonPt, mesonY , GetEventWeight()) ;
1686  fhPrimEtaAccYeta->Fill(mesonPt, mesonYeta, GetEventWeight()) ;
1687 
1688  if( IsHighMultiplicityAnalysisOn() )
1689  {
1690  fhPrimEtaAccPtCentrality->Fill(mesonPt, cen, GetEventWeight()) ;
1691  fhPrimEtaAccPtEventPlane->Fill(mesonPt, ep , GetEventWeight()) ;
1692  }
1693 
1694  if(fFillAngleHisto)
1695  {
1696  fhPrimEtaOpeningAngle ->Fill(mesonPt, angle, GetEventWeight());
1697  if(mesonPt > 5)fhPrimEtaOpeningAngleAsym->Fill(asym, angle, GetEventWeight());
1698  fhPrimEtaCosOpeningAngle ->Fill(mesonPt, TMath::Cos(angle), GetEventWeight());
1699  }
1700  }
1701  } // Accepted
1702  } // loop on primaries
1703 }
1704 
1705 //________________________________________________
1707 //________________________________________________
1708 void AliAnaPi0::FillArmenterosThetaStar(Int_t pdg)
1709 {
1710  // Get pTArm and AlphaArm
1711  Float_t momentumSquaredMother = fPi0Mom.P()*fPi0Mom.P();
1712  Float_t momentumDaughter1AlongMother = 0.;
1713  Float_t momentumDaughter2AlongMother = 0.;
1714 
1715  if (momentumSquaredMother > 0.)
1716  {
1717  momentumDaughter1AlongMother = (fPhotonMom1.Px()*fPi0Mom.Px() + fPhotonMom1.Py()*fPi0Mom.Py()+ fPhotonMom1.Pz()*fPi0Mom.Pz()) / sqrt(momentumSquaredMother);
1718  momentumDaughter2AlongMother = (fPhotonMom2.Px()*fPi0Mom.Px() + fPhotonMom2.Py()*fPi0Mom.Py()+ fPhotonMom2.Pz()*fPi0Mom.Pz()) / sqrt(momentumSquaredMother);
1719  }
1720 
1721  Float_t momentumSquaredDaughter1 = fPhotonMom1.P()*fPhotonMom1.P();
1722  Float_t ptArmSquared = momentumSquaredDaughter1 - momentumDaughter1AlongMother*momentumDaughter1AlongMother;
1723 
1724  Float_t pTArm = 0.;
1725  if (ptArmSquared > 0.)
1726  pTArm = sqrt(ptArmSquared);
1727 
1728  Float_t alphaArm = 0.;
1729  if(momentumDaughter1AlongMother + momentumDaughter2AlongMother > 0)
1730  alphaArm = (momentumDaughter1AlongMother -momentumDaughter2AlongMother) / (momentumDaughter1AlongMother + momentumDaughter2AlongMother);
1731 
1732  fPhotonMom1Boost = fPhotonMom1;
1733  fPhotonMom1Boost.Boost(-fPi0Mom.BoostVector());
1734  Float_t cosThStar=TMath::Cos(fPhotonMom1Boost.Vect().Angle(fPi0Mom.Vect()));
1735 
1736  Float_t en = fPi0Mom.Energy();
1737  Int_t ebin = -1;
1738  if(en > 8 && en <= 12) ebin = 0;
1739  if(en > 12 && en <= 16) ebin = 1;
1740  if(en > 16 && en <= 20) ebin = 2;
1741  if(en > 20) ebin = 3;
1742  if(ebin < 0 || ebin > 3) return ;
1743 
1744  if(pdg==111)
1745  {
1746  fhCosThStarPrimPi0->Fill(en, cosThStar, GetEventWeight());
1747  fhArmPrimPi0[ebin]->Fill(alphaArm, pTArm , GetEventWeight());
1748  }
1749  else
1750  {
1751  fhCosThStarPrimEta->Fill(en, cosThStar, GetEventWeight());
1752  fhArmPrimEta[ebin]->Fill(alphaArm,pTArm , GetEventWeight());
1753  }
1754 
1755  // if(GetDebug() > 2 )
1756  // {
1757  // Float_t asym = 0;
1758  // if(fPhotonMom1.E()+fPhotonMom2.E() > 0) asym = TMath::Abs(fPhotonMom1.E()-fPhotonMom2.E())/(fPhotonMom1.E()+fPhotonMom2.E());
1759  //
1760  // printf("AliAnaPi0::FillArmenterosThetaStar() - E %f, alphaArm %f, pTArm %f, cos(theta*) %f, asymmetry %1.3f\n",
1761  // en,alphaArm,pTArm,cosThStar,asym);
1762  // }
1763 }
1764 
1765 //_______________________________________________________________________________________
1770 //_______________________________________________________________________________________
1771 void AliAnaPi0::FillMCVersusRecDataHistograms(Int_t index1, Int_t index2,
1772  Float_t pt1, Float_t pt2,
1773  Int_t ncell1, Int_t ncell2,
1774  Double_t mass, Double_t pt, Double_t asym,
1775  Double_t deta, Double_t dphi)
1776 {
1777  Int_t ancPDG = 0;
1778  Int_t ancStatus = 0;
1779  Int_t ancLabel = GetMCAnalysisUtils()->CheckCommonAncestor(index1, index2,
1780  GetReader(), ancPDG, ancStatus,fPi0Mom, fProdVertex);
1781 
1782  Int_t momindex = -1;
1783  Int_t mompdg = -1;
1784  Int_t momstatus = -1;
1785  Int_t status = -1;
1786  Float_t prodR = -1;
1787  Int_t mcIndex = -1;
1788 
1789  if(ancLabel > -1)
1790  {
1791  AliDebug(1,Form("Common ancestor label %d, pdg %d, name %s, status %d",
1792  ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus));
1793 
1794  if(ancPDG==22)
1795  {//gamma
1796  mcIndex = 0;
1797  }
1798  else if(TMath::Abs(ancPDG)==11)
1799  {//e
1800  mcIndex = 1;
1801  }
1802  else if(ancPDG==111)
1803  {//Pi0
1804  mcIndex = 2;
1805  if(fMultiCutAnaSim)
1806  {
1807  for(Int_t ipt=0; ipt<fNPtCuts; ipt++)
1808  {
1809  for(Int_t icell=0; icell<fNCellNCuts; icell++)
1810  {
1811  for(Int_t iasym=0; iasym<fNAsymCuts; iasym++)
1812  {
1813  Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
1814  if(pt1 > fPtCuts[ipt] && pt2 > fPtCuts[ipt] &&
1815  asym < fAsymCuts[iasym] &&
1816  ncell1 >= fCellNCuts[icell] && ncell2 >= fCellNCuts[icell])
1817  {
1818  fhMCPi0MassPtRec [index]->Fill(pt, mass, GetEventWeight());
1819  fhMCPi0MassPtTrue[index]->Fill(fPi0Mom.Pt(),mass, GetEventWeight());
1820  if(mass < 0.17 && mass > 0.1)
1821  fhMCPi0PtTruePtRec[index]->Fill(fPi0Mom.Pt(), pt, GetEventWeight());
1822  }//pass the different cuts
1823  }// pid bit cut loop
1824  }// icell loop
1825  }// pt cut loop
1826  }// Multi cut ana sim
1827  else
1828  {
1829  fhMCPi0MassPtTrue[0]->Fill(fPi0Mom.Pt(), mass, GetEventWeight());
1830 
1831  if(mass < 0.17 && mass > 0.1)
1832  {
1833  fhMCPi0PtTruePtRec[0]->Fill(fPi0Mom.Pt(), pt, GetEventWeight());
1834 
1835  //Int_t uniqueId = -1;
1836  if(GetReader()->ReadStack())
1837  {
1838  TParticle* ancestor = GetMCStack()->Particle(ancLabel);
1839  status = ancestor->GetStatusCode();
1840  momindex = ancestor->GetFirstMother();
1841  if(momindex < 0) return;
1842  TParticle* mother = GetMCStack()->Particle(momindex);
1843  mompdg = TMath::Abs(mother->GetPdgCode());
1844  momstatus = mother->GetStatusCode();
1845  prodR = mother->R();
1846  //uniqueId = mother->GetUniqueID();
1847  }
1848  else
1849  {
1850  TClonesArray * mcparticles = GetReader()->GetAODMCParticles();
1851  AliAODMCParticle* ancestor = (AliAODMCParticle *) mcparticles->At(ancLabel);
1852  status = ancestor->GetStatus();
1853  momindex = ancestor->GetMother();
1854  if(momindex < 0) return;
1855  AliAODMCParticle* mother = (AliAODMCParticle *) mcparticles->At(momindex);
1856  mompdg = TMath::Abs(mother->GetPdgCode());
1857  momstatus = mother->GetStatus();
1858  prodR = TMath::Sqrt(mother->Xv()*mother->Xv()+mother->Yv()*mother->Yv());
1859  //uniqueId = mother->GetUniqueID();
1860  }
1861 
1862  // printf("Reco Pi0: pt %2.2f Prod vertex %3.3f, origin pdg %d, origin status %d, origin UI %d\n",
1863  // pt,prodR,mompdg,momstatus,uniqueId);
1864 
1865  fhMCPi0ProdVertex->Fill(pt, prodR , GetEventWeight());
1866  fhMCPi0PtStatus ->Fill(pt, status, GetEventWeight());
1867 
1868  if (momstatus == 21) fhMCPi0PtOrigin->Fill(pt, 0.5, GetEventWeight());//parton
1869  else if(mompdg < 22 ) fhMCPi0PtOrigin->Fill(pt, 1.5, GetEventWeight());//quark
1870  else if(mompdg > 2100 && mompdg < 2210)
1871  fhMCPi0PtOrigin->Fill(pt, 2.5, GetEventWeight());// resonances
1872  else if(mompdg == 221) fhMCPi0PtOrigin->Fill(pt, 8.5, GetEventWeight());//eta
1873  else if(mompdg == 331) fhMCPi0PtOrigin->Fill(pt, 9.5, GetEventWeight());//eta prime
1874  else if(mompdg == 213) fhMCPi0PtOrigin->Fill(pt, 4.5, GetEventWeight());//rho
1875  else if(mompdg == 223) fhMCPi0PtOrigin->Fill(pt, 5.5, GetEventWeight());//omega
1876  else if(mompdg >= 310 && mompdg <= 323)
1877  fhMCPi0PtOrigin->Fill(pt, 6.5, GetEventWeight());//k0S, k+-,k*
1878  else if(mompdg == 130) fhMCPi0PtOrigin->Fill(pt, 6.5, GetEventWeight());//k0L
1879  else if(momstatus == 11 || momstatus == 12 )
1880  fhMCPi0PtOrigin->Fill(pt, 3.5, GetEventWeight());//resonances
1881  else fhMCPi0PtOrigin->Fill(pt, 7.5, GetEventWeight());//other?
1882 
1883  if(status!=11)
1884  {
1885  if (momstatus == 21) fhMCNotResonancePi0PtOrigin->Fill(pt, 0.5, GetEventWeight());//parton
1886  else if(mompdg < 22 ) fhMCNotResonancePi0PtOrigin->Fill(pt, 1.5, GetEventWeight());//quark
1887  else if(mompdg > 2100 && mompdg < 2210)
1888  fhMCNotResonancePi0PtOrigin->Fill(pt, 2.5, GetEventWeight());// resonances
1889  else if(mompdg == 221) fhMCNotResonancePi0PtOrigin->Fill(pt, 8.5, GetEventWeight());//eta
1890  else if(mompdg == 331) fhMCNotResonancePi0PtOrigin->Fill(pt, 9.5, GetEventWeight());//eta prime
1891  else if(mompdg == 213) fhMCNotResonancePi0PtOrigin->Fill(pt, 4.5, GetEventWeight());//rho
1892  else if(mompdg == 223) fhMCNotResonancePi0PtOrigin->Fill(pt, 5.5, GetEventWeight());//omega
1893  else if(mompdg >= 310 && mompdg <= 323)
1894  fhMCNotResonancePi0PtOrigin->Fill(pt, 6.5, GetEventWeight());//k0S, k+-,k*
1895  else if(mompdg == 130) fhMCNotResonancePi0PtOrigin->Fill(pt, 6.5, GetEventWeight());//k0L
1896  else if(momstatus == 11 || momstatus == 12 )
1897  fhMCNotResonancePi0PtOrigin->Fill(pt, 3.5, GetEventWeight());//resonances
1898  else fhMCNotResonancePi0PtOrigin->Fill(pt, 7.5, GetEventWeight());//other?
1899  }
1900  }//pi0 mass region
1901  }
1902  }
1903  else if(ancPDG==221)
1904  {//Eta
1905  mcIndex = 3;
1906  if(fMultiCutAnaSim)
1907  {
1908  for(Int_t ipt=0; ipt<fNPtCuts; ipt++)
1909  {
1910  for(Int_t icell=0; icell<fNCellNCuts; icell++)
1911  {
1912  for(Int_t iasym=0; iasym<fNAsymCuts; iasym++)
1913  {
1914  Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
1915  if(pt1 > fPtCuts[ipt] && pt2 > fPtCuts[ipt] &&
1916  asym < fAsymCuts[iasym] &&
1917  ncell1 >= fCellNCuts[icell] && ncell2 >= fCellNCuts[icell])
1918  {
1919  fhMCEtaMassPtRec [index]->Fill(pt , mass, GetEventWeight());
1920  fhMCEtaMassPtTrue[index]->Fill(fPi0Mom.Pt(), mass, GetEventWeight());
1921  if(mass < 0.65 && mass > 0.45)
1922  fhMCEtaPtTruePtRec[index]->Fill(fPi0Mom.Pt(), pt, GetEventWeight());
1923  }//pass the different cuts
1924  }// pid bit cut loop
1925  }// icell loop
1926  }// pt cut loop
1927  } //Multi cut ana sim
1928  else
1929  {
1930  fhMCEtaMassPtTrue[0]->Fill(fPi0Mom.Pt(), mass, GetEventWeight());
1931  if(mass < 0.65 && mass > 0.45) fhMCEtaPtTruePtRec[0]->Fill(fPi0Mom.Pt(), pt, GetEventWeight());
1932 
1933  if(GetReader()->ReadStack())
1934  {
1935  TParticle* ancestor = GetMCStack()->Particle(ancLabel);
1936  momindex = ancestor->GetFirstMother();
1937  if(momindex < 0) return;
1938  TParticle* mother = GetMCStack()->Particle(momindex);
1939  mompdg = TMath::Abs(mother->GetPdgCode());
1940  momstatus = mother->GetStatusCode();
1941  prodR = mother->R();
1942  }
1943  else
1944  {
1945  TClonesArray * mcparticles = GetReader()->GetAODMCParticles();
1946  AliAODMCParticle* ancestor = (AliAODMCParticle *) mcparticles->At(ancLabel);
1947  momindex = ancestor->GetMother();
1948  if(momindex < 0) return;
1949  AliAODMCParticle* mother = (AliAODMCParticle *) mcparticles->At(momindex);
1950  mompdg = TMath::Abs(mother->GetPdgCode());
1951  momstatus = mother->GetStatus();
1952  prodR = TMath::Sqrt(mother->Xv()*mother->Xv()+mother->Yv()*mother->Yv());
1953  }
1954 
1955  fhMCEtaProdVertex->Fill(pt, prodR, GetEventWeight());
1956 
1957  if (momstatus == 21 ) fhMCEtaPtOrigin->Fill(pt, 0.5, GetEventWeight());//parton
1958  else if(mompdg < 22 ) fhMCEtaPtOrigin->Fill(pt, 1.5, GetEventWeight());//quark
1959  else if(mompdg > 2100 && mompdg < 2210)
1960  fhMCEtaPtOrigin->Fill(pt, 2.5, GetEventWeight());//qq resonances
1961  else if(mompdg == 331) fhMCEtaPtOrigin->Fill(pt, 5.5, GetEventWeight());//eta prime
1962  else if(momstatus == 11 || momstatus == 12 )
1963  fhMCEtaPtOrigin->Fill(pt, 3.5, GetEventWeight());//resonances
1964  else fhMCEtaPtOrigin->Fill(pt, 4.5, GetEventWeight());//stable, conversions?
1965  //printf("Other Meson pdg %d, Mother %s, pdg %d, status %d\n",pdg, TDatabasePDG::Instance()->GetParticle(mompdg)->GetName(),mompdg, momstatus );
1966 
1967  }// eta mass region
1968  }
1969  else if(ancPDG==-2212)
1970  {//AProton
1971  mcIndex = 4;
1972  }
1973  else if(ancPDG==-2112)
1974  {//ANeutron
1975  mcIndex = 5;
1976  }
1977  else if(TMath::Abs(ancPDG)==13)
1978  {//muons
1979  mcIndex = 6;
1980  }
1981  else if (TMath::Abs(ancPDG) > 100 && ancLabel > 7)
1982  {
1983  if(ancStatus==1)
1984  {//Stable particles, converted? not decayed resonances
1985  mcIndex = 6;
1986  }
1987  else
1988  {//resonances and other decays, more hadron conversions?
1989  mcIndex = 7;
1990  }
1991  }
1992  else
1993  {//Partons, colliding protons, strings, intermediate corrections
1994  if(ancStatus==11 || ancStatus==12)
1995  {//String fragmentation
1996  mcIndex = 8;
1997  }
1998  else if (ancStatus==21){
1999  if(ancLabel < 2)
2000  {//Colliding protons
2001  mcIndex = 11;
2002  }//colliding protons
2003  else if(ancLabel < 6)
2004  {//partonic initial states interactions
2005  mcIndex = 9;
2006  }
2007  else if(ancLabel < 8)
2008  {//Final state partons radiations?
2009  mcIndex = 10;
2010  }
2011  // else {
2012  // printf("AliAnaPi0::FillMCVersusRecDataHistograms() - Check ** Common ancestor label %d, pdg %d, name %s, status %d; \n",
2013  // ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus);
2014  // }
2015  }//status 21
2016  //else {
2017  // printf("AliAnaPi0::FillMCVersusRecDataHistograms() - Check *** Common ancestor label %d, pdg %d, name %s, status %d; \n",
2018  // ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus);
2019  // }
2020  }
2021  }//ancLabel > -1
2022  else
2023  { //ancLabel <= -1
2024  //printf("Not related at all label = %d\n",ancLabel);
2025  AliDebug(1,"Common ancestor not found");
2026 
2027  mcIndex = 12;
2028  }
2029 
2030  if(mcIndex >= 0 && mcIndex < 13)
2031  {
2032  fhMCOrgMass [mcIndex]->Fill(pt, mass, GetEventWeight());
2033  fhMCOrgAsym [mcIndex]->Fill(pt, asym, GetEventWeight());
2034  fhMCOrgDeltaEta[mcIndex]->Fill(pt, deta, GetEventWeight());
2035  fhMCOrgDeltaPhi[mcIndex]->Fill(pt, dphi, GetEventWeight());
2036  }
2037 }
2038 
2039 //__________________________________________
2043 //__________________________________________
2044 void AliAnaPi0::MakeAnalysisFillHistograms()
2045 {
2046  // In case of simulated data, fill acceptance histograms
2047  if(IsDataMC())
2048  {
2049  FillAcceptanceHistograms();
2050 
2051  if(fFillOnlyMCAcceptanceHisto) return;
2052  }
2053 
2054 //if (GetReader()->GetEventNumber()%10000 == 0)
2055 // printf("--- Event %d ---\n",GetReader()->GetEventNumber());
2056 
2057  if(!GetInputAODBranch())
2058  {
2059  AliFatal(Form("No input aod photons in AOD with name branch < %s >, STOP",GetInputAODName().Data()));
2060  return;
2061  }
2062 
2063  // Init some variables
2064  Int_t nPhot = GetInputAODBranch()->GetEntriesFast() ;
2065 
2066  AliDebug(1,Form("Photon entries %d", nPhot));
2067 
2068  // If less than photon 2 entries in the list, skip this event
2069  if ( nPhot < 2 )
2070  {
2071  AliDebug(1,Form("nPhotons %d, cent bin %d continue to next event",nPhot, GetEventCentralityBin()));
2072 
2073  if ( GetNCentrBin() > 1 && IsHighMultiplicityAnalysisOn() )
2074  fhCentralityNoPair->Fill(GetEventCentralityBin(), GetEventWeight());
2075 
2076  return ;
2077  }
2078 
2079  Int_t ncentr = GetNCentrBin();
2080  if(GetNCentrBin()==0) ncentr = 1;
2081 
2082  //Init variables
2083  Int_t module1 = -1;
2084  Int_t module2 = -1;
2085  Double_t vert[] = {0.0, 0.0, 0.0} ; //vertex
2086  Int_t evtIndex1 = 0 ;
2087  Int_t currentEvtIndex = -1;
2088  Int_t curCentrBin = GetEventCentralityBin();
2089  //Int_t curVzBin = GetEventVzBin();
2090  //Int_t curRPBin = GetEventRPBin();
2091  Int_t eventbin = GetEventMixBin();
2092 
2093  if(eventbin > GetNCentrBin()*GetNZvertBin()*GetNRPBin())
2094  {
2095  AliWarning(Form("Mix Bin not expected: cen bin %d, z bin %d, rp bin %d, total bin %d, Event Centrality %d, z vertex %2.3f, Reaction Plane %2.3f",
2096  GetEventCentralityBin(),GetEventVzBin(), GetEventRPBin(),eventbin,GetEventCentrality(),vert[2],GetEventPlaneAngle()));
2097  return;
2098  }
2099 
2100  // Get shower shape information of clusters
2101 // TObjArray *clusters = 0;
2102 // if (GetCalorimeter()==kEMCAL) clusters = GetEMCALClusters();
2103 // else if(GetCalorimeter()==kPHOS ) clusters = GetPHOSClusters() ;
2104 
2105  //---------------------------------
2106  // First loop on photons/clusters
2107  //---------------------------------
2108  for(Int_t i1=0; i1<nPhot-1; i1++)
2109  {
2110  AliAODPWG4Particle * p1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i1)) ;
2111 
2112  // Select photons within a pT range
2113  if ( p1->Pt() < GetMinPt() || p1->Pt() > GetMaxPt() ) continue ;
2114 
2115  //printf("AliAnaPi0::MakeAnalysisFillHistograms() : cluster1 id %d/%d\n",i1,nPhot-1);
2116 
2117  // get the event index in the mixed buffer where the photon comes from
2118  // in case of mixing with analysis frame, not own mixing
2119  evtIndex1 = GetEventIndex(p1, vert) ;
2120  if ( evtIndex1 == -1 )
2121  return ;
2122  if ( evtIndex1 == -2 )
2123  continue ;
2124 
2125  // Only effective in case of mixed event frame
2126  if(TMath::Abs(vert[2]) > GetZvertexCut()) continue ; //vertex cut
2127 
2128  if (evtIndex1 != currentEvtIndex)
2129  {
2130  // Fill event bin info
2131  if( DoOwnMix() ) fhEventBin->Fill(eventbin, GetEventWeight()) ;
2132 
2133  if( IsHighMultiplicityAnalysisOn() )
2134  {
2135  fhCentrality->Fill(curCentrBin, GetEventWeight());
2136 
2137  if( GetEventPlane() )
2138  fhEventPlaneResolution->Fill(curCentrBin, TMath::Cos(2.*GetEventPlane()->GetQsubRes()), GetEventWeight());
2139  }
2140 
2141  currentEvtIndex = evtIndex1 ;
2142  }
2143 
2144  //printf("AliAnaPi0::MakeAnalysisFillHistograms(): Photon 1 Evt %d Vertex : %f,%f,%f\n",evtIndex1, GetVertex(evtIndex1)[0] ,GetVertex(evtIndex1)[1],GetVertex(evtIndex1)[2]);
2145 
2146  //Get the momentum of this cluster
2147  fPhotonMom1.SetPxPyPzE(p1->Px(),p1->Py(),p1->Pz(),p1->E());
2148 
2149  //Get (Super)Module number of this cluster
2150  module1 = p1->GetSModNumber();// GetModuleNumber(p1);
2151 
2152  //------------------------------------------
2153  // Recover original cluster
2154 // Int_t iclus1 = -1 ;
2155 // AliVCluster * cluster1 = FindCluster(clusters,p1->GetCaloLabel(0),iclus1);
2156 // if(!cluster1) AliWarning("Cluster1 not found!");
2157 
2158  //---------------------------------
2159  // Second loop on photons/clusters
2160  //---------------------------------
2161  for(Int_t i2=i1+1; i2<nPhot; i2++)
2162  {
2163  AliAODPWG4Particle * p2 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i2)) ;
2164 
2165  // Select photons within a pT range
2166  if ( p2->Pt() < GetMinPt() || p2->Pt() > GetMaxPt() ) continue ;
2167 
2168  //printf("AliAnaPi0::MakeAnalysisFillHistograms() : cluster2 i %d/%d\n",i2,nPhot);
2169 
2170  //In case of mixing frame, check we are not in the same event as the first cluster
2171  Int_t evtIndex2 = GetEventIndex(p2, vert) ;
2172  if ( evtIndex2 == -1 )
2173  return ;
2174  if ( evtIndex2 == -2 )
2175  continue ;
2176  if (GetMixedEvent() && (evtIndex1 == evtIndex2))
2177  continue ;
2178 
2179 // //------------------------------------------
2180 // // Recover original cluster
2181 // Int_t iclus2 = -1;
2182 // AliVCluster * cluster2 = FindCluster(clusters,p2->GetCaloLabel(0),iclus2,iclus1+1);
2183 // // start new loop from iclus1+1 to gain some time
2184 // if(!cluster2) AliWarning("Cluster2 not found!");
2185 //
2186 // // Get the TOF,l0 and ncells from the clusters
2187 // Float_t tof1 = -1;
2188 // Float_t l01 = -1;
2189 // Int_t ncell1 = 0;
2190 // if(cluster1)
2191 // {
2192 // tof1 = cluster1->GetTOF()*1e9;
2193 // l01 = cluster1->GetM02();
2194 // ncell1 = cluster1->GetNCells();
2195 // //printf("cluster1: E %2.2f (%2.2f), l0 %2.2f, tof %2.2f\n",cluster1->E(),p1->E(),l01,tof1);
2196 // }
2197 // //else printf("cluster1 not available: calo label %d / %d, cluster ID %d\n",
2198 // // p1->GetCaloLabel(0),(GetReader()->GetInputEvent())->GetNumberOfCaloClusters()-1,cluster1->GetID());
2199 //
2200 // Float_t tof2 = -1;
2201 // Float_t l02 = -1;
2202 // Int_t ncell2 = 0;
2203 // if(cluster2)
2204 // {
2205 // tof2 = cluster2->GetTOF()*1e9;
2206 // l02 = cluster2->GetM02();
2207 // ncell2 = cluster2->GetNCells();
2208 // //printf("cluster2: E %2.2f (%2.2f), l0 %2.2f, tof %2.2f\n",cluster2->E(),p2->E(),l02,tof2);
2209 // }
2210 // //else printf("cluster2 not available: calo label %d / %d, cluster ID %d\n",
2211 // // p2->GetCaloLabel(0),(GetReader()->GetInputEvent())->GetNumberOfCaloClusters()-1,cluster2->GetID());
2212 //
2213 // if(cluster1 && cluster2)
2214 // {
2215 // Double_t t12diff = tof1-tof2;
2216 // if(TMath::Abs(t12diff) > GetPairTimeCut()) continue;
2217 // }
2218 
2219  Float_t tof1 = p1->GetTime();
2220  Float_t l01 = p1->GetM02();
2221  Int_t ncell1 = p1->GetNCells();
2222  //printf("cluster1: E %2.2f, l0 %2.2f, tof %2.2f\n",p1->E(),l01,tof1);
2223 
2224  Float_t tof2 = p2->GetTime();
2225  Float_t l02 = p2->GetM02();
2226  Int_t ncell2 = p2->GetNCells();
2227  //printf("cluster2: E %2.2f, l0 %2.2f, tof %2.2f\n",p2->E(),l02,tof2);
2228 
2229  Double_t t12diff = tof1-tof2;
2230  fhEPairDiffTime->Fill((fPhotonMom1 + fPhotonMom2).Pt(), t12diff, GetEventWeight());
2231  if(TMath::Abs(t12diff) > GetPairTimeCut()) continue;
2232 
2233  //------------------------------------------
2234 
2235  //printf("AliAnaPi0::MakeAnalysisFillHistograms(): Photon 2 Evt %d Vertex : %f,%f,%f\n",evtIndex2, GetVertex(evtIndex2)[0] ,GetVertex(evtIndex2)[1],GetVertex(evtIndex2)[2]);
2236 
2237  // Get the momentum of this cluster
2238  fPhotonMom2.SetPxPyPzE(p2->Px(),p2->Py(),p2->Pz(),p2->E());
2239 
2240  // Get module number
2241  module2 = p2->GetSModNumber(); //GetModuleNumber(p2);
2242 
2243  //---------------------------------
2244  // Get pair kinematics
2245  //---------------------------------
2246  Double_t m = (fPhotonMom1 + fPhotonMom2).M() ;
2247  Double_t pt = (fPhotonMom1 + fPhotonMom2).Pt();
2248  Double_t deta = fPhotonMom1.Eta() - fPhotonMom2.Eta();
2249  Double_t dphi = fPhotonMom1.Phi() - fPhotonMom2.Phi();
2250  Double_t a = TMath::Abs(p1->E()-p2->E())/(p1->E()+p2->E()) ;
2251 
2252  AliDebug(2,Form("E: fPhotonMom1 %f, fPhotonMom2 %f; Pair: pT %f, mass %f, a %f", p1->E(), p2->E(), (fPhotonMom1 + fPhotonMom2).E(),m,a));
2253 
2254  //--------------------------------
2255  // Opening angle selection
2256  //--------------------------------
2257  // Check if opening angle is too large or too small compared to what is expected
2258  Double_t angle = fPhotonMom1.Angle(fPhotonMom2.Vect());
2259  if(fUseAngleEDepCut && !GetNeutralMesonSelection()->IsAngleInWindow((fPhotonMom1+fPhotonMom2).E(),angle+0.05))
2260  {
2261  AliDebug(2,Form("Real pair angle %f (deg) not in E %f window",RadToDeg(angle), (fPhotonMom1+fPhotonMom2).E()));
2262  continue;
2263  }
2264 
2265  if(fUseAngleCut && (angle < fAngleCut || angle > fAngleMaxCut))
2266  {
2267  AliDebug(2,Form("Real pair cut %f < angle %f < cut %f (deg)",RadToDeg(fAngleCut), RadToDeg(angle), RadToDeg(fAngleMaxCut)));
2268  continue;
2269  }
2270 
2271  //-------------------------------------------------------------------------------------------------
2272  // Fill module dependent histograms, put a cut on assymmetry on the first available cut in the array
2273  //-------------------------------------------------------------------------------------------------
2274  if(a < fAsymCuts[0] && fFillSMCombinations)
2275  {
2276  if(module1==module2 && module1 >=0 && module1<fNModules)
2277  fhReMod[module1]->Fill(pt, m, GetEventWeight()) ;
2278 
2279  if (GetCalorimeter() == kEMCAL )
2280  {
2281  // Same sector
2282  Int_t j=0;
2283  for(Int_t i = 0; i < fNModules/2; i++)
2284  {
2285  j=2*i;
2286  if((module1==j && module2==j+1) || (module1==j+1 && module2==j)) fhReSameSectorEMCALMod[i]->Fill(pt, m, GetEventWeight()) ;
2287  }
2288 
2289  // Same side
2290  for(Int_t i = 0; i < fNModules-2; i++){
2291  if((module1==i && module2==i+2) || (module1==i+2 && module2==i)) fhReSameSideEMCALMod[i]->Fill(pt, m, GetEventWeight());
2292  }
2293  } // EMCAL
2294  else
2295  { // PHOS
2296  if((module1==0 && module2==1) || (module1==1 && module2==0)) fhReDiffPHOSMod[0]->Fill(pt, m, GetEventWeight()) ;
2297  if((module1==0 && module2==2) || (module1==2 && module2==0)) fhReDiffPHOSMod[1]->Fill(pt, m, GetEventWeight()) ;
2298  if((module1==1 && module2==2) || (module1==2 && module2==1)) fhReDiffPHOSMod[2]->Fill(pt, m, GetEventWeight()) ;
2299  } // PHOS
2300  }
2301 
2302  // In case we want only pairs in same (super) module, check their origin.
2303  Bool_t ok = kTRUE;
2304 
2305  if(fSameSM && module1!=module2) ok=kFALSE;
2306 
2307  if(ok)
2308  {
2309  // Check if one of the clusters comes from a conversion
2310  if(fCheckConversion)
2311  {
2312  if (p1->IsTagged() && p2->IsTagged()) fhReConv2->Fill(pt, m, GetEventWeight());
2313  else if(p1->IsTagged() || p2->IsTagged()) fhReConv ->Fill(pt, m, GetEventWeight());
2314  }
2315 
2316  // Fill shower shape cut histograms
2317  if(fFillSSCombinations)
2318  {
2319  if ( l01 > 0.01 && l01 < 0.4 &&
2320  l02 > 0.01 && l02 < 0.4 ) fhReSS[0]->Fill(pt, m, GetEventWeight()); // Tight
2321  else if( l01 > 0.4 && l02 > 0.4 ) fhReSS[1]->Fill(pt, m, GetEventWeight()); // Loose
2322  else if( l01 > 0.01 && l01 < 0.4 && l02 > 0.4 ) fhReSS[2]->Fill(pt, m, GetEventWeight()); // Both
2323  else if( l02 > 0.01 && l02 < 0.4 && l01 > 0.4 ) fhReSS[2]->Fill(pt, m, GetEventWeight()); // Both
2324  }
2325 
2326  // Fill histograms for different bad channel distance, centrality, assymmetry cut and pid bit
2327  for(Int_t ipid=0; ipid<fNPIDBits; ipid++)
2328  {
2329  if((p1->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton)) && (p2->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton)))
2330  {
2331  for(Int_t iasym=0; iasym < fNAsymCuts; iasym++)
2332  {
2333  if(a < fAsymCuts[iasym])
2334  {
2335  Int_t index = ((curCentrBin*fNPIDBits)+ipid)*fNAsymCuts + iasym;
2336  //printf("index %d :(cen %d * nPID %d + ipid %d)*nasym %d + iasym %d - max index %d\n",index,curCentrBin,fNPIDBits,ipid,fNAsymCuts,iasym, curCentrBin*fNPIDBits*fNAsymCuts);
2337 
2338  if(index < 0 || index >= ncentr*fNPIDBits*fNAsymCuts) continue ;
2339 
2340  fhRe1 [index]->Fill(pt, m, GetEventWeight());
2341  if(fMakeInvPtPlots)fhReInvPt1[index]->Fill(pt, m, 1./pt * GetEventWeight()) ;
2342 
2343  if(fFillBadDistHisto)
2344  {
2345  if(p1->DistToBad()>0 && p2->DistToBad()>0)
2346  {
2347  fhRe2 [index]->Fill(pt, m, GetEventWeight()) ;
2348  if(fMakeInvPtPlots)fhReInvPt2[index]->Fill(pt, m, 1./pt * GetEventWeight()) ;
2349 
2350  if(p1->DistToBad()>1 && p2->DistToBad()>1)
2351  {
2352  fhRe3 [index]->Fill(pt, m, GetEventWeight()) ;
2353  if(fMakeInvPtPlots)fhReInvPt3[index]->Fill(pt, m, 1./pt * GetEventWeight()) ;
2354  }// bad 3
2355  }// bad2
2356  }// Fill bad dist histos
2357  }//assymetry cut
2358  }// asymmetry cut loop
2359  }// bad 1
2360  }// pid bit loop
2361 
2362  // Fill histograms with opening angle
2363  if(fFillAngleHisto)
2364  {
2365  fhRealOpeningAngle ->Fill(pt, angle, GetEventWeight());
2366  fhRealCosOpeningAngle->Fill(pt, TMath::Cos(angle), GetEventWeight());
2367  }
2368 
2369  // Fill histograms with pair assymmetry
2370  if(fFillAsymmetryHisto)
2371  {
2372  fhRePtAsym->Fill(pt, a, GetEventWeight());
2373  if(m > 0.10 && m < 0.17) fhRePtAsymPi0->Fill(pt, a, GetEventWeight());
2374  if(m > 0.45 && m < 0.65) fhRePtAsymEta->Fill(pt, a, GetEventWeight());
2375  }
2376 
2377  //---------
2378  // MC data
2379  //---------
2380  // Do some MC checks on the origin of the pair, is there any common ancestor and if there is one, who?
2381  if(IsDataMC())
2382  {
2383  if(GetMCAnalysisUtils()->CheckTagBit(p1->GetTag(),AliMCAnalysisUtils::kMCConversion) &&
2384  GetMCAnalysisUtils()->CheckTagBit(p2->GetTag(),AliMCAnalysisUtils::kMCConversion))
2385  {
2386  fhReMCFromConversion->Fill(pt, m, GetEventWeight());
2387  }
2388  else if(!GetMCAnalysisUtils()->CheckTagBit(p1->GetTag(),AliMCAnalysisUtils::kMCConversion) &&
2389  !GetMCAnalysisUtils()->CheckTagBit(p2->GetTag(),AliMCAnalysisUtils::kMCConversion))
2390  {
2391  fhReMCFromNotConversion->Fill(pt, m, GetEventWeight());
2392  }
2393  else
2394  {
2395  fhReMCFromMixConversion->Fill(pt, m, GetEventWeight());
2396  }
2397 
2398  if(fFillOriginHisto)
2399  FillMCVersusRecDataHistograms(p1->GetLabel(), p2->GetLabel(),p1->Pt(), p2->Pt(),ncell1, ncell2, m, pt, a,deta, dphi);
2400  }
2401 
2402  //-----------------------
2403  // Multi cuts analysis
2404  //-----------------------
2405  if(fMultiCutAna)
2406  {
2407  // Histograms for different PID bits selection
2408  for(Int_t ipid=0; ipid<fNPIDBits; ipid++)
2409  {
2410  if(p1->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton) &&
2411  p2->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton)) fhRePIDBits[ipid]->Fill(pt, m, GetEventWeight()) ;
2412 
2413  //printf("ipt %d, ipid%d, name %s\n",ipt, ipid, fhRePtPIDCuts[ipt*fNPIDBitsBits+ipid]->GetName());
2414  } // pid bit cut loop
2415 
2416  // Several pt,ncell and asymmetry cuts
2417  for(Int_t ipt=0; ipt<fNPtCuts; ipt++)
2418  {
2419  for(Int_t icell=0; icell<fNCellNCuts; icell++)
2420  {
2421  for(Int_t iasym=0; iasym<fNAsymCuts; iasym++)
2422  {
2423  Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
2424  if(p1->E() > fPtCuts[ipt] && p2->E() > fPtCuts[ipt] &&
2425  a < fAsymCuts[iasym] &&
2426  ncell1 >= fCellNCuts[icell] && ncell2 >= fCellNCuts[icell])
2427  {
2428  fhRePtNCellAsymCuts[index]->Fill(pt, m, GetEventWeight()) ;
2429  //printf("ipt %d, icell%d, iasym %d, name %s\n",ipt, icell, iasym, fhRePtNCellAsymCuts[((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym]->GetName());
2430  if(fFillSMCombinations && module1==module2)
2431  {
2432  fhRePtNCellAsymCutsSM[module1][index]->Fill(pt, m, GetEventWeight()) ;
2433  }
2434  }
2435  }// pid bit cut loop
2436  }// icell loop
2437  }// pt cut loop
2438 
2439  if(GetHistogramRanges()->GetHistoTrackMultiplicityBins())
2440  {
2441  for(Int_t iasym = 0; iasym < fNAsymCuts; iasym++)
2442  {
2443  if(a < fAsymCuts[iasym]) fhRePtMult[iasym]->Fill(pt, GetTrackMultiplicity(), m, GetEventWeight()) ;
2444  }
2445  }
2446  }// multiple cuts analysis
2447  }// ok if same sm
2448  }// second same event particle
2449  }// first cluster
2450 
2451  //-------------------------------------------------------------
2452  // Mixing
2453  //-------------------------------------------------------------
2454  if(DoOwnMix())
2455  {
2456  // Recover events in with same characteristics as the current event
2457 
2458  // Check that the bin exists, if not (bad determination of RP, centrality or vz bin) do nothing
2459  if(eventbin < 0) return ;
2460 
2461  TList * evMixList=fEventsList[eventbin] ;
2462 
2463  if(!evMixList)
2464  {
2465  AliWarning(Form("Mix event list not available, bin %d",eventbin));
2466  return;
2467  }
2468 
2469  Int_t nMixed = evMixList->GetSize() ;
2470  for(Int_t ii=0; ii<nMixed; ii++)
2471  {
2472  TClonesArray* ev2= (TClonesArray*) (evMixList->At(ii));
2473  Int_t nPhot2=ev2->GetEntriesFast() ;
2474  Double_t m = -999;
2475  AliDebug(1,Form("Mixed event %d photon entries %d, centrality bin %d",ii, nPhot2, GetEventCentralityBin()));
2476 
2477  fhEventMixBin->Fill(eventbin, GetEventWeight()) ;
2478 
2479  //---------------------------------
2480  // First loop on photons/clusters
2481  //---------------------------------
2482  for(Int_t i1=0; i1<nPhot; i1++)
2483  {
2484  AliAODPWG4Particle * p1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i1)) ;
2485 
2486  if(fSameSM && GetModuleNumber(p1)!=module1) continue;
2487 
2488  //Get kinematics of cluster and (super) module of this cluster
2489  fPhotonMom1.SetPxPyPzE(p1->Px(),p1->Py(),p1->Pz(),p1->E());
2490  module1 = GetModuleNumber(p1);
2491 
2492  //---------------------------------
2493  // Second loop on other mixed event photons/clusters
2494  //---------------------------------
2495  for(Int_t i2=0; i2<nPhot2; i2++)
2496  {
2497  AliAODPWG4Particle * p2 = (AliAODPWG4Particle*) (ev2->At(i2)) ;
2498 
2499  // Get kinematics of second cluster and calculate those of the pair
2500  fPhotonMom2.SetPxPyPzE(p2->Px(),p2->Py(),p2->Pz(),p2->E());
2501  m = (fPhotonMom1+fPhotonMom2).M() ;
2502  Double_t pt = (fPhotonMom1 + fPhotonMom2).Pt();
2503  Double_t a = TMath::Abs(p1->E()-p2->E())/(p1->E()+p2->E()) ;
2504 
2505  // Check if opening angle is too large or too small compared to what is expected
2506  Double_t angle = fPhotonMom1.Angle(fPhotonMom2.Vect());
2507  if(fUseAngleEDepCut && !GetNeutralMesonSelection()->IsAngleInWindow((fPhotonMom1+fPhotonMom2).E(),angle+0.05))
2508  {
2509  AliDebug(2,Form("Mix pair angle %f (deg) not in E %f window",RadToDeg(angle), (fPhotonMom1+fPhotonMom2).E()));
2510  continue;
2511  }
2512 
2513  if(fUseAngleCut && (angle < fAngleCut || angle > fAngleMaxCut))
2514  {
2515  AliDebug(2,Form("Mix pair cut %f < angle %f < cut %f (deg)",RadToDeg(fAngleCut),RadToDeg(angle),RadToDeg(fAngleMaxCut)));
2516  continue;
2517  }
2518 
2519  AliDebug(2,Form("Mixed Event: pT: fPhotonMom1 %2.2f, fPhotonMom2 %2.2f; Pair: pT %2.2f, mass %2.3f, a %2.3f",p1->Pt(), p2->Pt(), pt,m,a));
2520 
2521  // In case we want only pairs in same (super) module, check their origin.
2522  module2 = GetModuleNumber(p2);
2523 
2524  //-------------------------------------------------------------------------------------------------
2525  //Fill module dependent histograms, put a cut on assymmetry on the first available cut in the array
2526  //-------------------------------------------------------------------------------------------------
2527  if(a < fAsymCuts[0] && fFillSMCombinations)
2528  {
2529  if(module1==module2 && module1 >=0 && module1<fNModules)
2530  fhMiMod[module1]->Fill(pt, m, GetEventWeight()) ;
2531 
2532  if(GetCalorimeter()==kEMCAL)
2533  {
2534  // Same sector
2535  Int_t j=0;
2536  for(Int_t i = 0; i < fNModules/2; i++)
2537  {
2538  j=2*i;
2539  if((module1==j && module2==j+1) || (module1==j+1 && module2==j)) fhMiSameSectorEMCALMod[i]->Fill(pt, m, GetEventWeight()) ;
2540  }
2541 
2542  // Same side
2543  for(Int_t i = 0; i < fNModules-2; i++)
2544  {
2545  if((module1==i && module2==i+2) || (module1==i+2 && module2==i)) fhMiSameSideEMCALMod[i]->Fill(pt, m, GetEventWeight());
2546  }
2547  } // EMCAL
2548  else
2549  { // PHOS
2550  if((module1==0 && module2==1) || (module1==1 && module2==0)) fhMiDiffPHOSMod[0]->Fill(pt, m, GetEventWeight()) ;
2551  if((module1==0 && module2==2) || (module1==2 && module2==0)) fhMiDiffPHOSMod[1]->Fill(pt, m, GetEventWeight()) ;
2552  if((module1==1 && module2==2) || (module1==2 && module2==1)) fhMiDiffPHOSMod[2]->Fill(pt, m, GetEventWeight()) ;
2553  } // PHOS
2554  }
2555 
2556  Bool_t ok = kTRUE;
2557  if(fSameSM && module1!=module2) ok=kFALSE;
2558  if(ok){
2559 
2560  // Check if one of the clusters comes from a conversion
2561  if(fCheckConversion)
2562  {
2563  if (p1->IsTagged() && p2->IsTagged()) fhMiConv2->Fill(pt, m, GetEventWeight());
2564  else if(p1->IsTagged() || p2->IsTagged()) fhMiConv ->Fill(pt, m, GetEventWeight());
2565  }
2566 
2567  // Fill histograms for different bad channel distance, centrality, assymmetry cut and pid bit
2568  for(Int_t ipid=0; ipid<fNPIDBits; ipid++)
2569  {
2570  if((p1->IsPIDOK(ipid,AliCaloPID::kPhoton)) && (p2->IsPIDOK(ipid,AliCaloPID::kPhoton)))
2571  {
2572  for(Int_t iasym=0; iasym < fNAsymCuts; iasym++)
2573  {
2574  if(a < fAsymCuts[iasym])
2575  {
2576  Int_t index = ((curCentrBin*fNPIDBits)+ipid)*fNAsymCuts + iasym;
2577 
2578  if(index < 0 || index >= ncentr*fNPIDBits*fNAsymCuts) continue ;
2579 
2580  fhMi1[index]->Fill(pt, m, GetEventWeight()) ;
2581 
2582  if(fMakeInvPtPlots)fhMiInvPt1[index]->Fill(pt, m, 1./pt * GetEventWeight()) ;
2583 
2584  if(fFillBadDistHisto)
2585  {
2586  if(p1->DistToBad()>0 && p2->DistToBad()>0)
2587  {
2588  fhMi2[index]->Fill(pt, m, GetEventWeight()) ;
2589  if(fMakeInvPtPlots)fhMiInvPt2[index]->Fill(pt, m, 1./pt * GetEventWeight()) ;
2590 
2591  if(p1->DistToBad()>1 && p2->DistToBad()>1)
2592  {
2593  fhMi3[index]->Fill(pt, m, GetEventWeight()) ;
2594  if(fMakeInvPtPlots)fhMiInvPt3[index]->Fill(pt, m, 1./pt * GetEventWeight()) ;
2595  }
2596  }
2597  }// Fill bad dist histo
2598  }//Asymmetry cut
2599  }// Asymmetry loop
2600  }//PID cut
2601  }//loop for histograms
2602 
2603  //-----------------------
2604  // Multi cuts analysis
2605  //-----------------------
2606  if(fMultiCutAna)
2607  {
2608  // Several pt,ncell and asymmetry cuts
2609 
2610  for(Int_t ipt=0; ipt<fNPtCuts; ipt++)
2611  {
2612  for(Int_t icell=0; icell<fNCellNCuts; icell++)
2613  {
2614  for(Int_t iasym=0; iasym<fNAsymCuts; iasym++)
2615  {
2616  Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
2617  if(p1->Pt() > fPtCuts[ipt] && p2->Pt() > fPtCuts[ipt] &&
2618  a < fAsymCuts[iasym] //&&
2619  //p1->GetBtag() >= fCellNCuts[icell] && p2->GetBtag() >= fCellNCuts[icell] // trick, correct it.
2620  )
2621  {
2622  fhMiPtNCellAsymCuts[index]->Fill(pt, m, GetEventWeight()) ;
2623  //printf("ipt %d, icell%d, iasym %d, name %s\n",ipt, icell, iasym, fhRePtNCellAsymCuts[((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym]->GetName());
2624  }
2625  }// pid bit cut loop
2626  }// icell loop
2627  }// pt cut loop
2628  } // Multi cut ana
2629 
2630  // Fill histograms with opening angle
2631  if(fFillAngleHisto)
2632  {
2633  fhMixedOpeningAngle ->Fill(pt, angle, GetEventWeight());
2634  fhMixedCosOpeningAngle->Fill(pt, TMath::Cos(angle), GetEventWeight());
2635  }
2636  }//ok
2637  }// second cluster loop
2638  }//first cluster loop
2639  }//loop on mixed events
2640 
2641  //--------------------------------------------------------
2642  // Add the current event to the list of events for mixing
2643  //--------------------------------------------------------
2644 
2645  TClonesArray *currentEvent = new TClonesArray(*GetInputAODBranch());
2646  // Add current event to buffer and Remove redundant events
2647  if( currentEvent->GetEntriesFast() > 0 )
2648  {
2649  evMixList->AddFirst(currentEvent) ;
2650  currentEvent=0 ; //Now list of particles belongs to buffer and it will be deleted with buffer
2651  if( evMixList->GetSize() >= GetNMaxEvMix() )
2652  {
2653  TClonesArray * tmp = (TClonesArray*) (evMixList->Last()) ;
2654  evMixList->RemoveLast() ;
2655  delete tmp ;
2656  }
2657  }
2658  else
2659  { // empty event
2660  delete currentEvent ;
2661  currentEvent=0 ;
2662  }
2663  }// DoOwnMix
2664 
2665  AliDebug(1,"End fill histograms");
2666 }
2667 
2668 //________________________________________________________________________
2672 //________________________________________________________________________
2673 Int_t AliAnaPi0::GetEventIndex(AliAODPWG4Particle * part, Double_t * vert)
2674 {
2675  Int_t evtIndex = -1 ;
2676 
2677  if(GetReader()->GetDataType()!=AliCaloTrackReader::kMC)
2678  {
2679  if (GetMixedEvent())
2680  {
2681  evtIndex = GetMixedEvent()->EventIndexForCaloCluster(part->GetCaloLabel(0)) ;
2682  GetVertex(vert,evtIndex);
2683 
2684  if(TMath::Abs(vert[2])> GetZvertexCut())
2685  evtIndex = -2 ; //Event can not be used (vertex, centrality,... cuts not fulfilled)
2686  }
2687  else
2688  {
2689  // Single event
2690  GetVertex(vert);
2691 
2692  if(TMath::Abs(vert[2])> GetZvertexCut())
2693  evtIndex = -1 ; //Event can not be used (vertex, centrality,... cuts not fulfilled)
2694  else
2695  evtIndex = 0 ;
2696  }
2697  } // No MC reader
2698  else
2699  {
2700  evtIndex = 0;
2701  vert[0] = 0. ;
2702  vert[1] = 0. ;
2703  vert[2] = 0. ;
2704  }
2705 
2706  return evtIndex ;
2707 }
2708 
AliAnaPi0::fhPrimEtaY
TH2F * fhPrimEtaY
! Rapidity distribution of primary particles vs pT
Definition: AliAnaPi0.h:359
AliHistogramRanges::GetHistoPtMax
Float_t GetHistoPtMax() const
Definition: AliHistogramRanges.h:39
AliAnaPi0::fhRe3
TH2F ** fhRe3
[GetNCentrBin()*fNPIDBits*fNAsymCuts]
Definition: AliAnaPi0.h:260
AliAnaPi0::fhPrimEtaAccPtCentrality
TH2F * fhPrimEtaAccPtCentrality
! primary eta with accepted daughters reconstructed centrality vs pT
Definition: AliAnaPi0.h:371
pdg
Int_t pdg
Definition: charmCutsOptimization.C:34
ClassImp
ClassImp(AliAnalysisTaskTriggerRates) AliAnalysisTaskTriggerRates
Definition: AliAnalysisTaskTriggerRates.cxx:34
AliAnaPi0::fhPrimPi0Y
TH2F * fhPrimPi0Y
! Rapidity distribution of primary particles vs pT
Definition: AliAnaPi0.h:338
AliAnaPi0::fPhotonMom1Boost
TLorentzVector fPhotonMom1Boost
! Photon cluster momentum, temporary array
Definition: AliAnaPi0.h:195
AliAnaPi0::fNModules
Int_t fNModules
[GetNCentrBin()*GetNZvertBin()*GetNRPBin()]
Definition: AliAnaPi0.h:160
AliAnaPi0::fhPrimEtaYetaYcut
TH2F * fhPrimEtaYetaYcut
! PseudoRapidity distribution of primary particles vs pT, Y<1
Definition: AliAnaPi0.h:362
AliAnaPi0::fhPrimPi0AccPt
TH1F * fhPrimPi0AccPt
! Spectrum of primary with accepted daughters
Definition: AliAnaPi0.h:337
AliHistogramRanges::GetHistoPtMin
Float_t GetHistoPtMin() const
Definition: AliHistogramRanges.h:38
AliAnaPi0::fhPrimEtaAccPhi
TH2F * fhPrimEtaAccPhi
! Azimutal distribution of primary with accepted daughters vs pT
Definition: AliAnaPi0.h:365
AliAnaPi0::fPi0Mom
TLorentzVector fPi0Mom
! Pi0 cluster momentum, temporary array
Definition: AliAnaPi0.h:197
AliAnaPi0::fhReMCFromMixConversion
TH2F * fhReMCFromMixConversion
! Invariant mass of 2 clusters one from conversion and the other not
Definition: AliAnaPi0.h:422
AliCaloTrackReader.h
AliAnaPi0::fPIDBits
Int_t fPIDBits[10]
Array with different PID bits.
Definition: AliAnaPi0.h:177
AliAnaCaloTrackCorrBaseClass::AddToHistogramsName
virtual void AddToHistogramsName(TString add)
Definition: AliAnaCaloTrackCorrBaseClass.h:85
AliAnaPi0::fhReConv
TH2F * fhReConv
[fNModules-1]
Definition: AliAnaPi0.h:239
AliAnaPi0::fFillArmenterosThetaStar
Bool_t fFillArmenterosThetaStar
Fill armenteros histograms.
Definition: AliAnaPi0.h:189
AliAnaPi0::fNPIDBits
Int_t fNPIDBits
Number of possible PID bit combinations.
Definition: AliAnaPi0.h:176
AliAnaPi0::fPhotonMom2
TLorentzVector fPhotonMom2
! Photon cluster momentum, temporary array
Definition: AliAnaPi0.h:196
title
const char * title
Definition: MakeQAPdf.C:26
AliAnaPi0::fCellNCuts
Int_t fCellNCuts[10]
Array with different cell number cluster cuts.
Definition: AliAnaPi0.h:175
AliAnaPi0::fhMCOrgDeltaEta
TH2F * fhMCOrgDeltaEta[13]
! Delta Eta vs pt of real pairs, check common origin of pair
Definition: AliAnaPi0.h:387
AliCaloPID.h
AliAnaPi0::fhReSS
TH2F * fhReSS[3]
[fNAsymCuts]
Definition: AliAnaPi0.h:307
AliAnaPi0::fhPrimEtaCosOpeningAngle
TH2F * fhPrimEtaCosOpeningAngle
! Cosinus of opening angle of pair version pair energy, eta primaries
Definition: AliAnaPi0.h:368
AliAnaPi0::fhPrimPi0PtStatus
TH2F * fhPrimPi0PtStatus
! Spectrum of generated pi0 vs pi0 status
Definition: AliAnaPi0.h:379
AliAnaPi0::fCheckConversion
Bool_t fCheckConversion
Fill histograms with tagged photons as conversion.
Definition: AliAnaPi0.h:183
AliAnaCaloTrackCorrBaseClass::GetPairTimeCut
virtual Float_t GetPairTimeCut() const
Time cut in ns.
Definition: AliAnaCaloTrackCorrBaseClass.h:221
AliAnaCaloTrackCorrBaseClass::GetVertex
virtual void GetVertex(Double_t vertex[3]) const
Definition: AliAnaCaloTrackCorrBaseClass.h:271
AliHistogramRanges::GetHistoTrackMultiplicityMax
Int_t GetHistoTrackMultiplicityMax() const
Definition: AliHistogramRanges.h:133
AliAnaPi0::fhPrimPi0E
TH1F * fhPrimPi0E
! Spectrum of Primary
Definition: AliAnaPi0.h:334
AliAnaPi0::fUseAngleEDepCut
Bool_t fUseAngleEDepCut
Select pairs depending on their opening angle.
Definition: AliAnaPi0.h:163
AliAnaPi0::fhMiSameSectorEMCALMod
TH2F ** fhMiSameSectorEMCALMod
[fNModules-2]
Definition: AliAnaPi0.h:232
AliAnaPi0
Selected photon clusters invariant mass analysis.
Definition: AliAnaPi0.h:38
AliAnaPi0::fhCentrality
TH1F * fhCentrality
! Histogram with centrality bins with at least one pare
Definition: AliAnaPi0.h:319
AliAnaCaloTrackCorrBaseClass::SetInputAODName
virtual void SetInputAODName(TString name)
Definition: AliAnaCaloTrackCorrBaseClass.h:124
mass
Double_t mass
Definition: charmCutsOptimization.C:35
AliAnaCaloTrackCorrBaseClass::DegToRad
Float_t DegToRad(Float_t deg) const
Definition: AliAnaCaloTrackCorrBaseClass.h:289
AliAnaPi0::fhPrimPi0PtCentrality
TH2F * fhPrimPi0PtCentrality
! primary pi0 reconstructed centrality vs pT
Definition: AliAnaPi0.h:348
AliAnaPi0::fhPrimEtaPtCentrality
TH2F * fhPrimEtaPtCentrality
! primary eta reconstructed centrality vs pT
Definition: AliAnaPi0.h:369
AliAnaCaloTrackCorrBaseClass::GetZvertexCut
virtual Float_t GetZvertexCut() const
Maximal number of events for mixin.
Definition: AliAnaCaloTrackCorrBaseClass.h:235
AliAnaCaloTrackCorrBaseClass::GetEventPlaneAngle
virtual Double_t GetEventPlaneAngle() const
Definition: AliAnaCaloTrackCorrBaseClass.h:112
AliAnaPi0::fhReMCFromConversion
TH2F * fhReMCFromConversion
! Invariant mass of 2 clusters originated in conversions
Definition: AliAnaPi0.h:420
AliAnaPi0::fhReSameSideEMCALMod
TH2F ** fhReSameSideEMCALMod
[fNModules]
Definition: AliAnaPi0.h:217
AliAnaPi0::fhPrimPi0OpeningAngle
TH2F * fhPrimPi0OpeningAngle
! Opening angle of pair versus pair energy, primaries
Definition: AliAnaPi0.h:345
AliAnaCaloTrackCorrBaseClass::GetNeutralMesonSelection
virtual AliNeutralMesonSelection * GetNeutralMesonSelection()
Definition: AliAnaCaloTrackCorrBaseClass.h:332
AliAnaPi0::fhRealOpeningAngle
TH2F * fhRealOpeningAngle
! Opening angle of pair versus pair energy
Definition: AliAnaPi0.h:326
AliAnaPi0::fhMCOrgMass
TH2F * fhMCOrgMass[13]
! Mass vs pt of real pairs, check common origin of pair
Definition: AliAnaPi0.h:385
AliHistogramRanges::GetHistoMassBins
Int_t GetHistoMassBins() const
Definition: AliHistogramRanges.h:97
AliHistogramRanges::GetHistoPhiBins
Int_t GetHistoPhiBins() const
Definition: AliHistogramRanges.h:60
AliAnaPi0::fhRePtAsymPi0
TH2F * fhRePtAsymPi0
! REAL two-photon pt vs asymmetry, close to pi0 mass
Definition: AliAnaPi0.h:312
AliFiducialCut.h
AliAnaPi0::fPhotonMom1
TLorentzVector fPhotonMom1
! Photon cluster momentum, temporary array
Definition: AliAnaPi0.h:194
AliAnaPi0::fhRealCosOpeningAngle
TH2F * fhRealCosOpeningAngle
! Cosinus of opening angle of pair version pair energy
Definition: AliAnaPi0.h:327
AliAnaPi0::fhRe1
TH2F ** fhRe1
REAL two-photon invariant mass distribution for different centralities and Asymmetry.
Definition: AliAnaPi0.h:245
AliHistogramRanges::GetHistoMassMin
Float_t GetHistoMassMin() const
Definition: AliHistogramRanges.h:98
AliAnaPi0::fhPrimPi0ProdVertex
TH2F * fhPrimPi0ProdVertex
! Spectrum of primary pi0 vs production vertex
Definition: AliAnaPi0.h:417
AliAnaPi0::fhReDiffPHOSMod
TH2F ** fhReDiffPHOSMod
[fNModules/2]
Definition: AliAnaPi0.h:223
AliAnaCaloTrackCorrBaseClass::RadToDeg
Float_t RadToDeg(Float_t rad) const
Definition: AliAnaCaloTrackCorrBaseClass.h:291
AliAnaPi0::fhMiConv
TH2F * fhMiConv
! MIXED two-photon invariant mass distribution one of the pair was 2 clusters with small mass ...
Definition: AliAnaPi0.h:240
AliHistogramRanges::GetHistoTrackMultiplicityMin
Int_t GetHistoTrackMultiplicityMin() const
Definition: AliHistogramRanges.h:132
AliAnaPi0::fAsymCuts
Float_t fAsymCuts[10]
Array with different assymetry cuts.
Definition: AliAnaPi0.h:173
AliAnaPi0::fPtCuts
Float_t fPtCuts[10]
Array with different pt cuts.
Definition: AliAnaPi0.h:171
AliAnaPi0::fhMCOrgAsym
TH2F * fhMCOrgAsym[13]
! Asymmetry vs pt of real pairs, check common origin of pair
Definition: AliAnaPi0.h:386
AliAnaPi0::fhRePtNCellAsymCuts
TH2F ** fhRePtNCellAsymCuts
[GetNCentrBin()*fNPIDBits*fNAsymCuts]
Definition: AliAnaPi0.h:293
AliAnaPi0::fhPrimPi0PtOrigin
TH2F * fhPrimPi0PtOrigin
! Spectrum of generated pi0 vs mother
Definition: AliAnaPi0.h:376
AliHistogramRanges::GetHistoDiffTimeMin
Float_t GetHistoDiffTimeMin() const
Definition: AliHistogramRanges.h:278
AliCalorimeterUtils::IsMCParticleInCalorimeterAcceptance
Bool_t IsMCParticleInCalorimeterAcceptance(Int_t calo, TParticle *particle)
Definition: AliCalorimeterUtils.cxx:1405
AliAnaPi0::fhMCPi0MassPtRec
TH2F ** fhMCPi0MassPtRec
Real pi0 pairs, reconstructed mass vs reconstructed pt of original pair.
Definition: AliAnaPi0.h:393
AliAnaPi0::fhPrimPi0AccPtCentrality
TH2F * fhPrimPi0AccPtCentrality
! primary pi0 with accepted daughters reconstructed centrality vs pT
Definition: AliAnaPi0.h:350
AliHistogramRanges::GetHistoPhiMin
Float_t GetHistoPhiMin() const
Definition: AliHistogramRanges.h:61
AliAnaPi0::fNAsymCuts
Int_t fNAsymCuts
Number of assymmetry cuts.
Definition: AliAnaPi0.h:172
AliHistogramRanges::GetHistoDiffTimeMax
Float_t GetHistoDiffTimeMax() const
Definition: AliHistogramRanges.h:279
AliAnaPi0::fhMixedCosOpeningAngle
TH2F * fhMixedCosOpeningAngle
! Cosinus of opening angle of pair version pair energy
Definition: AliAnaPi0.h:329
AliAnaPi0::fProdVertex
TVector3 fProdVertex
! Production vertex, temporary array
Definition: AliAnaPi0.h:198
AliAnaPi0::fhReInvPt2
TH2F ** fhReInvPt2
[GetNCentrBin()*fNPIDBits*fNAsymCuts]
Definition: AliAnaPi0.h:276
AliAnaPi0::fhMiDiffPHOSMod
TH2F ** fhMiDiffPHOSMod
[fNModules/2]
Definition: AliAnaPi0.h:235
AliAnaPi0::fhMi2
TH2F ** fhMi2
[GetNCentrBin()*fNPIDBits*fNAsymCuts]
Definition: AliAnaPi0.h:256
AliAnaPi0::GetCreateOutputObjects
TList * GetCreateOutputObjects()
Definition: AliAnaPi0.cxx:235
AliAnaPi0::fhReInvPt3
TH2F ** fhReInvPt3
[GetNCentrBin()*fNPIDBits*fNAsymCuts]
Definition: AliAnaPi0.h:284
AliAnaPi0::fFillOriginHisto
Bool_t fFillOriginHisto
Fill histograms depending on their origin.
Definition: AliAnaPi0.h:188
AliAnaPi0::fhPrimPi0Phi
TH2F * fhPrimPi0Phi
! Azimutal distribution of primary particles vs pT
Definition: AliAnaPi0.h:343
AliNeutralMesonSelection::IsAngleInWindow
Bool_t IsAngleInWindow(Float_t angle, Float_t e) const
Definition: AliNeutralMesonSelection.cxx:194
etamin
const Double_t etamin
Definition: AddTaskCFDStar.C:40
AliHistogramRanges::GetHistoAsymmetryMax
Float_t GetHistoAsymmetryMax() const
Definition: AliHistogramRanges.h:108
AliHistogramRanges::GetHistoMassMax
Float_t GetHistoMassMax() const
Definition: AliHistogramRanges.h:99
AliAnaCaloTrackCorrBaseClass
Base class for CaloTrackCorr analysis algorithms.
Definition: AliAnaCaloTrackCorrBaseClass.h:58
AliNeutralMesonSelection.h
AliAnaPi0::fFillAngleHisto
Bool_t fFillAngleHisto
Fill histograms with pair opening angle.
Definition: AliAnaPi0.h:186
AliAnaCaloTrackCorrBaseClass::GetCalorimeterString
virtual TString GetCalorimeterString() const
Definition: AliAnaCaloTrackCorrBaseClass.h:169
AliAnaPi0::fFillBadDistHisto
Bool_t fFillBadDistHisto
Do plots for different distances to bad channels.
Definition: AliAnaPi0.h:184
AliAnaCaloTrackCorrBaseClass::IsRealCaloAcceptanceOn
virtual Bool_t IsRealCaloAcceptanceOn() const
Definition: AliAnaCaloTrackCorrBaseClass.h:183
AliAnaCaloTrackCorrBaseClass::GetFiducialCut
virtual AliFiducialCut * GetFiducialCut()
Definition: AliAnaCaloTrackCorrBaseClass.h:324
AliAnaCaloTrackCorrBaseClass::GetInputAODBranch
virtual TClonesArray * GetInputAODBranch() const
Definition: AliAnaCaloTrackCorrBaseClass.h:137
AliAnaPi0::fhMCEtaPtTruePtRec
TH2F ** fhMCEtaPtTruePtRec
[fNPtCuts*fNAsymCuts*fNCellNCuts]
Definition: AliAnaPi0.h:408
AliAnaPi0::fMultiCutAna
Bool_t fMultiCutAna
Do analysis with several or fixed cut.
Definition: AliAnaPi0.h:168
AliCaloTrackReader::GetAODMCParticles
virtual TClonesArray * GetAODMCParticles() const
Definition: AliCaloTrackReader.cxx:699
AliAnaCaloTrackCorrBaseClass::GetHistogramRanges
virtual AliHistogramRanges * GetHistogramRanges()
Definition: AliAnaCaloTrackCorrBaseClass.h:326
AliAnaPi0::fhReMod
TH2F ** fhReMod
REAL two-photon invariant mass distribution for different calorimeter modules.
Definition: AliAnaPi0.h:214
AliCaloTrackReader::ReadStack
Bool_t ReadStack() const
Definition: AliCaloTrackReader.h:603
AliAnaPi0::fhCentralityNoPair
TH1F * fhCentralityNoPair
! Histogram with centrality bins with no pair
Definition: AliAnaPi0.h:320
AliAnaPi0::fFillAsymmetryHisto
Bool_t fFillAsymmetryHisto
Fill histograms with asymmetry vs pt.
Definition: AliAnaPi0.h:187
AliHistogramRanges::GetHistoDiffTimeBins
Int_t GetHistoDiffTimeBins() const
Definition: AliHistogramRanges.h:277
AliAnaPi0::GetEventIndex
Int_t GetEventIndex(AliAODPWG4Particle *part, Double_t *vert)
Definition: AliAnaPi0.cxx:2673
AliAnaPi0::FillAcceptanceHistograms
void FillAcceptanceHistograms()
Fill acceptance histograms if MC data is available.
Definition: AliAnaPi0.cxx:1312
AliAnaPi0::fhReSameSectorEMCALMod
TH2F ** fhReSameSectorEMCALMod
[fNModules-2]
Definition: AliAnaPi0.h:220
AliAnaPi0::fhMCEtaMassPtTrue
TH2F ** fhMCEtaMassPtTrue
[fNPtCuts*fNAsymCuts*fNCellNCuts]
Definition: AliAnaPi0.h:405
AliAnaPi0::fhReInvPt1
TH2F ** fhReInvPt1
[GetNCentrBin()*fNPIDBits*fNAsymCuts]
Definition: AliAnaPi0.h:269
AliAnaPi0::fNCellNCuts
Int_t fNCellNCuts
Number of cuts with number of cells in cluster.
Definition: AliAnaPi0.h:174
AliAnaPi0::fhMiInvPt1
TH2F ** fhMiInvPt1
[GetNCentrBin()*fNPIDBits*fNAsymCuts]
Definition: AliAnaPi0.h:272
AliHistogramRanges::GetHistoAsymmetryBins
Int_t GetHistoAsymmetryBins() const
Definition: AliHistogramRanges.h:106
ptmax
const Double_t ptmax
Definition: AddTaskCFDStar.C:39
AliAnaPi0::fhPrimEtaAccYeta
TH2F * fhPrimEtaAccYeta
! PseudoRapidity distribution of primary with accepted daughters vs pT
Definition: AliAnaPi0.h:363
AliAnaPi0::fAngleMaxCut
Float_t fAngleMaxCut
Select pairs with opening angle smaller than a threshold.
Definition: AliAnaPi0.h:165
AliAnaPi0::fhEventMixBin
TH1I * fhEventMixBin
! Number of mixed pairs in a particular bin (cen,vz,rp)
Definition: AliAnaPi0.h:318
AliAnaPi0::FillArmenterosThetaStar
void FillArmenterosThetaStar(Int_t pdg)
Fill armenteros plots.
Definition: AliAnaPi0.cxx:1708
AliAnaCaloTrackCorrBaseClass::GetEMCALGeometry
virtual AliEMCALGeometry * GetEMCALGeometry() const
Definition: AliAnaCaloTrackCorrBaseClass.h:336
AliAnaPi0::fhPrimPi0CosOpeningAngle
TH2F * fhPrimPi0CosOpeningAngle
! Cosinus of opening angle of pair version pair energy, pi0 primaries
Definition: AliAnaPi0.h:347
AliAnaPi0::fhMCPi0PtOrigin
TH2F * fhMCPi0PtOrigin
[fNPtCuts*fNAsymCuts*fNCellNCuts]
Definition: AliAnaPi0.h:410
AliAnaPi0::fhMiSameSideEMCALMod
TH2F ** fhMiSameSideEMCALMod
[fNModules]
Definition: AliAnaPi0.h:229
AliAnaPi0::fFillOnlyMCAcceptanceHisto
Bool_t fFillOnlyMCAcceptanceHisto
Do analysis only of MC kinematics input.
Definition: AliAnaPi0.h:190
AliFiducialCut::IsInFiducialCut
Bool_t IsInFiducialCut(Float_t eta, Float_t phi, Int_t det) const
Definition: AliFiducialCut.cxx:78
AliAnaPi0::fMakeInvPtPlots
Bool_t fMakeInvPtPlots
Do plots with inverse pt weight.
Definition: AliAnaPi0.h:180
AliAnaPi0.h
AliAnaPi0::fhPrimEtaPtEventPlane
TH2F * fhPrimEtaPtEventPlane
! primary eta reconstructed event plane vs pT
Definition: AliAnaPi0.h:370
AliAnaPi0::fhPrimPi0AccE
TH1F * fhPrimPi0AccE
! Spectrum of primary with accepted daughters
Definition: AliAnaPi0.h:336
AliAnaPi0::fNPtCuts
Int_t fNPtCuts
Number of pt cuts.
Definition: AliAnaPi0.h:170
AliAnaPi0::fhPrimEtaE
TH1F * fhPrimEtaE
! Spectrum of Primary
Definition: AliAnaPi0.h:355
AliAnaPi0::fhPrimEtaAccPtEventPlane
TH2F * fhPrimEtaAccPtEventPlane
! primary eta with accepted daughters reconstructed event plane vs pT
Definition: AliAnaPi0.h:372
AliAnaCaloTrackCorrBaseClass::GetCaloUtils
virtual AliCalorimeterUtils * GetCaloUtils() const
Definition: AliAnaCaloTrackCorrBaseClass.h:322
AliAnaPi0::fhPrimNotResonancePi0PtOrigin
TH2F * fhPrimNotResonancePi0PtOrigin
! Spectrum of generated pi0 vs mother
Definition: AliAnaPi0.h:378
AliAnaPi0::fhPrimEtaAccY
TH2F * fhPrimEtaAccY
! Rapidity distribution of primary with accepted daughters vs pT
Definition: AliAnaPi0.h:360
AliAnaPi0::fhPrimEtaPt
TH1F * fhPrimEtaPt
! Spectrum of Primary
Definition: AliAnaPi0.h:356
AliAnaPi0::fMultiCutAnaSim
Bool_t fMultiCutAnaSim
Do analysis with several or fixed cut, in the simulation related part.
Definition: AliAnaPi0.h:169
AliAnaPi0::fhPrimEtaAccPt
TH1F * fhPrimEtaAccPt
! Spectrum of primary with accepted daughters
Definition: AliAnaPi0.h:358
AliAnaPi0::fhArmPrimPi0
TH2F * fhArmPrimPi0[4]
! Armenteros plots for primary pi0 in 6 energy bins
Definition: AliAnaPi0.h:424
AliAnaPi0::Print
void Print(const Option_t *opt) const
Print some relevant parameters set for the analysis.
Definition: AliAnaPi0.cxx:1269
AliAnaPi0::fhPrimPi0AccPtEventPlane
TH2F * fhPrimPi0AccPtEventPlane
! primary pi0 with accepted daughters reconstructed event plane vs pT
Definition: AliAnaPi0.h:351
AliAnaCaloTrackCorrBaseClass::GetEventPlane
virtual AliEventplane * GetEventPlane() const
Definition: AliAnaCaloTrackCorrBaseClass.h:111
AliCalorimeterUtils::GetNumberOfSuperModulesUsed
Int_t GetNumberOfSuperModulesUsed() const
Definition: AliCalorimeterUtils.h:350
ptmin
const Double_t ptmin
Definition: AddTaskCFDStar.C:38
AliAnaCaloTrackCorrBaseClass::IsHighMultiplicityAnalysisOn
virtual Bool_t IsHighMultiplicityAnalysisOn() const
Definition: AliAnaCaloTrackCorrBaseClass.h:200
AliAnaPi0::fhPrimEtaOpeningAngle
TH2F * fhPrimEtaOpeningAngle
! Opening angle of pair versus pair energy, eta primaries
Definition: AliAnaPi0.h:366
AliAnaPi0::fhMiPtNCellAsymCuts
TH2F ** fhMiPtNCellAsymCuts
[fNPtCuts*fNAsymCuts*fNCellNCuts]
Definition: AliAnaPi0.h:296
AliAnaPi0::fhPrimPi0OpeningAngleAsym
TH2F * fhPrimPi0OpeningAngleAsym
! Opening angle of pair versus pair E asymmetry, pi0 primaries
Definition: AliAnaPi0.h:346
AliAnaPi0::fhPrimPi0AccY
TH2F * fhPrimPi0AccY
! Rapidity distribution of primary with accepted daughters vs pT
Definition: AliAnaPi0.h:339
AliAnaPi0::fhPrimPi0AccPhi
TH2F * fhPrimPi0AccPhi
! Azimutal distribution of primary with accepted daughters vs pT
Definition: AliAnaPi0.h:344
AliAnaPi0::fAngleCut
Float_t fAngleCut
Select pairs with opening angle larger than a threshold.
Definition: AliAnaPi0.h:164
AliAnaCaloTrackCorrBaseClass::GetEventWeight
virtual Double_t GetEventWeight() const
Definition: AliAnaCaloTrackCorrBaseClass.h:243
AliAnaPi0::fhPrimPi0AccYeta
TH2F * fhPrimPi0AccYeta
! PseudoRapidity distribution of primary with accepted daughters vs pT
Definition: AliAnaPi0.h:342
AliAnaCaloTrackCorrBaseClass::GetNCentrBin
virtual Int_t GetNCentrBin() const
Number of bins in reaction plain.
Definition: AliAnaCaloTrackCorrBaseClass.h:232
AliHistogramRanges::GetHistoTrackMultiplicityBins
Int_t GetHistoTrackMultiplicityBins() const
Definition: AliHistogramRanges.h:131
AliAnaPi0::fhMi1
TH2F ** fhMi1
[GetNCentrBin()*fNPIDBits*fNAsymCuts]
Definition: AliAnaPi0.h:248
AliAnaPi0::fhCosThStarPrimEta
TH2F * fhCosThStarPrimEta
! cos(theta*) plots vs E for primary eta, same as asymmetry ...
Definition: AliAnaPi0.h:427
AliAnaPi0::GetAnalysisCuts
TObjString * GetAnalysisCuts()
Save parameters used for analysis.
Definition: AliAnaPi0.cxx:187
AliAnaPi0::MakeAnalysisFillHistograms
void MakeAnalysisFillHistograms()
Definition: AliAnaPi0.cxx:2044
AliHistogramRanges::GetHistoEtaMin
Float_t GetHistoEtaMin() const
Definition: AliHistogramRanges.h:89
AliAnaPi0::fFillSMCombinations
Bool_t fFillSMCombinations
Fill histograms with different cluster pairs in SM combinations.
Definition: AliAnaPi0.h:182
Data
Bool_t Data(TH1F *h, Double_t *rangefit, Bool_t writefit, Double_t &sgn, Double_t &errsgn, Double_t &bkg, Double_t &errbkg, Double_t &sgnf, Double_t &errsgnf, Double_t &sigmafit, Int_t &status)
Definition: charmCutsOptimization.C:371
AliAnaPi0::fhReMCFromNotConversion
TH2F * fhReMCFromNotConversion
! Invariant mass of 2 clusters not originated in conversions
Definition: AliAnaPi0.h:421
AliAnaPi0::fhPrimPi0YetaYcut
TH2F * fhPrimPi0YetaYcut
! PseudoRapidity distribution of primary particles vs pT, Y<1
Definition: AliAnaPi0.h:341
AliAnaCaloTrackCorrBaseClass::GetNMaxEvMix
virtual Int_t GetNMaxEvMix() const
Number of bins in track multiplicity.
Definition: AliAnaCaloTrackCorrBaseClass.h:234
AliAnaPi0::fhPrimPi0Yeta
TH2F * fhPrimPi0Yeta
! PseudoRapidity distribution of primary particles vs pT
Definition: AliAnaPi0.h:340
AliAnaPi0::fhRePtMult
TH3F ** fhRePtMult
[fNPIDBits]
Definition: AliAnaPi0.h:305
AliAnaCaloTrackCorrBaseClass::GetModuleNumber
virtual Int_t GetModuleNumber(AliAODPWG4Particle *part) const
Definition: AliAnaCaloTrackCorrBaseClass.h:301
AliAnaPi0::fhMCNotResonancePi0PtOrigin
TH2F * fhMCNotResonancePi0PtOrigin
! Mass of reconstructed pi0 pairs in calorimeter vs mother origin ID, pi0 status 1.
Definition: AliAnaPi0.h:412
AliAnaPi0::fhRePtAsym
TH2F * fhRePtAsym
! REAL two-photon pt vs asymmetry
Definition: AliAnaPi0.h:311
AliHistogramRanges::GetHistoEtaMax
Float_t GetHistoEtaMax() const
Definition: AliHistogramRanges.h:90
AliAnaPi0::fhEPairDiffTime
TH2F * fhEPairDiffTime
! E pair vs Pair of clusters time difference vs E
Definition: AliAnaPi0.h:429
AliAnaPi0::InitParameters
void InitParameters()
Definition: AliAnaPi0.cxx:153
AliAnaPi0::fhMCEtaPtOrigin
TH2F * fhMCEtaPtOrigin
! Mass of reconstructed eta pairs in calorimeter vs mother origin ID.
Definition: AliAnaPi0.h:411
AliAnaPi0::fhMCPi0PtStatus
TH2F * fhMCPi0PtStatus
! Mass of reconstructed pi0 pairs in calorimeter vs mother status.
Definition: AliAnaPi0.h:413
AliHistogramRanges::GetHistoPtBins
Int_t GetHistoPtBins() const
Definition: AliHistogramRanges.h:37
AliAnaPi0::fhMiConv2
TH2F * fhMiConv2
! MIXED two-photon invariant mass distribution both pair photons recombined from 2 clusters with smal...
Definition: AliAnaPi0.h:242
AliAnaPi0::~AliAnaPi0
virtual ~AliAnaPi0()
Destructor.
Definition: AliAnaPi0.cxx:127
AliAnaPi0::fhPrimEtaProdVertex
TH2F * fhPrimEtaProdVertex
! Spectrum of primary eta vs production vertex
Definition: AliAnaPi0.h:418
AliAnaPi0::fSameSM
Bool_t fSameSM
Select only pairs in same SM;.
Definition: AliAnaPi0.h:181
AliAnaPi0::fCheckAccInSector
Bool_t fCheckAccInSector
Check that the decay pi0 falls in the same SM or sector.
Definition: AliAnaPi0.h:192
AliAnaPi0::fhEventBin
TH1I * fhEventBin
! Number of real pairs in a particular bin (cen,vz,rp)
Definition: AliAnaPi0.h:317
AliAnaPi0::fhPrimPi0Pt
TH1F * fhPrimPi0Pt
! Spectrum of Primary
Definition: AliAnaPi0.h:335
AliAnaPi0::fhRe2
TH2F ** fhRe2
[GetNCentrBin()*fNPIDBits*fNAsymCuts]
Definition: AliAnaPi0.h:252
AliAnaPi0::fhPrimEtaAccE
TH1F * fhPrimEtaAccE
! Spectrum of primary with accepted daughters
Definition: AliAnaPi0.h:357
AliAnaPi0::fhMCEtaMassPtRec
TH2F ** fhMCEtaMassPtRec
[fNPtCuts*fNAsymCuts*fNCellNCuts]
Definition: AliAnaPi0.h:402
AliAnaPi0::fUseAngleCut
Bool_t fUseAngleCut
Select pairs depending on their opening angle.
Definition: AliAnaPi0.h:162
AliAnaPi0::fhRePtAsymEta
TH2F * fhRePtAsymEta
! REAL two-photon pt vs asymmetry, close to eta mass
Definition: AliAnaPi0.h:313
AliAnaCaloTrackCorrBaseClass::GetNRPBin
virtual Int_t GetNRPBin() const
Number of bins in vertex.
Definition: AliAnaCaloTrackCorrBaseClass.h:231
AliAnaPi0::fhEventPlaneResolution
TH2F * fhEventPlaneResolution
! Histogram with Event plane resolution vs centrality
Definition: AliAnaPi0.h:322
AliAnaPi0::fhMixedOpeningAngle
TH2F * fhMixedOpeningAngle
! Opening angle of pair versus pair energy
Definition: AliAnaPi0.h:328
AliAnaPi0::FillMCVersusRecDataHistograms
void FillMCVersusRecDataHistograms(Int_t index1, Int_t index2, Float_t pt1, Float_t pt2, Int_t ncells1, Int_t ncells2, Double_t mass, Double_t pt, Double_t asym, Double_t deta, Double_t dphi)
Definition: AliAnaPi0.cxx:1771
etamax
const Double_t etamax
Definition: AddTaskCFDStar.C:41
AliAnaPi0::fhMCPi0PtTruePtRec
TH2F ** fhMCPi0PtTruePtRec
[fNPtCuts*fNAsymCuts*fNCellNCuts]
Definition: AliAnaPi0.h:399
AliAnaPi0::fhPrimEtaYeta
TH2F * fhPrimEtaYeta
! PseudoRapidity distribution of primary particles vs pT
Definition: AliAnaPi0.h:361
AliAnaCaloTrackCorrBaseClass::GetMCAnalysisUtils
virtual AliMCAnalysisUtils * GetMCAnalysisUtils()
Definition: AliAnaCaloTrackCorrBaseClass.h:330
AliHistogramRanges::GetHistoAsymmetryMin
Float_t GetHistoAsymmetryMin() const
Definition: AliHistogramRanges.h:107
AliAnaPi0::fhCosThStarPrimPi0
TH2F * fhCosThStarPrimPi0
! cos(theta*) plots vs E for primary pi0, same as asymmetry ...
Definition: AliAnaPi0.h:426
AliAnaCaloTrackCorrBaseClass::Print
virtual void Print(const Option_t *) const
Print some relevant parameters set for the analysis.
Definition: AliAnaCaloTrackCorrBaseClass.cxx:611
AliAnaPi0::fhRePIDBits
TH2F ** fhRePIDBits
[fNPtCuts*fNAsymCuts*fNCellNCutsfNModules]
Definition: AliAnaPi0.h:302
AliAnaPi0::fhMCPi0ProdVertex
TH2F * fhMCPi0ProdVertex
! Spectrum of selected pi0 vs production vertex
Definition: AliAnaPi0.h:415
AliAnaPi0::fhPrimPi0PtEventPlane
TH2F * fhPrimPi0PtEventPlane
! primary pi0 reconstructed event plane vs pT
Definition: AliAnaPi0.h:349
AliAnaPi0::fhRePtNCellAsymCutsSM
TH2F ** fhRePtNCellAsymCutsSM[12]
[fNPtCuts*fNAsymCuts*fNCellNCuts]
Definition: AliAnaPi0.h:299
AliAnaPi0::fEventsList
TList ** fEventsList
Containers for photons in stored events.
Definition: AliAnaPi0.h:158
AliAnaPi0::fhMCEtaProdVertex
TH2F * fhMCEtaProdVertex
! Spectrum of selected eta vs production vertex
Definition: AliAnaPi0.h:416
AliAnaPi0::fhPrimEtaOpeningAngleAsym
TH2F * fhPrimEtaOpeningAngleAsym
! Opening angle of pair versus pair E asymmetry, eta primaries
Definition: AliAnaPi0.h:367
AliHistogramRanges::GetHistoPhiMax
Float_t GetHistoPhiMax() const
Definition: AliHistogramRanges.h:62
AliAnaCaloTrackCorrBaseClass::GetReader
virtual AliCaloTrackReader * GetReader() const
Definition: AliAnaCaloTrackCorrBaseClass.h:334
AliAnaPi0::fhMiInvPt3
TH2F ** fhMiInvPt3
[GetNCentrBin()*fNPIDBits*fNAsymCuts]
Definition: AliAnaPi0.h:288
AliAnaPi0::AliAnaPi0
AliAnaPi0()
Default Constructor. Initialized parameters with default values.
Definition: AliAnaPi0.cxx:53
AliHistogramRanges::GetHistoEtaBins
Int_t GetHistoEtaBins() const
Definition: AliHistogramRanges.h:88
AliAnaPi0::fhMi3
TH2F ** fhMi3
[GetNCentrBin()*fNPIDBits*fNAsymCuts]
Definition: AliAnaPi0.h:264
AliAnaPi0::fhArmPrimEta
TH2F * fhArmPrimEta[4]
! Armenteros plots for primary eta in 6 energy bins
Definition: AliAnaPi0.h:425
AliAnaPi0::fFillSSCombinations
Bool_t fFillSSCombinations
Do invariant mass for different combination of shower shape clusters.
Definition: AliAnaPi0.h:185
AliAnaPi0::fhPrimEtaPtOrigin
TH2F * fhPrimEtaPtOrigin
! Spectrum of generated eta vs mother
Definition: AliAnaPi0.h:377
AliAnaPi0::fhMiInvPt2
TH2F ** fhMiInvPt2
[GetNCentrBin()*fNPIDBits*fNAsymCuts]
Definition: AliAnaPi0.h:280
AliMCAnalysisUtils::CheckCommonAncestor
Int_t CheckCommonAncestor(Int_t index1, Int_t index2, const AliCaloTrackReader *reader, Int_t &ancPDG, Int_t &ancStatus, TLorentzVector &momentum, TVector3 &prodVertex)
Definition: AliMCAnalysisUtils.cxx:65
AliMCAnalysisUtils::CheckTagBit
Bool_t CheckTagBit(Int_t tag, UInt_t test) const
Definition: AliMCAnalysisUtils.h:104
AliAnaPi0::fhMCOrgDeltaPhi
TH2F * fhMCOrgDeltaPhi[13]
! Delta Phi vs pt of real pairs, check common origin of pair
Definition: AliAnaPi0.h:388
nptbins
Int_t nptbins
Definition: CalculateAveragePt.C:56
AliAnaPi0::fhMiMod
TH2F ** fhMiMod
[fNModules]
Definition: AliAnaPi0.h:226
AliAnaCaloTrackCorrBaseClass::GetMixedEvent
virtual AliMixedEvent * GetMixedEvent() const
Definition: AliAnaCaloTrackCorrBaseClass.h:265
AliAnaPi0::fhReConv2
TH2F * fhReConv2
! REAL two-photon invariant mass distribution both pair photons recombined from 2 clusters with small...
Definition: AliAnaPi0.h:241
phimin
const Double_t phimin
Definition: AddTaskCFDStar.C:32
AliAnaPi0::fhPrimEtaPhi
TH2F * fhPrimEtaPhi
! Azimutal distribution of primary particles vs pT
Definition: AliAnaPi0.h:364
AliAnaPi0::fhMCPi0MassPtTrue
TH2F ** fhMCPi0MassPtTrue
[fNPtCuts*fNAsymCuts*fNCellNCuts]
Definition: AliAnaPi0.h:396