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