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