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