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