AliPhysics  ec707b8 (ec707b8)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator 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  if(primStack->Energy() == TMath::Abs(primStack->Pz())) continue ; //Protection against floating point exception
1486 
1487  //printf("i %d, %s %d %s %d \n",i, stack->Particle(i)->GetName(), stack->Particle(i)->GetPdgCode(),
1488  // prim->GetName(), prim->GetPdgCode());
1489 
1490  //Photon kinematics
1491  primStack->Momentum(fPi0Mom);
1492 
1493  mesonY = 0.5*TMath::Log((primStack->Energy()+primStack->Pz())/(primStack->Energy()-primStack->Pz())) ;
1494  }
1495  else
1496  {
1497  primAOD = (AliAODMCParticle *) mcparticles->At(i);
1498  if(!primAOD)
1499  {
1500  AliWarning("AOD primaries pointer not available!!");
1501  continue;
1502  }
1503 
1504  // If too small skip
1505  if( primAOD->E() < 0.4 ) continue;
1506 
1507  pdg = primAOD->GetPdgCode();
1508  nDaught = primAOD->GetNDaughters();
1509  iphot1 = primAOD->GetFirstDaughter() ;
1510  iphot2 = primAOD->GetLastDaughter() ;
1511 
1512  if(primAOD->E() == TMath::Abs(primAOD->Pz())) continue ; //Protection against floating point exception
1513 
1514  //Photon kinematics
1515  fPi0Mom.SetPxPyPzE(primAOD->Px(),primAOD->Py(),primAOD->Pz(),primAOD->E());
1516 
1517  mesonY = 0.5*TMath::Log((primAOD->E()+primAOD->Pz())/(primAOD->E()-primAOD->Pz())) ;
1518  }
1519 
1520  // Select only pi0 or eta
1521  if( pdg != 111 && pdg != 221) continue ;
1522 
1523  mesonPt = fPi0Mom.Pt () ;
1524  mesonE = fPi0Mom.E () ;
1525  mesonYeta= fPi0Mom.Eta() ;
1526  mesonPhi = fPi0Mom.Phi() ;
1527  if( mesonPhi < 0 ) mesonPhi+=TMath::TwoPi();
1528  mesonPhi *= TMath::RadToDeg();
1529 
1530  if(pdg == 111)
1531  {
1532  if(TMath::Abs(mesonY) < 1.0)
1533  {
1534  fhPrimPi0E ->Fill(mesonE , GetEventWeight()) ;
1535  fhPrimPi0Pt ->Fill(mesonPt, GetEventWeight()) ;
1536  fhPrimPi0Phi->Fill(mesonPt, mesonPhi, GetEventWeight()) ;
1537 
1538  fhPrimPi0YetaYcut->Fill(mesonPt, mesonYeta, GetEventWeight()) ;
1539 
1540  if( IsHighMultiplicityAnalysisOn() )
1541  {
1542  fhPrimPi0PtCentrality->Fill(mesonPt, cen, GetEventWeight()) ;
1543  fhPrimPi0PtEventPlane->Fill(mesonPt, ep , GetEventWeight()) ;
1544  }
1545  }
1546 
1547  fhPrimPi0Y ->Fill(mesonPt, mesonY , GetEventWeight()) ;
1548  fhPrimPi0Yeta->Fill(mesonPt, mesonYeta, GetEventWeight()) ;
1549  }
1550  else if(pdg == 221)
1551  {
1552  if(TMath::Abs(mesonY) < 1.0)
1553  {
1554  fhPrimEtaE ->Fill(mesonE , GetEventWeight()) ;
1555  fhPrimEtaPt ->Fill(mesonPt, GetEventWeight()) ;
1556  fhPrimEtaPhi->Fill(mesonPt, mesonPhi, GetEventWeight()) ;
1557 
1558  fhPrimEtaYetaYcut->Fill(mesonPt, mesonYeta, GetEventWeight()) ;
1559 
1560  if( IsHighMultiplicityAnalysisOn() )
1561  {
1562  fhPrimEtaPtCentrality->Fill(mesonPt, cen, GetEventWeight()) ;
1563  fhPrimEtaPtEventPlane->Fill(mesonPt, ep , GetEventWeight()) ;
1564  }
1565  }
1566 
1567  fhPrimEtaY ->Fill(mesonPt, mesonY , GetEventWeight()) ;
1568  fhPrimEtaYeta->Fill(mesonPt, mesonYeta, GetEventWeight()) ;
1569  }
1570 
1571  //Origin of meson
1572  if(fFillOriginHisto && TMath::Abs(mesonY) < 0.7)
1573  {
1574  Int_t momindex = -1;
1575  Int_t mompdg = -1;
1576  Int_t momstatus = -1;
1577  Int_t status = -1;
1578  Float_t momR = 0;
1579  if(GetReader()->ReadStack()) momindex = primStack->GetFirstMother();
1580  if(GetReader()->ReadAODMCParticles()) momindex = primAOD ->GetMother();
1581 
1582  if(momindex >= 0 && momindex < nprim)
1583  {
1584  if(GetReader()->ReadStack())
1585  {
1586  status = primStack->GetStatusCode();
1587  TParticle* mother = stack->Particle(momindex);
1588  mompdg = TMath::Abs(mother->GetPdgCode());
1589  momstatus = mother->GetStatusCode();
1590  momR = mother->R();
1591  }
1592 
1593  if(GetReader()->ReadAODMCParticles())
1594  {
1595  status = primAOD->GetStatus();
1596  AliAODMCParticle* mother = (AliAODMCParticle*) mcparticles->At(momindex);
1597  mompdg = TMath::Abs(mother->GetPdgCode());
1598  momstatus = mother->GetStatus();
1599  momR = TMath::Sqrt(mother->Xv()*mother->Xv()+mother->Yv()*mother->Yv());
1600  }
1601 
1602  if(pdg == 111)
1603  {
1604  fhPrimPi0ProdVertex->Fill(mesonPt, momR , GetEventWeight());
1605  fhPrimPi0PtStatus ->Fill(mesonPt, status, GetEventWeight());
1606 
1607 
1608  if (momstatus == 21) fhPrimPi0PtOrigin->Fill(mesonPt, 0.5, GetEventWeight());//parton
1609  else if(mompdg < 22 ) fhPrimPi0PtOrigin->Fill(mesonPt, 1.5, GetEventWeight());//quark
1610  else if(mompdg > 2100 && mompdg < 2210)
1611  fhPrimPi0PtOrigin->Fill(mesonPt, 2.5, GetEventWeight());// resonances
1612  else if(mompdg == 221) fhPrimPi0PtOrigin->Fill(mesonPt, 8.5, GetEventWeight());//eta
1613  else if(mompdg == 331) fhPrimPi0PtOrigin->Fill(mesonPt, 9.5, GetEventWeight());//eta prime
1614  else if(mompdg == 213) fhPrimPi0PtOrigin->Fill(mesonPt, 4.5, GetEventWeight());//rho
1615  else if(mompdg == 223) fhPrimPi0PtOrigin->Fill(mesonPt, 5.5, GetEventWeight());//omega
1616  else if(mompdg >= 310 && mompdg <= 323)
1617  fhPrimPi0PtOrigin->Fill(mesonPt, 6.5, GetEventWeight());//k0S, k+-,k*
1618  else if(mompdg == 130) fhPrimPi0PtOrigin->Fill(mesonPt, 6.5, GetEventWeight());//k0L
1619  else if(momstatus == 11 || momstatus == 12 )
1620  fhPrimPi0PtOrigin->Fill(mesonPt, 3.5, GetEventWeight());//resonances
1621  else fhPrimPi0PtOrigin->Fill(mesonPt, 7.5, GetEventWeight());//other?
1622 
1623  //printf("Prim Pi0: index %d, pt %2.2f Prod vertex %3.3f, origin pdg %d, origin status %d, origin UI %d\n",
1624  // momindex, mesonPt,mother->R(),mompdg,momstatus,mother->GetUniqueID());
1625 
1626  if(status!=11)
1627  {
1628  if (momstatus == 21) fhPrimNotResonancePi0PtOrigin->Fill(mesonPt, 0.5, GetEventWeight());//parton
1629  else if(mompdg < 22 ) fhPrimNotResonancePi0PtOrigin->Fill(mesonPt, 1.5, GetEventWeight());//quark
1630  else if(mompdg > 2100 && mompdg < 2210)
1631  fhPrimNotResonancePi0PtOrigin->Fill(mesonPt, 2.5, GetEventWeight());// resonances
1632  else if(mompdg == 221) fhPrimNotResonancePi0PtOrigin->Fill(mesonPt, 8.5, GetEventWeight());//eta
1633  else if(mompdg == 331) fhPrimNotResonancePi0PtOrigin->Fill(mesonPt, 9.5, GetEventWeight());//eta prime
1634  else if(mompdg == 213) fhPrimNotResonancePi0PtOrigin->Fill(mesonPt, 4.5, GetEventWeight());//rho
1635  else if(mompdg == 223) fhPrimNotResonancePi0PtOrigin->Fill(mesonPt, 5.5, GetEventWeight());//omega
1636  else if(mompdg >= 310 && mompdg <= 323)
1637  fhPrimNotResonancePi0PtOrigin->Fill(mesonPt, 6.5, GetEventWeight());//k0S, k+-,k*
1638  else if(mompdg == 130) fhPrimNotResonancePi0PtOrigin->Fill(mesonPt, 6.5, GetEventWeight());//k0L
1639  else if(momstatus == 11 || momstatus == 12 )
1640  fhPrimNotResonancePi0PtOrigin->Fill(mesonPt, 3.5, GetEventWeight());//resonances
1641  else fhPrimNotResonancePi0PtOrigin->Fill(mesonPt, 7.5, GetEventWeight());//other?
1642  }
1643 
1644  }//pi0
1645  else
1646  {
1647  if (momstatus == 21 ) fhPrimEtaPtOrigin->Fill(mesonPt, 0.5, GetEventWeight());//parton
1648  else if(mompdg < 22 ) fhPrimEtaPtOrigin->Fill(mesonPt, 1.5, GetEventWeight());//quark
1649  else if(mompdg > 2100 && mompdg < 2210)
1650  fhPrimEtaPtOrigin->Fill(mesonPt, 2.5, GetEventWeight());//qq resonances
1651  else if(mompdg == 331) fhPrimEtaPtOrigin->Fill(mesonPt, 5.5, GetEventWeight());//eta prime
1652  else if(momstatus == 11 || momstatus == 12 )
1653  fhPrimEtaPtOrigin->Fill(mesonPt, 3.5, GetEventWeight());//resonances
1654  else fhPrimEtaPtOrigin->Fill(mesonPt, 4.5, GetEventWeight());//stable, conversions?
1655  //printf("Other Meson pdg %d, Mother %s, pdg %d, status %d\n",pdg, TDatabasePDG::Instance()->GetParticle(mompdg)->GetName(),mompdg, momstatus );
1656 
1657  fhPrimEtaProdVertex->Fill(mesonPt, momR, GetEventWeight());
1658  }
1659  } // pi0 has mother
1660  }
1661 
1662  //Check if both photons hit Calorimeter
1663  if(nDaught != 2 ) continue; //Only interested in 2 gamma decay
1664 
1665  if(iphot1 < 0 || iphot1 >= nprim || iphot2 < 0 || iphot2 >= nprim) continue ;
1666 
1667  Int_t pdg1 = 0;
1668  Int_t pdg2 = 0;
1669  Bool_t inacceptance1 = kTRUE;
1670  Bool_t inacceptance2 = kTRUE;
1671 
1672  if(GetReader()->ReadStack())
1673  {
1674  TParticle * phot1 = stack->Particle(iphot1) ;
1675  TParticle * phot2 = stack->Particle(iphot2) ;
1676 
1677  if(!phot1 || !phot2) continue ;
1678 
1679  pdg1 = phot1->GetPdgCode();
1680  pdg2 = phot2->GetPdgCode();
1681 
1682  phot1->Momentum(fPhotonMom1);
1683  phot2->Momentum(fPhotonMom2);
1684 
1685  // Check if photons hit the Calorimeter acceptance
1686  if(IsRealCaloAcceptanceOn())
1687  {
1688  if( !GetCaloUtils()->IsMCParticleInCalorimeterAcceptance( GetCalorimeter(), phot1 )) inacceptance1 = kFALSE ;
1689  if( !GetCaloUtils()->IsMCParticleInCalorimeterAcceptance( GetCalorimeter(), phot2 )) inacceptance2 = kFALSE ;
1690  }
1691  }
1692 
1693  if(GetReader()->ReadAODMCParticles())
1694  {
1695  AliAODMCParticle * phot1 = (AliAODMCParticle *) mcparticles->At(iphot1) ;
1696  AliAODMCParticle * phot2 = (AliAODMCParticle *) mcparticles->At(iphot2) ;
1697 
1698  if(!phot1 || !phot2) continue ;
1699 
1700  pdg1 = phot1->GetPdgCode();
1701  pdg2 = phot2->GetPdgCode();
1702 
1703  fPhotonMom1.SetPxPyPzE(phot1->Px(),phot1->Py(),phot1->Pz(),phot1->E());
1704  fPhotonMom2.SetPxPyPzE(phot2->Px(),phot2->Py(),phot2->Pz(),phot2->E());
1705 
1706  // Check if photons hit the Calorimeter acceptance
1707  if(IsRealCaloAcceptanceOn())
1708  {
1709  if( !GetCaloUtils()->IsMCParticleInCalorimeterAcceptance( GetCalorimeter(), phot1 )) inacceptance1 = kFALSE ;
1710  if( !GetCaloUtils()->IsMCParticleInCalorimeterAcceptance( GetCalorimeter(), phot2 )) inacceptance2 = kFALSE ;
1711  }
1712  }
1713 
1714  if( pdg1 != 22 || pdg2 !=22) continue ;
1715 
1716  // Check if photons hit desired acceptance in the fidutial borders fixed in the analysis
1717  if(IsFiducialCutOn())
1718  {
1719  if( inacceptance1 && !GetFiducialCut()->IsInFiducialCut(fPhotonMom1.Eta(), fPhotonMom1.Phi(), GetCalorimeter()) ) inacceptance1 = kFALSE ;
1720  if( inacceptance2 && !GetFiducialCut()->IsInFiducialCut(fPhotonMom2.Eta(), fPhotonMom2.Phi(), GetCalorimeter()) ) inacceptance2 = kFALSE ;
1721  }
1722 
1723  if(fFillArmenterosThetaStar) FillArmenterosThetaStar(pdg);
1724 
1725  if(GetCalorimeter()==kEMCAL && inacceptance1 && inacceptance2 && fCheckAccInSector)
1726  {
1727  Int_t absID1=0;
1728  Int_t absID2=0;
1729 
1730  Float_t photonPhi1 = fPhotonMom1.Phi();
1731  Float_t photonPhi2 = fPhotonMom2.Phi();
1732 
1733  if(photonPhi1 < 0) photonPhi1+=TMath::TwoPi();
1734  if(photonPhi2 < 0) photonPhi2+=TMath::TwoPi();
1735 
1736  GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(fPhotonMom1.Eta(),photonPhi1,absID1);
1737  GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(fPhotonMom2.Eta(),photonPhi2,absID2);
1738 
1739  Int_t sm1 = GetEMCALGeometry()->GetSuperModuleNumber(absID1);
1740  Int_t sm2 = GetEMCALGeometry()->GetSuperModuleNumber(absID2);
1741 
1742  Int_t j=0;
1743  Bool_t sameSector = kFALSE;
1744  for(Int_t isector = 0; isector < fNModules/2; isector++)
1745  {
1746  j=2*isector;
1747  if((sm1==j && sm2==j+1) || (sm1==j+1 && sm2==j)) sameSector = kTRUE;
1748  }
1749 
1750  if(sm1!=sm2 && !sameSector)
1751  {
1752  inacceptance1 = kFALSE;
1753  inacceptance2 = kFALSE;
1754  }
1755  //if(sm1!=sm2)printf("sm1 %d, sm2 %d, same sector %d, in acceptance %d\n",sm1,sm2,sameSector,inacceptance);
1756  // if(GetEMCALGeometry()->Impact(phot1) && GetEMCALGeometry()->Impact(phot2))
1757  // inacceptance = kTRUE;
1758  }
1759 
1760  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",
1761  GetCalorimeterString().Data(),
1762  mesonPt,mesonYeta,mesonPhi,
1763  fPhotonMom1.Pt(),fPhotonMom1.Eta(),fPhotonMom1.Phi()*TMath::RadToDeg(),
1764  fPhotonMom2.Pt(),fPhotonMom2.Eta(),fPhotonMom2.Phi()*TMath::RadToDeg(),
1765  inacceptance1, inacceptance2));
1766 
1767  if(inacceptance1 && inacceptance2)
1768  {
1769  Float_t asym = TMath::Abs((fPhotonMom1.E()-fPhotonMom2.E()) / (fPhotonMom1.E()+fPhotonMom2.E()));
1770  Double_t angle = fPhotonMom1.Angle(fPhotonMom2.Vect());
1771 
1772  Bool_t cutAngle = kFALSE;
1773  if(fUseAngleCut && (angle < fAngleCut || angle > fAngleMaxCut)) cutAngle = kTRUE;
1774 
1775  AliDebug(2,Form("\t ACCEPTED pdg %d: pt %2.2f, phi %2.2f, eta %2.2f",pdg,mesonPt,mesonPhi,mesonYeta));
1776 
1777  if(pdg==111)
1778  {
1779  fhPrimPi0AccE ->Fill(mesonE , GetEventWeight()) ;
1780  fhPrimPi0AccPt ->Fill(mesonPt, GetEventWeight()) ;
1781  fhPrimPi0AccPhi ->Fill(mesonPt, mesonPhi , GetEventWeight()) ;
1782  fhPrimPi0AccY ->Fill(mesonPt, mesonY , GetEventWeight()) ;
1783  fhPrimPi0AccYeta->Fill(mesonPt, mesonYeta, GetEventWeight()) ;
1784 
1785  if( IsHighMultiplicityAnalysisOn() )
1786  {
1787  fhPrimPi0AccPtCentrality->Fill(mesonPt, cen, GetEventWeight()) ;
1788  fhPrimPi0AccPtEventPlane->Fill(mesonPt, ep , GetEventWeight()) ;
1789  }
1790 
1791  if(fFillAngleHisto)
1792  {
1793  fhPrimPi0OpeningAngle ->Fill(mesonPt, angle, GetEventWeight());
1794  if(mesonPt > 5)fhPrimPi0OpeningAngleAsym->Fill(asym, angle, GetEventWeight());
1795  fhPrimPi0CosOpeningAngle ->Fill(mesonPt, TMath::Cos(angle), GetEventWeight());
1796  }
1797 
1798  if(fPhotonMom1.Pt() > GetMinPt() && fPhotonMom2.Pt() > GetMinPt() && !cutAngle)
1799  {
1800  fhPrimPi0AccPtPhotonCuts->Fill(mesonPt, GetEventWeight()) ;
1801 
1802  if(fFillAngleHisto)
1803  fhPrimPi0OpeningAnglePhotonCuts->Fill(mesonPt, angle, GetEventWeight());
1804  }
1805  }
1806  else if(pdg==221)
1807  {
1808  fhPrimEtaAccE ->Fill(mesonE , GetEventWeight()) ;
1809  fhPrimEtaAccPt ->Fill(mesonPt, GetEventWeight()) ;
1810  fhPrimEtaAccPhi ->Fill(mesonPt, mesonPhi , GetEventWeight()) ;
1811  fhPrimEtaAccY ->Fill(mesonPt, mesonY , GetEventWeight()) ;
1812  fhPrimEtaAccYeta->Fill(mesonPt, mesonYeta, GetEventWeight()) ;
1813 
1814  if( IsHighMultiplicityAnalysisOn() )
1815  {
1816  fhPrimEtaAccPtCentrality->Fill(mesonPt, cen, GetEventWeight()) ;
1817  fhPrimEtaAccPtEventPlane->Fill(mesonPt, ep , GetEventWeight()) ;
1818  }
1819 
1820  if(fFillAngleHisto)
1821  {
1822  fhPrimEtaOpeningAngle ->Fill(mesonPt, angle, GetEventWeight());
1823  if(mesonPt > 5)fhPrimEtaOpeningAngleAsym->Fill(asym, angle, GetEventWeight());
1824  fhPrimEtaCosOpeningAngle ->Fill(mesonPt, TMath::Cos(angle), GetEventWeight());
1825 
1826  if(fPhotonMom1.Pt() > GetMinPt() && fPhotonMom2.Pt() > GetMinPt() && !cutAngle)
1827  {
1828  if(fUseAngleCut && angle > fAngleCut && angle < fAngleMaxCut)
1829  fhPrimEtaAccPtPhotonCuts->Fill(mesonPt, GetEventWeight()) ;
1830 
1831  if(fFillAngleHisto)
1832  fhPrimEtaOpeningAnglePhotonCuts->Fill(mesonPt, angle, GetEventWeight());
1833  }
1834  }
1835  }
1836  } // Accepted
1837  } // loop on primaries
1838 }
1839 
1840 //________________________________________________
1842 //________________________________________________
1843 void AliAnaPi0::FillArmenterosThetaStar(Int_t pdg)
1844 {
1845  // Get pTArm and AlphaArm
1846  Float_t momentumSquaredMother = fPi0Mom.P()*fPi0Mom.P();
1847  Float_t momentumDaughter1AlongMother = 0.;
1848  Float_t momentumDaughter2AlongMother = 0.;
1849 
1850  if (momentumSquaredMother > 0.)
1851  {
1852  momentumDaughter1AlongMother = (fPhotonMom1.Px()*fPi0Mom.Px() + fPhotonMom1.Py()*fPi0Mom.Py()+ fPhotonMom1.Pz()*fPi0Mom.Pz()) / sqrt(momentumSquaredMother);
1853  momentumDaughter2AlongMother = (fPhotonMom2.Px()*fPi0Mom.Px() + fPhotonMom2.Py()*fPi0Mom.Py()+ fPhotonMom2.Pz()*fPi0Mom.Pz()) / sqrt(momentumSquaredMother);
1854  }
1855 
1856  Float_t momentumSquaredDaughter1 = fPhotonMom1.P()*fPhotonMom1.P();
1857  Float_t ptArmSquared = momentumSquaredDaughter1 - momentumDaughter1AlongMother*momentumDaughter1AlongMother;
1858 
1859  Float_t pTArm = 0.;
1860  if (ptArmSquared > 0.)
1861  pTArm = sqrt(ptArmSquared);
1862 
1863  Float_t alphaArm = 0.;
1864  if(momentumDaughter1AlongMother + momentumDaughter2AlongMother > 0)
1865  alphaArm = (momentumDaughter1AlongMother -momentumDaughter2AlongMother) / (momentumDaughter1AlongMother + momentumDaughter2AlongMother);
1866 
1867  fPhotonMom1Boost = fPhotonMom1;
1868  fPhotonMom1Boost.Boost(-fPi0Mom.BoostVector());
1869  Float_t cosThStar=TMath::Cos(fPhotonMom1Boost.Vect().Angle(fPi0Mom.Vect()));
1870 
1871  Float_t en = fPi0Mom.Energy();
1872  Int_t ebin = -1;
1873  if(en > 8 && en <= 12) ebin = 0;
1874  if(en > 12 && en <= 16) ebin = 1;
1875  if(en > 16 && en <= 20) ebin = 2;
1876  if(en > 20) ebin = 3;
1877  if(ebin < 0 || ebin > 3) return ;
1878 
1879  if(pdg==111)
1880  {
1881  fhCosThStarPrimPi0->Fill(en, cosThStar, GetEventWeight());
1882  fhArmPrimPi0[ebin]->Fill(alphaArm, pTArm , GetEventWeight());
1883  }
1884  else
1885  {
1886  fhCosThStarPrimEta->Fill(en, cosThStar, GetEventWeight());
1887  fhArmPrimEta[ebin]->Fill(alphaArm,pTArm , GetEventWeight());
1888  }
1889 
1890  // if(GetDebug() > 2 )
1891  // {
1892  // Float_t asym = 0;
1893  // if(fPhotonMom1.E()+fPhotonMom2.E() > 0) asym = TMath::Abs(fPhotonMom1.E()-fPhotonMom2.E())/(fPhotonMom1.E()+fPhotonMom2.E());
1894  //
1895  // printf("AliAnaPi0::FillArmenterosThetaStar() - E %f, alphaArm %f, pTArm %f, cos(theta*) %f, asymmetry %1.3f\n",
1896  // en,alphaArm,pTArm,cosThStar,asym);
1897  // }
1898 }
1899 
1900 //_______________________________________________________________________________________
1905 //_______________________________________________________________________________________
1906 void AliAnaPi0::FillMCVersusRecDataHistograms(Int_t index1, Int_t index2,
1907  Float_t pt1, Float_t pt2,
1908  Int_t ncell1, Int_t ncell2,
1909  Double_t mass, Double_t pt, Double_t asym,
1910  Double_t deta, Double_t dphi)
1911 {
1912  Int_t ancPDG = 0;
1913  Int_t ancStatus = 0;
1914  Int_t ancLabel = GetMCAnalysisUtils()->CheckCommonAncestor(index1, index2,
1915  GetReader(), ancPDG, ancStatus,fPi0Mom, fProdVertex);
1916 
1917  Int_t momindex = -1;
1918  Int_t mompdg = -1;
1919  Int_t momstatus = -1;
1920  Int_t status = -1;
1921  Float_t prodR = -1;
1922  Int_t mcIndex = -1;
1923 
1924  if(ancLabel > -1)
1925  {
1926  AliDebug(1,Form("Common ancestor label %d, pdg %d, name %s, status %d",
1927  ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus));
1928 
1929  if(ancPDG==22)
1930  {//gamma
1931  mcIndex = 0;
1932  }
1933  else if(TMath::Abs(ancPDG)==11)
1934  {//e
1935  mcIndex = 1;
1936  }
1937  else if(ancPDG==111)
1938  {//Pi0
1939  mcIndex = 2;
1940  if(fMultiCutAnaSim)
1941  {
1942  for(Int_t ipt=0; ipt<fNPtCuts; ipt++)
1943  {
1944  for(Int_t icell=0; icell<fNCellNCuts; icell++)
1945  {
1946  for(Int_t iasym=0; iasym<fNAsymCuts; iasym++)
1947  {
1948  Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
1949  if(pt1 > fPtCuts[ipt] && pt2 > fPtCuts[ipt] &&
1950  asym < fAsymCuts[iasym] &&
1951  ncell1 >= fCellNCuts[icell] && ncell2 >= fCellNCuts[icell])
1952  {
1953  fhMCPi0MassPtRec [index]->Fill(pt, mass, GetEventWeight());
1954  fhMCPi0MassPtTrue[index]->Fill(fPi0Mom.Pt(),mass, GetEventWeight());
1955  if(mass < 0.17 && mass > 0.1)
1956  fhMCPi0PtTruePtRec[index]->Fill(fPi0Mom.Pt(), pt, GetEventWeight());
1957  }//pass the different cuts
1958  }// pid bit cut loop
1959  }// icell loop
1960  }// pt cut loop
1961  }// Multi cut ana sim
1962  else
1963  {
1964  fhMCPi0MassPtTrue[0]->Fill(fPi0Mom.Pt(), mass, GetEventWeight());
1965 
1966  if(mass < 0.17 && mass > 0.1)
1967  {
1968  fhMCPi0PtTruePtRec[0]->Fill(fPi0Mom.Pt(), pt, GetEventWeight());
1969 
1970  //Int_t uniqueId = -1;
1971  if(GetReader()->ReadStack())
1972  {
1973  TParticle* ancestor = GetMCStack()->Particle(ancLabel);
1974  status = ancestor->GetStatusCode();
1975  momindex = ancestor->GetFirstMother();
1976  if(momindex < 0) return;
1977  TParticle* mother = GetMCStack()->Particle(momindex);
1978  mompdg = TMath::Abs(mother->GetPdgCode());
1979  momstatus = mother->GetStatusCode();
1980  prodR = mother->R();
1981  //uniqueId = mother->GetUniqueID();
1982  }
1983  else
1984  {
1985  TClonesArray * mcparticles = GetReader()->GetAODMCParticles();
1986  AliAODMCParticle* ancestor = (AliAODMCParticle *) mcparticles->At(ancLabel);
1987  status = ancestor->GetStatus();
1988  momindex = ancestor->GetMother();
1989  if(momindex < 0) return;
1990  AliAODMCParticle* mother = (AliAODMCParticle *) mcparticles->At(momindex);
1991  mompdg = TMath::Abs(mother->GetPdgCode());
1992  momstatus = mother->GetStatus();
1993  prodR = TMath::Sqrt(mother->Xv()*mother->Xv()+mother->Yv()*mother->Yv());
1994  //uniqueId = mother->GetUniqueID();
1995  }
1996 
1997  // printf("Reco Pi0: pt %2.2f Prod vertex %3.3f, origin pdg %d, origin status %d, origin UI %d\n",
1998  // pt,prodR,mompdg,momstatus,uniqueId);
1999 
2000  fhMCPi0ProdVertex->Fill(pt, prodR , GetEventWeight());
2001  fhMCPi0PtStatus ->Fill(pt, status, GetEventWeight());
2002 
2003  if (momstatus == 21) fhMCPi0PtOrigin->Fill(pt, 0.5, GetEventWeight());//parton
2004  else if(mompdg < 22 ) fhMCPi0PtOrigin->Fill(pt, 1.5, GetEventWeight());//quark
2005  else if(mompdg > 2100 && mompdg < 2210)
2006  fhMCPi0PtOrigin->Fill(pt, 2.5, GetEventWeight());// resonances
2007  else if(mompdg == 221) fhMCPi0PtOrigin->Fill(pt, 8.5, GetEventWeight());//eta
2008  else if(mompdg == 331) fhMCPi0PtOrigin->Fill(pt, 9.5, GetEventWeight());//eta prime
2009  else if(mompdg == 213) fhMCPi0PtOrigin->Fill(pt, 4.5, GetEventWeight());//rho
2010  else if(mompdg == 223) fhMCPi0PtOrigin->Fill(pt, 5.5, GetEventWeight());//omega
2011  else if(mompdg >= 310 && mompdg <= 323)
2012  fhMCPi0PtOrigin->Fill(pt, 6.5, GetEventWeight());//k0S, k+-,k*
2013  else if(mompdg == 130) fhMCPi0PtOrigin->Fill(pt, 6.5, GetEventWeight());//k0L
2014  else if(momstatus == 11 || momstatus == 12 )
2015  fhMCPi0PtOrigin->Fill(pt, 3.5, GetEventWeight());//resonances
2016  else fhMCPi0PtOrigin->Fill(pt, 7.5, GetEventWeight());//other?
2017 
2018  if(status!=11)
2019  {
2020  if (momstatus == 21) fhMCNotResonancePi0PtOrigin->Fill(pt, 0.5, GetEventWeight());//parton
2021  else if(mompdg < 22 ) fhMCNotResonancePi0PtOrigin->Fill(pt, 1.5, GetEventWeight());//quark
2022  else if(mompdg > 2100 && mompdg < 2210)
2023  fhMCNotResonancePi0PtOrigin->Fill(pt, 2.5, GetEventWeight());// resonances
2024  else if(mompdg == 221) fhMCNotResonancePi0PtOrigin->Fill(pt, 8.5, GetEventWeight());//eta
2025  else if(mompdg == 331) fhMCNotResonancePi0PtOrigin->Fill(pt, 9.5, GetEventWeight());//eta prime
2026  else if(mompdg == 213) fhMCNotResonancePi0PtOrigin->Fill(pt, 4.5, GetEventWeight());//rho
2027  else if(mompdg == 223) fhMCNotResonancePi0PtOrigin->Fill(pt, 5.5, GetEventWeight());//omega
2028  else if(mompdg >= 310 && mompdg <= 323)
2029  fhMCNotResonancePi0PtOrigin->Fill(pt, 6.5, GetEventWeight());//k0S, k+-,k*
2030  else if(mompdg == 130) fhMCNotResonancePi0PtOrigin->Fill(pt, 6.5, GetEventWeight());//k0L
2031  else if(momstatus == 11 || momstatus == 12 )
2032  fhMCNotResonancePi0PtOrigin->Fill(pt, 3.5, GetEventWeight());//resonances
2033  else fhMCNotResonancePi0PtOrigin->Fill(pt, 7.5, GetEventWeight());//other?
2034  }
2035  }//pi0 mass region
2036  }
2037  }
2038  else if(ancPDG==221)
2039  {//Eta
2040  mcIndex = 3;
2041  if(fMultiCutAnaSim)
2042  {
2043  for(Int_t ipt=0; ipt<fNPtCuts; ipt++)
2044  {
2045  for(Int_t icell=0; icell<fNCellNCuts; icell++)
2046  {
2047  for(Int_t iasym=0; iasym<fNAsymCuts; iasym++)
2048  {
2049  Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
2050  if(pt1 > fPtCuts[ipt] && pt2 > fPtCuts[ipt] &&
2051  asym < fAsymCuts[iasym] &&
2052  ncell1 >= fCellNCuts[icell] && ncell2 >= fCellNCuts[icell])
2053  {
2054  fhMCEtaMassPtRec [index]->Fill(pt , mass, GetEventWeight());
2055  fhMCEtaMassPtTrue[index]->Fill(fPi0Mom.Pt(), mass, GetEventWeight());
2056  if(mass < 0.65 && mass > 0.45)
2057  fhMCEtaPtTruePtRec[index]->Fill(fPi0Mom.Pt(), pt, GetEventWeight());
2058  }//pass the different cuts
2059  }// pid bit cut loop
2060  }// icell loop
2061  }// pt cut loop
2062  } //Multi cut ana sim
2063  else
2064  {
2065  fhMCEtaMassPtTrue[0]->Fill(fPi0Mom.Pt(), mass, GetEventWeight());
2066  if(mass < 0.65 && mass > 0.45) fhMCEtaPtTruePtRec[0]->Fill(fPi0Mom.Pt(), pt, GetEventWeight());
2067 
2068  if(GetReader()->ReadStack())
2069  {
2070  TParticle* ancestor = GetMCStack()->Particle(ancLabel);
2071  momindex = ancestor->GetFirstMother();
2072  if(momindex < 0) return;
2073  TParticle* mother = GetMCStack()->Particle(momindex);
2074  mompdg = TMath::Abs(mother->GetPdgCode());
2075  momstatus = mother->GetStatusCode();
2076  prodR = mother->R();
2077  }
2078  else
2079  {
2080  TClonesArray * mcparticles = GetReader()->GetAODMCParticles();
2081  AliAODMCParticle* ancestor = (AliAODMCParticle *) mcparticles->At(ancLabel);
2082  momindex = ancestor->GetMother();
2083  if(momindex < 0) return;
2084  AliAODMCParticle* mother = (AliAODMCParticle *) mcparticles->At(momindex);
2085  mompdg = TMath::Abs(mother->GetPdgCode());
2086  momstatus = mother->GetStatus();
2087  prodR = TMath::Sqrt(mother->Xv()*mother->Xv()+mother->Yv()*mother->Yv());
2088  }
2089 
2090  fhMCEtaProdVertex->Fill(pt, prodR, GetEventWeight());
2091 
2092  if (momstatus == 21 ) fhMCEtaPtOrigin->Fill(pt, 0.5, GetEventWeight());//parton
2093  else if(mompdg < 22 ) fhMCEtaPtOrigin->Fill(pt, 1.5, GetEventWeight());//quark
2094  else if(mompdg > 2100 && mompdg < 2210)
2095  fhMCEtaPtOrigin->Fill(pt, 2.5, GetEventWeight());//qq resonances
2096  else if(mompdg == 331) fhMCEtaPtOrigin->Fill(pt, 5.5, GetEventWeight());//eta prime
2097  else if(momstatus == 11 || momstatus == 12 )
2098  fhMCEtaPtOrigin->Fill(pt, 3.5, GetEventWeight());//resonances
2099  else fhMCEtaPtOrigin->Fill(pt, 4.5, GetEventWeight());//stable, conversions?
2100  //printf("Other Meson pdg %d, Mother %s, pdg %d, status %d\n",pdg, TDatabasePDG::Instance()->GetParticle(mompdg)->GetName(),mompdg, momstatus );
2101 
2102  }// eta mass region
2103  }
2104  else if(ancPDG==-2212)
2105  {//AProton
2106  mcIndex = 4;
2107  }
2108  else if(ancPDG==-2112)
2109  {//ANeutron
2110  mcIndex = 5;
2111  }
2112  else if(TMath::Abs(ancPDG)==13)
2113  {//muons
2114  mcIndex = 6;
2115  }
2116  else if (TMath::Abs(ancPDG) > 100 && ancLabel > 7)
2117  {
2118  if(ancStatus==1)
2119  {//Stable particles, converted? not decayed resonances
2120  mcIndex = 6;
2121  }
2122  else
2123  {//resonances and other decays, more hadron conversions?
2124  mcIndex = 7;
2125  }
2126  }
2127  else
2128  {//Partons, colliding protons, strings, intermediate corrections
2129  if(ancStatus==11 || ancStatus==12)
2130  {//String fragmentation
2131  mcIndex = 8;
2132  }
2133  else if (ancStatus==21){
2134  if(ancLabel < 2)
2135  {//Colliding protons
2136  mcIndex = 11;
2137  }//colliding protons
2138  else if(ancLabel < 6)
2139  {//partonic initial states interactions
2140  mcIndex = 9;
2141  }
2142  else if(ancLabel < 8)
2143  {//Final state partons radiations?
2144  mcIndex = 10;
2145  }
2146  // else {
2147  // printf("AliAnaPi0::FillMCVersusRecDataHistograms() - Check ** Common ancestor label %d, pdg %d, name %s, status %d; \n",
2148  // ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus);
2149  // }
2150  }//status 21
2151  //else {
2152  // printf("AliAnaPi0::FillMCVersusRecDataHistograms() - Check *** Common ancestor label %d, pdg %d, name %s, status %d; \n",
2153  // ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus);
2154  // }
2155  }
2156  }//ancLabel > -1
2157  else
2158  { //ancLabel <= -1
2159  //printf("Not related at all label = %d\n",ancLabel);
2160  AliDebug(1,"Common ancestor not found");
2161 
2162  mcIndex = 12;
2163  }
2164 
2165  if(mcIndex >= 0 && mcIndex < 13)
2166  {
2167  fhMCOrgMass [mcIndex]->Fill(pt, mass, GetEventWeight());
2168  fhMCOrgAsym [mcIndex]->Fill(pt, asym, GetEventWeight());
2169  fhMCOrgDeltaEta[mcIndex]->Fill(pt, deta, GetEventWeight());
2170  fhMCOrgDeltaPhi[mcIndex]->Fill(pt, dphi, GetEventWeight());
2171  }
2172 }
2173 
2174 //__________________________________________
2178 //__________________________________________
2179 void AliAnaPi0::MakeAnalysisFillHistograms()
2180 {
2181  // In case of simulated data, fill acceptance histograms
2182  if(IsDataMC())
2183  {
2184  FillAcceptanceHistograms();
2185 
2186  if(fFillOnlyMCAcceptanceHisto) return;
2187  }
2188 
2189 //if (GetReader()->GetEventNumber()%10000 == 0)
2190 // printf("--- Event %d ---\n",GetReader()->GetEventNumber());
2191 
2192  if(!GetInputAODBranch())
2193  {
2194  AliFatal(Form("No input aod photons in AOD with name branch < %s >, STOP",GetInputAODName().Data()));
2195  return;
2196  }
2197 
2198  //
2199  // Init some variables
2200  //
2201 
2202  // Analysis within the same detector:
2203  TClonesArray* secondLoopInputData = GetInputAODBranch();
2204 
2205  Int_t nPhot = GetInputAODBranch()->GetEntriesFast() ;
2206  Int_t nPhot2 = nPhot; // second loop
2207  Int_t minEntries = 2;
2208  Int_t last = 1; // last entry in first loop to be removed
2209 
2210  // Combine 2 detectors:
2211  if(fPairWithOtherDetector)
2212  {
2213  // Input from CaloTrackCorr tasks
2214  secondLoopInputData = (TClonesArray *) GetReader()->GetAODBranchList()->FindObject(fOtherDetectorInputName);
2215 
2216  // In case input from PCM or other external source
2217  if(!secondLoopInputData) secondLoopInputData = (TClonesArray *) GetReader()->GetOutputEvent()->FindObject(fOtherDetectorInputName);
2218  if(!secondLoopInputData) secondLoopInputData = (TClonesArray *) GetReader()->GetInputEvent ()->FindObject(fOtherDetectorInputName);
2219 
2220  if(!secondLoopInputData)
2221  {
2222  AliFatal(Form("Input for other detector not found in branch %s!",fOtherDetectorInputName.Data()));
2223  return ; // coverity
2224  }
2225 
2226  nPhot2 = secondLoopInputData->GetEntriesFast() ; // add one since we remove one for the normal case
2227  minEntries = 1;
2228  last = 0;
2229  }
2230 
2231  AliDebug(1,Form("Photon entries %d", nPhot));
2232 
2233  // If less than photon 2 (1) entries in the list, skip this event
2234  if ( nPhot < minEntries )
2235  {
2236  AliDebug(1,Form("nPhotons %d, cent bin %d continue to next event",nPhot, GetEventCentralityBin()));
2237 
2238  if ( GetNCentrBin() > 1 && IsHighMultiplicityAnalysisOn() )
2239  fhCentralityNoPair->Fill(GetEventCentralityBin(), GetEventWeight());
2240 
2241  return ;
2242  }
2243  else if(fPairWithOtherDetector && nPhot2 < minEntries)
2244  {
2245  AliDebug(1,Form("nPhotons %d, cent bin %d continue to next event",nPhot, GetEventCentralityBin()));
2246 
2247  if ( GetNCentrBin() > 1 && IsHighMultiplicityAnalysisOn() )
2248  fhCentralityNoPair->Fill(GetEventCentralityBin(), GetEventWeight());
2249 
2250  return ;
2251  }
2252 
2253  Int_t ncentr = GetNCentrBin();
2254  if(GetNCentrBin()==0) ncentr = 1;
2255 
2256  //Init variables
2257  Int_t module1 = -1;
2258  Int_t module2 = -1;
2259  Double_t vert[] = {0.0, 0.0, 0.0} ; //vertex
2260  Int_t evtIndex1 = 0 ;
2261  Int_t currentEvtIndex = -1;
2262  Int_t curCentrBin = GetEventCentralityBin();
2263  //Int_t curVzBin = GetEventVzBin();
2264  //Int_t curRPBin = GetEventRPBin();
2265  Int_t eventbin = GetEventMixBin();
2266 
2267  if(eventbin > GetNCentrBin()*GetNZvertBin()*GetNRPBin())
2268  {
2269  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",
2270  GetEventCentralityBin(),GetEventVzBin(), GetEventRPBin(),eventbin,GetEventCentrality(),vert[2],GetEventPlaneAngle()));
2271  return;
2272  }
2273 
2274  // Get shower shape information of clusters
2275 // TObjArray *clusters = 0;
2276 // if (GetCalorimeter()==kEMCAL) clusters = GetEMCALClusters();
2277 // else if(GetCalorimeter()==kPHOS ) clusters = GetPHOSClusters() ;
2278 
2279  //---------------------------------
2280  // First loop on photons/clusters
2281  //---------------------------------
2282  for(Int_t i1 = 0; i1 < nPhot-last; i1++)
2283  {
2284  AliAODPWG4Particle * p1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i1)) ;
2285 
2286  // Select photons within a pT range
2287  if ( p1->Pt() < GetMinPt() || p1->Pt() > GetMaxPt() ) continue ;
2288 
2289  //printf("AliAnaPi0::MakeAnalysisFillHistograms() : cluster1 id %d/%d\n",i1,nPhot-1);
2290 
2291  // get the event index in the mixed buffer where the photon comes from
2292  // in case of mixing with analysis frame, not own mixing
2293  evtIndex1 = GetEventIndex(p1, vert) ;
2294  if ( evtIndex1 == -1 )
2295  return ;
2296  if ( evtIndex1 == -2 )
2297  continue ;
2298 
2299  // Only effective in case of mixed event frame
2300  if(TMath::Abs(vert[2]) > GetZvertexCut()) continue ; //vertex cut
2301 
2302  if (evtIndex1 != currentEvtIndex)
2303  {
2304  // Fill event bin info
2305  if( DoOwnMix() ) fhEventBin->Fill(eventbin, GetEventWeight()) ;
2306 
2307  if( IsHighMultiplicityAnalysisOn() )
2308  {
2309  fhCentrality->Fill(curCentrBin, GetEventWeight());
2310 
2311  if( GetEventPlane() )
2312  fhEventPlaneResolution->Fill(curCentrBin, TMath::Cos(2.*GetEventPlane()->GetQsubRes()), GetEventWeight());
2313  }
2314 
2315  currentEvtIndex = evtIndex1 ;
2316  }
2317 
2318  //printf("AliAnaPi0::MakeAnalysisFillHistograms(): Photon 1 Evt %d Vertex : %f,%f,%f\n",evtIndex1, GetVertex(evtIndex1)[0] ,GetVertex(evtIndex1)[1],GetVertex(evtIndex1)[2]);
2319 
2320  //Get the momentum of this cluster
2321  fPhotonMom1.SetPxPyPzE(p1->Px(),p1->Py(),p1->Pz(),p1->E());
2322 
2323  //Get (Super)Module number of this cluster
2324  module1 = p1->GetSModNumber();// GetModuleNumber(p1);
2325 
2326  //------------------------------------------
2327  // Recover original cluster
2328  // Int_t iclus1 = -1 ;
2329  // AliVCluster * cluster1 = FindCluster(clusters,p1->GetCaloLabel(0),iclus1);
2330  // if(!cluster1) AliWarning("Cluster1 not found!");
2331 
2332  //---------------------------------
2333  // Second loop on photons/clusters
2334  //---------------------------------
2335  Int_t first = i1+1;
2336  if(fPairWithOtherDetector) first = 0;
2337 
2338  for(Int_t i2 = first; i2 < nPhot2; i2++)
2339  {
2340  //AliAODPWG4Particle * p2 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i2)) ;
2341  AliAODPWG4Particle * p2 = (AliAODPWG4Particle*) (secondLoopInputData->At(i2)) ;
2342 
2343  // Select photons within a pT range
2344  if ( p2->Pt() < GetMinPt() || p2->Pt() > GetMaxPt() ) continue ;
2345 
2346  //printf("AliAnaPi0::MakeAnalysisFillHistograms() : cluster2 i %d/%d\n",i2,nPhot);
2347 
2348  //In case of mixing frame, check we are not in the same event as the first cluster
2349  Int_t evtIndex2 = GetEventIndex(p2, vert) ;
2350  if ( evtIndex2 == -1 )
2351  return ;
2352  if ( evtIndex2 == -2 )
2353  continue ;
2354  if (GetMixedEvent() && (evtIndex1 == evtIndex2))
2355  continue ;
2356 
2357 // //------------------------------------------
2358 // // Recover original cluster
2359 // Int_t iclus2 = -1;
2360 // AliVCluster * cluster2 = FindCluster(clusters,p2->GetCaloLabel(0),iclus2,iclus1+1);
2361 // // start new loop from iclus1+1 to gain some time
2362 // if(!cluster2) AliWarning("Cluster2 not found!");
2363 //
2364 // // Get the TOF,l0 and ncells from the clusters
2365 // Float_t tof1 = -1;
2366 // Float_t l01 = -1;
2367 // Int_t ncell1 = 0;
2368 // if(cluster1)
2369 // {
2370 // tof1 = cluster1->GetTOF()*1e9;
2371 // l01 = cluster1->GetM02();
2372 // ncell1 = cluster1->GetNCells();
2373 // //printf("cluster1: E %2.2f (%2.2f), l0 %2.2f, tof %2.2f\n",cluster1->E(),p1->E(),l01,tof1);
2374 // }
2375 // //else printf("cluster1 not available: calo label %d / %d, cluster ID %d\n",
2376 // // p1->GetCaloLabel(0),(GetReader()->GetInputEvent())->GetNumberOfCaloClusters()-1,cluster1->GetID());
2377 //
2378 // Float_t tof2 = -1;
2379 // Float_t l02 = -1;
2380 // Int_t ncell2 = 0;
2381 // if(cluster2)
2382 // {
2383 // tof2 = cluster2->GetTOF()*1e9;
2384 // l02 = cluster2->GetM02();
2385 // ncell2 = cluster2->GetNCells();
2386 // //printf("cluster2: E %2.2f (%2.2f), l0 %2.2f, tof %2.2f\n",cluster2->E(),p2->E(),l02,tof2);
2387 // }
2388 // //else printf("cluster2 not available: calo label %d / %d, cluster ID %d\n",
2389 // // p2->GetCaloLabel(0),(GetReader()->GetInputEvent())->GetNumberOfCaloClusters()-1,cluster2->GetID());
2390 //
2391 // if(cluster1 && cluster2)
2392 // {
2393 // Double_t t12diff = tof1-tof2;
2394 // if(TMath::Abs(t12diff) > GetPairTimeCut()) continue;
2395 // }
2396 
2397  Float_t tof1 = p1->GetTime();
2398  Float_t l01 = p1->GetM02();
2399  Int_t ncell1 = p1->GetNCells();
2400  //printf("cluster1: E %2.2f, l0 %2.2f, tof %2.2f\n",p1->E(),l01,tof1);
2401 
2402  Float_t tof2 = p2->GetTime();
2403  Float_t l02 = p2->GetM02();
2404  Int_t ncell2 = p2->GetNCells();
2405  //printf("cluster2: E %2.2f, l0 %2.2f, tof %2.2f\n",p2->E(),l02,tof2);
2406 
2407  Double_t t12diff = tof1-tof2;
2408  fhEPairDiffTime->Fill((fPhotonMom1 + fPhotonMom2).Pt(), t12diff, GetEventWeight());
2409  if(TMath::Abs(t12diff) > GetPairTimeCut()) continue;
2410 
2411  //------------------------------------------
2412 
2413  //printf("AliAnaPi0::MakeAnalysisFillHistograms(): Photon 2 Evt %d Vertex : %f,%f,%f\n",evtIndex2, GetVertex(evtIndex2)[0] ,GetVertex(evtIndex2)[1],GetVertex(evtIndex2)[2]);
2414 
2415  // Get the momentum of this cluster
2416  fPhotonMom2.SetPxPyPzE(p2->Px(),p2->Py(),p2->Pz(),p2->E());
2417 
2418  // Get module number
2419  module2 = p2->GetSModNumber(); //GetModuleNumber(p2);
2420 
2421  //---------------------------------
2422  // Get pair kinematics
2423  //---------------------------------
2424  Double_t m = (fPhotonMom1 + fPhotonMom2).M() ;
2425  Double_t pt = (fPhotonMom1 + fPhotonMom2).Pt();
2426  Double_t deta = fPhotonMom1.Eta() - fPhotonMom2.Eta();
2427  Double_t dphi = fPhotonMom1.Phi() - fPhotonMom2.Phi();
2428  Double_t a = TMath::Abs(p1->E()-p2->E())/(p1->E()+p2->E()) ;
2429 
2430  AliDebug(2,Form("E: fPhotonMom1 %f, fPhotonMom2 %f; Pair: pT %f, mass %f, a %f", p1->E(), p2->E(), (fPhotonMom1 + fPhotonMom2).E(),m,a));
2431 
2432  //--------------------------------
2433  // Opening angle selection
2434  //--------------------------------
2435  // Check if opening angle is too large or too small compared to what is expected
2436  Double_t angle = fPhotonMom1.Angle(fPhotonMom2.Vect());
2437  if(fUseAngleEDepCut && !GetNeutralMesonSelection()->IsAngleInWindow((fPhotonMom1+fPhotonMom2).E(),angle+0.05))
2438  {
2439  AliDebug(2,Form("Real pair angle %f (deg) not in E %f window",RadToDeg(angle), (fPhotonMom1+fPhotonMom2).E()));
2440  continue;
2441  }
2442 
2443  if(fUseAngleCut && (angle < fAngleCut || angle > fAngleMaxCut))
2444  {
2445  AliDebug(2,Form("Real pair cut %f < angle %f < cut %f (deg)",RadToDeg(fAngleCut), RadToDeg(angle), RadToDeg(fAngleMaxCut)));
2446  continue;
2447  }
2448 
2449  //-------------------------------------------------------------------------------------------------
2450  // Fill module dependent histograms, put a cut on assymmetry on the first available cut in the array
2451  //-------------------------------------------------------------------------------------------------
2452  if(a < fAsymCuts[0] && fFillSMCombinations)
2453  {
2454  if(!fPairWithOtherDetector)
2455  {
2456  if(module1==module2 && module1 >=0 && module1<fNModules)
2457  fhReMod[module1]->Fill(pt, m, GetEventWeight()) ;
2458 
2459  if (GetCalorimeter() == kEMCAL )
2460  {
2461  // Same sector
2462  Int_t j=0;
2463  for(Int_t i = 0; i < fNModules/2; i++)
2464  {
2465  j=2*i;
2466  if((module1==j && module2==j+1) || (module1==j+1 && module2==j)) fhReSameSectorEMCALMod[i]->Fill(pt, m, GetEventWeight()) ;
2467  }
2468 
2469  // Same side
2470  for(Int_t i = 0; i < fNModules-2; i++){
2471  if((module1==i && module2==i+2) || (module1==i+2 && module2==i)) fhReSameSideEMCALMod[i]->Fill(pt, m, GetEventWeight());
2472  }
2473  } // EMCAL
2474  else
2475  { // PHOS
2476  if((module1==0 && module2==1) || (module1==1 && module2==0)) fhReDiffPHOSMod[0]->Fill(pt, m, GetEventWeight()) ;
2477  if((module1==0 && module2==2) || (module1==2 && module2==0)) fhReDiffPHOSMod[1]->Fill(pt, m, GetEventWeight()) ;
2478  if((module1==1 && module2==2) || (module1==2 && module2==1)) fhReDiffPHOSMod[2]->Fill(pt, m, GetEventWeight()) ;
2479  } // PHOS
2480  }
2481  else
2482  {
2483  Float_t phi1 = GetPhi(fPhotonMom1.Phi());
2484  Float_t phi2 = GetPhi(fPhotonMom2.Phi());
2485  Bool_t etaside = 0;
2486  if( (p1->GetDetectorTag()==kEMCAL && fPhotonMom1.Eta() < 0)
2487  || (p2->GetDetectorTag()==kEMCAL && fPhotonMom2.Eta() < 0)) etaside = 1;
2488 
2489  if ( phi1 > DegToRad(260) && phi2 > DegToRad(260) && phi1 < DegToRad(280) && phi2 < DegToRad(280)) fhReSameSectorDCALPHOSMod[0+etaside]->Fill(pt, m, GetEventWeight());
2490  else if ( phi1 > DegToRad(280) && phi2 > DegToRad(280) && phi1 < DegToRad(300) && phi2 < DegToRad(300)) fhReSameSectorDCALPHOSMod[2+etaside]->Fill(pt, m, GetEventWeight());
2491  else if ( phi1 > DegToRad(300) && phi2 > DegToRad(300) && phi1 < DegToRad(320) && phi2 < DegToRad(320)) fhReSameSectorDCALPHOSMod[4+etaside]->Fill(pt, m, GetEventWeight());
2492  else if ( (phi1 > DegToRad(260) && phi2 > DegToRad(280) && phi1 < DegToRad(280) && phi2 < DegToRad(300))
2493  || (phi1 > DegToRad(280) && phi2 > DegToRad(260) && phi1 < DegToRad(300) && phi2 < DegToRad(280))) fhReDiffSectorDCALPHOSMod[0+etaside]->Fill(pt, m, GetEventWeight());
2494  else if ( (phi1 > DegToRad(280) && phi2 > DegToRad(300) && phi1 < DegToRad(300) && phi2 < DegToRad(320))
2495  || (phi1 > DegToRad(300) && phi2 > DegToRad(280) && phi1 < DegToRad(320) && phi2 < DegToRad(300))) fhReDiffSectorDCALPHOSMod[2+etaside]->Fill(pt, m, GetEventWeight());
2496  else if ( (phi1 > DegToRad(260) && phi2 > DegToRad(300) && phi1 < DegToRad(280) && phi2 < DegToRad(320))
2497  || (phi1 > DegToRad(300) && phi2 > DegToRad(260) && phi1 < DegToRad(320) && phi2 < DegToRad(280))) fhReDiffSectorDCALPHOSMod[4+etaside]->Fill(pt, m, GetEventWeight());
2498  else fhReDiffSectorDCALPHOSMod[6+etaside]->Fill(pt, m, GetEventWeight());
2499  }
2500  } // Fill SM combinations
2501 
2502  // In case we want only pairs in same (super) module, check their origin.
2503  Bool_t ok = kTRUE;
2504 
2505  if(fSameSM)
2506  {
2507  if(!fPairWithOtherDetector)
2508  {
2509  if(module1!=module2) ok=kFALSE;
2510  }
2511  else // PHOS and DCal in same sector
2512  {
2513  Float_t phi1 = GetPhi(fPhotonMom1.Phi());
2514  Float_t phi2 = GetPhi(fPhotonMom2.Phi());
2515  ok=kFALSE;
2516  if ( phi1 > DegToRad(260) && phi2 > DegToRad(260) && phi1 < DegToRad(280) && phi2 < DegToRad(280)) ok = kTRUE;
2517  else if ( phi1 > DegToRad(280) && phi2 > DegToRad(280) && phi1 < DegToRad(300) && phi2 < DegToRad(300)) ok = kTRUE;
2518  else if ( phi1 > DegToRad(300) && phi2 > DegToRad(300) && phi1 < DegToRad(320) && phi2 < DegToRad(320)) ok = kTRUE;
2519  }
2520  } // Pair only in same SM
2521 
2522  if(!ok) continue;
2523 
2524  //
2525  // Fill histograms with selected cluster pairs
2526  //
2527 
2528  // Check if one of the clusters comes from a conversion
2529  if(fCheckConversion)
2530  {
2531  if (p1->IsTagged() && p2->IsTagged()) fhReConv2->Fill(pt, m, GetEventWeight());
2532  else if(p1->IsTagged() || p2->IsTagged()) fhReConv ->Fill(pt, m, GetEventWeight());
2533  }
2534 
2535  // Fill shower shape cut histograms
2536  if(fFillSSCombinations)
2537  {
2538  if ( l01 > 0.01 && l01 < 0.4 &&
2539  l02 > 0.01 && l02 < 0.4 ) fhReSS[0]->Fill(pt, m, GetEventWeight()); // Tight
2540  else if( l01 > 0.4 && l02 > 0.4 ) fhReSS[1]->Fill(pt, m, GetEventWeight()); // Loose
2541  else if( l01 > 0.01 && l01 < 0.4 && l02 > 0.4 ) fhReSS[2]->Fill(pt, m, GetEventWeight()); // Both
2542  else if( l02 > 0.01 && l02 < 0.4 && l01 > 0.4 ) fhReSS[2]->Fill(pt, m, GetEventWeight()); // Both
2543  }
2544 
2545  // Fill histograms for different bad channel distance, centrality, assymmetry cut and pid bit
2546  for(Int_t ipid=0; ipid<fNPIDBits; ipid++)
2547  {
2548  if((p1->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton)) && (p2->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton)))
2549  {
2550  for(Int_t iasym=0; iasym < fNAsymCuts; iasym++)
2551  {
2552  if(a < fAsymCuts[iasym])
2553  {
2554  Int_t index = ((curCentrBin*fNPIDBits)+ipid)*fNAsymCuts + iasym;
2555  //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);
2556 
2557  if(index < 0 || index >= ncentr*fNPIDBits*fNAsymCuts) continue ;
2558 
2559  fhRe1 [index]->Fill(pt, m, GetEventWeight());
2560  if(fMakeInvPtPlots)fhReInvPt1[index]->Fill(pt, m, 1./pt * GetEventWeight()) ;
2561 
2562  if(fFillBadDistHisto)
2563  {
2564  if(p1->DistToBad()>0 && p2->DistToBad()>0)
2565  {
2566  fhRe2 [index]->Fill(pt, m, GetEventWeight()) ;
2567  if(fMakeInvPtPlots)fhReInvPt2[index]->Fill(pt, m, 1./pt * GetEventWeight()) ;
2568 
2569  if(p1->DistToBad()>1 && p2->DistToBad()>1)
2570  {
2571  fhRe3 [index]->Fill(pt, m, GetEventWeight()) ;
2572  if(fMakeInvPtPlots)fhReInvPt3[index]->Fill(pt, m, 1./pt * GetEventWeight()) ;
2573  }// bad 3
2574  }// bad2
2575  }// Fill bad dist histos
2576  }//assymetry cut
2577  }// asymmetry cut loop
2578  }// bad 1
2579  }// pid bit loop
2580 
2581  // Fill histograms with opening angle
2582  if(fFillAngleHisto)
2583  {
2584  fhRealOpeningAngle ->Fill(pt, angle, GetEventWeight());
2585  fhRealCosOpeningAngle->Fill(pt, TMath::Cos(angle), GetEventWeight());
2586  }
2587 
2588  // Fill histograms with pair assymmetry
2589  if(fFillAsymmetryHisto)
2590  {
2591  fhRePtAsym->Fill(pt, a, GetEventWeight());
2592  if(m > 0.10 && m < 0.17) fhRePtAsymPi0->Fill(pt, a, GetEventWeight());
2593  if(m > 0.45 && m < 0.65) fhRePtAsymEta->Fill(pt, a, GetEventWeight());
2594  }
2595 
2596  // Check cell time content in cluster
2597  if ( fFillSecondaryCellTiming)
2598  {
2599  if ( p1->GetFiducialArea() == 0 && p2->GetFiducialArea() == 0 )
2600  fhReSecondaryCellInTimeWindow ->Fill(pt, m, GetEventWeight());
2601 
2602  else if ( p1->GetFiducialArea() != 0 && p2->GetFiducialArea() != 0 )
2603  fhReSecondaryCellOutTimeWindow->Fill(pt, m, GetEventWeight());
2604  }
2605 
2606  //---------
2607  // MC data
2608  //---------
2609  // Do some MC checks on the origin of the pair, is there any common ancestor and if there is one, who?
2610  if(IsDataMC())
2611  {
2612  if(GetMCAnalysisUtils()->CheckTagBit(p1->GetTag(),AliMCAnalysisUtils::kMCConversion) &&
2613  GetMCAnalysisUtils()->CheckTagBit(p2->GetTag(),AliMCAnalysisUtils::kMCConversion))
2614  {
2615  fhReMCFromConversion->Fill(pt, m, GetEventWeight());
2616  }
2617  else if(!GetMCAnalysisUtils()->CheckTagBit(p1->GetTag(),AliMCAnalysisUtils::kMCConversion) &&
2618  !GetMCAnalysisUtils()->CheckTagBit(p2->GetTag(),AliMCAnalysisUtils::kMCConversion))
2619  {
2620  fhReMCFromNotConversion->Fill(pt, m, GetEventWeight());
2621  }
2622  else
2623  {
2624  fhReMCFromMixConversion->Fill(pt, m, GetEventWeight());
2625  }
2626 
2627  if(fFillOriginHisto)
2628  FillMCVersusRecDataHistograms(p1->GetLabel(), p2->GetLabel(),p1->Pt(), p2->Pt(),ncell1, ncell2, m, pt, a,deta, dphi);
2629  }
2630 
2631  //-----------------------
2632  // Multi cuts analysis
2633  //-----------------------
2634  if(fMultiCutAna)
2635  {
2636  // Histograms for different PID bits selection
2637  for(Int_t ipid=0; ipid<fNPIDBits; ipid++)
2638  {
2639  if(p1->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton) &&
2640  p2->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton)) fhRePIDBits[ipid]->Fill(pt, m, GetEventWeight()) ;
2641 
2642  //printf("ipt %d, ipid%d, name %s\n",ipt, ipid, fhRePtPIDCuts[ipt*fNPIDBitsBits+ipid]->GetName());
2643  } // pid bit cut loop
2644 
2645  // Several pt,ncell and asymmetry cuts
2646  for(Int_t ipt = 0; ipt < fNPtCuts; ipt++)
2647  {
2648  for(Int_t icell = 0; icell < fNCellNCuts; icell++)
2649  {
2650  for(Int_t iasym = 0; iasym < fNAsymCuts; iasym++)
2651  {
2652  Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
2653  if(p1->E() > fPtCuts[ipt] && p2->E() > fPtCuts[ipt] &&
2654  a < fAsymCuts[iasym] &&
2655  ncell1 >= fCellNCuts[icell] && ncell2 >= fCellNCuts[icell])
2656  {
2657  fhRePtNCellAsymCuts[index]->Fill(pt, m, GetEventWeight()) ;
2658  //printf("ipt %d, icell%d, iasym %d, name %s\n",ipt, icell, iasym, fhRePtNCellAsymCuts[((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym]->GetName());
2659  if(fFillSMCombinations && module1==module2)
2660  {
2661  fhRePtNCellAsymCutsSM[module1][index]->Fill(pt, m, GetEventWeight()) ;
2662  }
2663  }
2664  }// pid bit cut loop
2665  }// icell loop
2666  }// pt cut loop
2667 
2668  if(GetHistogramRanges()->GetHistoTrackMultiplicityBins())
2669  {
2670  for(Int_t iasym = 0; iasym < fNAsymCuts; iasym++)
2671  {
2672  if(a < fAsymCuts[iasym]) fhRePtMult[iasym]->Fill(pt, GetTrackMultiplicity(), m, GetEventWeight()) ;
2673  }
2674  }
2675  }// multiple cuts analysis
2676  }// second same event particle
2677  }// first cluster
2678 
2679  //-------------------------------------------------------------
2680  // Mixing
2681  //-------------------------------------------------------------
2682  if(DoOwnMix())
2683  {
2684  // Recover events in with same characteristics as the current event
2685 
2686  // Check that the bin exists, if not (bad determination of RP, centrality or vz bin) do nothing
2687  if(eventbin < 0) return ;
2688 
2689  TList * evMixList=fEventsList[eventbin] ;
2690 
2691  if(!evMixList)
2692  {
2693  AliWarning(Form("Mix event list not available, bin %d",eventbin));
2694  return;
2695  }
2696 
2697  Int_t nMixed = evMixList->GetSize() ;
2698  for(Int_t ii=0; ii<nMixed; ii++)
2699  {
2700  TClonesArray* ev2= (TClonesArray*) (evMixList->At(ii));
2701  Int_t nPhot2=ev2->GetEntriesFast() ;
2702  Double_t m = -999;
2703  AliDebug(1,Form("Mixed event %d photon entries %d, centrality bin %d",ii, nPhot2, GetEventCentralityBin()));
2704 
2705  fhEventMixBin->Fill(eventbin, GetEventWeight()) ;
2706 
2707  //---------------------------------
2708  // First loop on photons/clusters
2709  //---------------------------------
2710  for(Int_t i1 = 0; i1 < nPhot; i1++)
2711  {
2712  AliAODPWG4Particle * p1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i1)) ;
2713 
2714  // Select photons within a pT range
2715  if ( p1->Pt() < GetMinPt() || p1->Pt() > GetMaxPt() ) continue ;
2716 
2717  // Not sure why this line is here
2718  //if(fSameSM && GetModuleNumber(p1)!=module1) continue;
2719 
2720  //Get kinematics of cluster and (super) module of this cluster
2721  fPhotonMom1.SetPxPyPzE(p1->Px(),p1->Py(),p1->Pz(),p1->E());
2722  module1 = GetModuleNumber(p1);
2723 
2724  //---------------------------------
2725  // Second loop on other mixed event photons/clusters
2726  //---------------------------------
2727  for(Int_t i2 = 0; i2 < nPhot2; i2++)
2728  {
2729  AliAODPWG4Particle * p2 = (AliAODPWG4Particle*) (ev2->At(i2)) ;
2730 
2731  // Select photons within a pT range
2732  if ( p2->Pt() < GetMinPt() || p2->Pt() > GetMaxPt() ) continue ;
2733 
2734  // Get kinematics of second cluster and calculate those of the pair
2735  fPhotonMom2.SetPxPyPzE(p2->Px(),p2->Py(),p2->Pz(),p2->E());
2736  m = (fPhotonMom1+fPhotonMom2).M() ;
2737  Double_t pt = (fPhotonMom1 + fPhotonMom2).Pt();
2738  Double_t a = TMath::Abs(p1->E()-p2->E())/(p1->E()+p2->E()) ;
2739 
2740  // Check if opening angle is too large or too small compared to what is expected
2741  Double_t angle = fPhotonMom1.Angle(fPhotonMom2.Vect());
2742  if(fUseAngleEDepCut && !GetNeutralMesonSelection()->IsAngleInWindow((fPhotonMom1+fPhotonMom2).E(),angle+0.05))
2743  {
2744  AliDebug(2,Form("Mix pair angle %f (deg) not in E %f window",RadToDeg(angle), (fPhotonMom1+fPhotonMom2).E()));
2745  continue;
2746  }
2747 
2748  if(fUseAngleCut && (angle < fAngleCut || angle > fAngleMaxCut))
2749  {
2750  AliDebug(2,Form("Mix pair cut %f < angle %f < cut %f (deg)",RadToDeg(fAngleCut),RadToDeg(angle),RadToDeg(fAngleMaxCut)));
2751  continue;
2752  }
2753 
2754  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));
2755 
2756  // In case we want only pairs in same (super) module, check their origin.
2757  module2 = GetModuleNumber(p2);
2758 
2759  //-------------------------------------------------------------------------------------------------
2760  // Fill module dependent histograms, put a cut on assymmetry on the first available cut in the array
2761  //-------------------------------------------------------------------------------------------------
2762  if(a < fAsymCuts[0] && fFillSMCombinations)
2763  {
2764  if(!fPairWithOtherDetector)
2765  {
2766  if(module1==module2 && module1 >=0 && module1<fNModules)
2767  fhMiMod[module1]->Fill(pt, m, GetEventWeight()) ;
2768 
2769  if(GetCalorimeter()==kEMCAL)
2770  {
2771  // Same sector
2772  Int_t j=0;
2773  for(Int_t i = 0; i < fNModules/2; i++)
2774  {
2775  j=2*i;
2776  if((module1==j && module2==j+1) || (module1==j+1 && module2==j)) fhMiSameSectorEMCALMod[i]->Fill(pt, m, GetEventWeight()) ;
2777  }
2778 
2779  // Same side
2780  for(Int_t i = 0; i < fNModules-2; i++)
2781  {
2782  if((module1==i && module2==i+2) || (module1==i+2 && module2==i)) fhMiSameSideEMCALMod[i]->Fill(pt, m, GetEventWeight());
2783  }
2784  } // EMCAL
2785  else
2786  { // PHOS
2787  if((module1==0 && module2==1) || (module1==1 && module2==0)) fhMiDiffPHOSMod[0]->Fill(pt, m, GetEventWeight()) ;
2788  if((module1==0 && module2==2) || (module1==2 && module2==0)) fhMiDiffPHOSMod[1]->Fill(pt, m, GetEventWeight()) ;
2789  if((module1==1 && module2==2) || (module1==2 && module2==1)) fhMiDiffPHOSMod[2]->Fill(pt, m, GetEventWeight()) ;
2790  } // PHOS
2791  }
2792  else
2793  {
2794  Float_t phi1 = GetPhi(fPhotonMom1.Phi());
2795  Float_t phi2 = GetPhi(fPhotonMom2.Phi());
2796  Bool_t etaside = 0;
2797  if( (p1->GetDetectorTag()==kEMCAL && fPhotonMom1.Eta() < 0)
2798  || (p2->GetDetectorTag()==kEMCAL && fPhotonMom2.Eta() < 0)) etaside = 1;
2799 
2800  if ( phi1 > DegToRad(260) && phi2 > DegToRad(260) && phi1 < DegToRad(280) && phi2 < DegToRad(280)) fhMiSameSectorDCALPHOSMod[0+etaside]->Fill(pt, m, GetEventWeight());
2801  else if ( phi1 > DegToRad(280) && phi2 > DegToRad(280) && phi1 < DegToRad(300) && phi2 < DegToRad(300)) fhMiSameSectorDCALPHOSMod[2+etaside]->Fill(pt, m, GetEventWeight());
2802  else if ( phi1 > DegToRad(300) && phi2 > DegToRad(300) && phi1 < DegToRad(320) && phi2 < DegToRad(320)) fhMiSameSectorDCALPHOSMod[4+etaside]->Fill(pt, m, GetEventWeight());
2803  else if ( (phi1 > DegToRad(260) && phi2 > DegToRad(280) && phi1 < DegToRad(280) && phi2 < DegToRad(300))
2804  || (phi1 > DegToRad(280) && phi2 > DegToRad(260) && phi1 < DegToRad(300) && phi2 < DegToRad(280))) fhMiDiffSectorDCALPHOSMod[0+etaside]->Fill(pt, m, GetEventWeight());
2805  else if ( (phi1 > DegToRad(280) && phi2 > DegToRad(300) && phi1 < DegToRad(300) && phi2 < DegToRad(320))
2806  || (phi1 > DegToRad(300) && phi2 > DegToRad(280) && phi1 < DegToRad(320) && phi2 < DegToRad(300))) fhMiDiffSectorDCALPHOSMod[2+etaside]->Fill(pt, m, GetEventWeight());
2807  else if ( (phi1 > DegToRad(260) && phi2 > DegToRad(300) && phi1 < DegToRad(280) && phi2 < DegToRad(320))
2808  || (phi1 > DegToRad(300) && phi2 > DegToRad(260) && phi1 < DegToRad(320) && phi2 < DegToRad(280))) fhMiDiffSectorDCALPHOSMod[4+etaside]->Fill(pt, m, GetEventWeight());
2809  else fhMiDiffSectorDCALPHOSMod[6+etaside]->Fill(pt, m, GetEventWeight());
2810  }
2811  }
2812 
2813  Bool_t ok = kTRUE;
2814  if(fSameSM)
2815  {
2816  if(!fPairWithOtherDetector)
2817  {
2818  if(module1!=module2) ok=kFALSE;
2819  }
2820  else // PHOS and DCal in same sector
2821  {
2822  Float_t phi1 = GetPhi(fPhotonMom1.Phi());
2823  Float_t phi2 = GetPhi(fPhotonMom2.Phi());
2824  ok=kFALSE;
2825  if ( phi1 > DegToRad(260) && phi2 > DegToRad(260) && phi1 < DegToRad(280) && phi2 < DegToRad(280)) ok = kTRUE;
2826  else if ( phi1 > DegToRad(280) && phi2 > DegToRad(280) && phi1 < DegToRad(300) && phi2 < DegToRad(300)) ok = kTRUE;
2827  else if ( phi1 > DegToRad(300) && phi2 > DegToRad(300) && phi1 < DegToRad(320) && phi2 < DegToRad(320)) ok = kTRUE;
2828  }
2829  } // Pair only in same SM
2830 
2831  if(!ok) continue ;
2832 
2833  //
2834  // Do the final histograms with the selected clusters
2835  //
2836 
2837  // Check if one of the clusters comes from a conversion
2838  if(fCheckConversion)
2839  {
2840  if (p1->IsTagged() && p2->IsTagged()) fhMiConv2->Fill(pt, m, GetEventWeight());
2841  else if(p1->IsTagged() || p2->IsTagged()) fhMiConv ->Fill(pt, m, GetEventWeight());
2842  }
2843 
2844  // Fill histograms for different bad channel distance, centrality, assymmetry cut and pid bit
2845  for(Int_t ipid=0; ipid<fNPIDBits; ipid++)
2846  {
2847  if((p1->IsPIDOK(ipid,AliCaloPID::kPhoton)) && (p2->IsPIDOK(ipid,AliCaloPID::kPhoton)))
2848  {
2849  for(Int_t iasym=0; iasym < fNAsymCuts; iasym++)
2850  {
2851  if(a < fAsymCuts[iasym])
2852  {
2853  Int_t index = ((curCentrBin*fNPIDBits)+ipid)*fNAsymCuts + iasym;
2854 
2855  if(index < 0 || index >= ncentr*fNPIDBits*fNAsymCuts) continue ;
2856 
2857  fhMi1[index]->Fill(pt, m, GetEventWeight()) ;
2858 
2859  if(fMakeInvPtPlots)fhMiInvPt1[index]->Fill(pt, m, 1./pt * GetEventWeight()) ;
2860 
2861  if(fFillBadDistHisto)
2862  {
2863  if(p1->DistToBad()>0 && p2->DistToBad()>0)
2864  {
2865  fhMi2[index]->Fill(pt, m, GetEventWeight()) ;
2866  if(fMakeInvPtPlots)fhMiInvPt2[index]->Fill(pt, m, 1./pt * GetEventWeight()) ;
2867 
2868  if(p1->DistToBad()>1 && p2->DistToBad()>1)
2869  {
2870  fhMi3[index]->Fill(pt, m, GetEventWeight()) ;
2871  if(fMakeInvPtPlots)fhMiInvPt3[index]->Fill(pt, m, 1./pt * GetEventWeight()) ;
2872  }
2873  }
2874  }// Fill bad dist histo
2875  }//Asymmetry cut
2876  }// Asymmetry loop
2877  }//PID cut
2878  }//loop for histograms
2879 
2880  //-----------------------
2881  // Multi cuts analysis
2882  //-----------------------
2883  if(fMultiCutAna)
2884  {
2885  // Several pt,ncell and asymmetry cuts
2886 
2887  for(Int_t ipt=0; ipt<fNPtCuts; ipt++)
2888  {
2889  for(Int_t icell=0; icell<fNCellNCuts; icell++)
2890  {
2891  for(Int_t iasym=0; iasym<fNAsymCuts; iasym++)
2892  {
2893  Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
2894  if(p1->Pt() > fPtCuts[ipt] && p2->Pt() > fPtCuts[ipt] &&
2895  a < fAsymCuts[iasym] //&&
2896  //p1->GetBtag() >= fCellNCuts[icell] && p2->GetBtag() >= fCellNCuts[icell] // trick, correct it.
2897  )
2898  {
2899  fhMiPtNCellAsymCuts[index]->Fill(pt, m, GetEventWeight()) ;
2900  //printf("ipt %d, icell%d, iasym %d, name %s\n",ipt, icell, iasym, fhRePtNCellAsymCuts[((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym]->GetName());
2901  }
2902  }// pid bit cut loop
2903  }// icell loop
2904  }// pt cut loop
2905  } // Multi cut ana
2906 
2907  // Fill histograms with opening angle
2908  if(fFillAngleHisto)
2909  {
2910  fhMixedOpeningAngle ->Fill(pt, angle, GetEventWeight());
2911  fhMixedCosOpeningAngle->Fill(pt, TMath::Cos(angle), GetEventWeight());
2912  }
2913 
2914  // Check cell time content in cluster
2915  if ( fFillSecondaryCellTiming)
2916  {
2917  if ( p1->GetFiducialArea() == 0 && p2->GetFiducialArea() == 0 )
2918  fhMiSecondaryCellInTimeWindow ->Fill(pt, m, GetEventWeight());
2919 
2920  else if ( p1->GetFiducialArea() != 0 && p2->GetFiducialArea() != 0 )
2921  fhMiSecondaryCellOutTimeWindow->Fill(pt, m, GetEventWeight());
2922  }
2923 
2924  }// second cluster loop
2925  }//first cluster loop
2926  }//loop on mixed events
2927 
2928  //--------------------------------------------------------
2929  // Add the current event to the list of events for mixing
2930  //--------------------------------------------------------
2931 
2932  //TClonesArray *currentEvent = new TClonesArray(*GetInputAODBranch());
2933  TClonesArray *currentEvent = new TClonesArray(*secondLoopInputData);
2934 
2935  // Add current event to buffer and Remove redundant events
2936  if( currentEvent->GetEntriesFast() > 0 )
2937  {
2938  evMixList->AddFirst(currentEvent) ;
2939  currentEvent=0 ; //Now list of particles belongs to buffer and it will be deleted with buffer
2940  if( evMixList->GetSize() >= GetNMaxEvMix() )
2941  {
2942  TClonesArray * tmp = (TClonesArray*) (evMixList->Last()) ;
2943  evMixList->RemoveLast() ;
2944  delete tmp ;
2945  }
2946  }
2947  else
2948  { // empty event
2949  delete currentEvent ;
2950  currentEvent=0 ;
2951  }
2952  }// DoOwnMix
2953 
2954  AliDebug(1,"End fill histograms");
2955 }
2956 
2957 //________________________________________________________________________
2961 //________________________________________________________________________
2962 Int_t AliAnaPi0::GetEventIndex(AliAODPWG4Particle * part, Double_t * vert)
2963 {
2964  Int_t evtIndex = -1 ;
2965 
2966  if(GetReader()->GetDataType()!=AliCaloTrackReader::kMC)
2967  {
2968  if (GetMixedEvent())
2969  {
2970  evtIndex = GetMixedEvent()->EventIndexForCaloCluster(part->GetCaloLabel(0)) ;
2971  GetVertex(vert,evtIndex);
2972 
2973  if(TMath::Abs(vert[2])> GetZvertexCut())
2974  evtIndex = -2 ; //Event can not be used (vertex, centrality,... cuts not fulfilled)
2975  }
2976  else
2977  {
2978  // Single event
2979  GetVertex(vert);
2980 
2981  if(TMath::Abs(vert[2])> GetZvertexCut())
2982  evtIndex = -1 ; //Event can not be used (vertex, centrality,... cuts not fulfilled)
2983  else
2984  evtIndex = 0 ;
2985  }
2986  } // No MC reader
2987  else
2988  {
2989  evtIndex = 0;
2990  vert[0] = 0. ;
2991  vert[1] = 0. ;
2992  vert[2] = 0. ;
2993  }
2994 
2995  return evtIndex ;
2996 }
2997 
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
ClassImp
ClassImp(AliAnalysisTaskTriggerRates) AliAnalysisTaskTriggerRates
Definition: AliAnalysisTaskTriggerRates.cxx: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:602
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:1581
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:700
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:613
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:101
AliAnaPi0::GetEventIndex
Int_t GetEventIndex(AliAODPWG4Particle *part, Double_t *vert)
Definition: AliAnaPi0.cxx:2962
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:1843
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:604
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:2179
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
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:1906
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