AliPhysics  a34469b (a34469b)
 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), fMultiCutAnaAcc(kFALSE),
58 fNPtCuts(0), fNAsymCuts(0), fNCellNCuts(0), fNPIDBits(0), fNAngleCutBins(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), fFillOpAngleCutHisto(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 fhRePtNCellAsymCutsOpAngle(0x0), fhMiPtNCellAsymCutsOpAngle(0x0), fhRePtNCellAsymCutsSMOpAngle(),
82 /*fhRePIDBits(0x0), fhRePtMult(0x0),*/ fhReSS(),
83 fhRePtAsym(0x0), fhRePtAsymPi0(0x0), fhRePtAsymEta(0x0),
84 fhMiPtAsym(0x0), fhMiPtAsymPi0(0x0), fhMiPtAsymEta(0x0),
85 fhEventBin(0), fhEventMixBin(0),
86 fhCentrality(0x0), fhCentralityNoPair(0x0),
87 fhEventPlaneResolution(0x0),
88 fhRealOpeningAngle(0x0), fhRealCosOpeningAngle(0x0), fhMixedOpeningAngle(0x0), fhMixedCosOpeningAngle(0x0),
89 fhRealOpeningAnglePerSM(),fhMixedOpeningAnglePerSM(),
90 // MC histograms
91 fhPrimPi0E(0x0), fhPrimPi0Pt(0x0),
92 fhPrimPi0AccE(0x0), fhPrimPi0AccPt(0x0), fhPrimPi0AccPtPhotonCuts(0x0),
93 fhPrimPi0Y(0x0), fhPrimPi0AccY(0x0),
94 fhPrimPi0Yeta(0x0), fhPrimPi0YetaYcut(0x0), fhPrimPi0AccYeta(0x0),
95 fhPrimPi0Phi(0x0), fhPrimPi0AccPhi(0x0),
96 fhPrimPi0OpeningAngle(0x0), fhPrimPi0OpeningAnglePhotonCuts(0x0),
97 fhPrimPi0OpeningAngleAsym(0x0),fhPrimPi0CosOpeningAngle(0x0),
98 fhPrimPi0PtCentrality(0), fhPrimPi0PtEventPlane(0),
99 fhPrimPi0AccPtCentrality(0), fhPrimPi0AccPtEventPlane(0),
100 fhPrimEtaE(0x0), fhPrimEtaPt(0x0),
101 fhPrimEtaAccE(0x0), fhPrimEtaAccPt(0x0), fhPrimEtaAccPtPhotonCuts(0x0),
102 fhPrimEtaY(0x0), fhPrimEtaAccY(0x0),
103 fhPrimEtaYeta(0x0), fhPrimEtaYetaYcut(0x0), fhPrimEtaAccYeta(0x0),
104 fhPrimEtaPhi(0x0), fhPrimEtaAccPhi(0x0),
105 fhPrimEtaOpeningAngle(0x0), fhPrimEtaOpeningAnglePhotonCuts(0x0),
106 fhPrimEtaOpeningAngleAsym(0x0),fhPrimEtaCosOpeningAngle(0x0),
107 fhPrimEtaPtCentrality(0), fhPrimEtaPtEventPlane(0),
108 fhPrimEtaAccPtCentrality(0), fhPrimEtaAccPtEventPlane(0),
109 fhPrimPi0PtOrigin(0x0), fhPrimEtaPtOrigin(0x0),
110 fhPrimNotResonancePi0PtOrigin(0x0), fhPrimPi0PtStatus(0x0),
111 fhMCOrgMass(), fhMCOrgAsym(), fhMCOrgDeltaEta(), fhMCOrgDeltaPhi(),
112 fhMCPi0MassPtRec(), fhMCPi0MassPtTrue(), fhMCPi0PtTruePtRec(),
113 fhMCEtaMassPtRec(), fhMCEtaMassPtTrue(), fhMCEtaPtTruePtRec(),
114 fhMCPi0PtOrigin(0x0), fhMCEtaPtOrigin(0x0),
115 fhMCNotResonancePi0PtOrigin(0),fhMCPi0PtStatus(0x0),
116 fhMCPi0ProdVertex(0), fhMCEtaProdVertex(0),
117 fhPrimPi0ProdVertex(0), fhPrimEtaProdVertex(0),
118 fhReMCFromConversion(0), fhReMCFromNotConversion(0), fhReMCFromMixConversion(0),
119 fhCosThStarPrimPi0(0), fhCosThStarPrimEta(0),
120 fhEPairDiffTime(0),
121 fhReSecondaryCellInTimeWindow(0), fhMiSecondaryCellInTimeWindow(0),
122 fhReSecondaryCellOutTimeWindow(0), fhMiSecondaryCellOutTimeWindow(0)
123 
124 {
125  InitParameters();
126 
127  for(Int_t i = 0; i < 4; i++)
128  {
129  fhArmPrimEta[i] = 0;
130  fhArmPrimPi0[i] = 0;
131  }
132 
133  for(Int_t ism = 0; ism < 20; ism++)
134  {
135  fhRealOpeningAnglePerSM [ism] = 0;
136  fhMixedOpeningAnglePerSM[ism] = 0;
137  }
138 
139  for(Int_t icut = 0; icut < 10; icut++)
140  {
141  fhReOpAngleBinMinClusterEtaPhi [icut] = 0;
142  fhReOpAngleBinMaxClusterEtaPhi [icut] = 0;
143  fhReOpAngleBinMinClusterColRow [icut] = 0;
144  fhReOpAngleBinMaxClusterColRow [icut] = 0;
145  fhReOpAngleBinMinClusterEPerSM [icut] = 0;
146  fhReOpAngleBinMaxClusterEPerSM [icut] = 0;
152  fhReOpAngleBinPairClusterMass [icut] = 0;
154 // fhReOpAngleBinPairClusterAbsIdMaxCell[icut] = 0;
155 
158  fhPrimEtaAccPtOpAngCuts [icut] = 0;
159  fhPrimPi0AccPtOpAngCuts [icut] = 0;
160 
161  fhMiOpAngleBinMinClusterEtaPhi [icut] = 0;
162  fhMiOpAngleBinMaxClusterEtaPhi [icut] = 0;
163 // fhMiColRowClusterMinOpAngleBin [icut] = 0;
164 // fhMiOpAngleBinMaxClusterColRow [icut] = 0;
165  fhMiOpAngleBinMinClusterEPerSM [icut] = 0;
166  fhMiOpAngleBinMaxClusterEPerSM [icut] = 0;
172  fhMiOpAngleBinPairClusterMass [icut] = 0;
174 // fhMiOpAngleBinPairClusterAbsIdMaxCell[icut] = 0;
175 
176  fhPtBinClusterEtaPhi [icut] = 0;
177  fhPtBinClusterColRow [icut] = 0;
178  }
179 }
180 
181 //_____________________
183 //_____________________
185 {
186  // Remove event containers
187 
188  if(DoOwnMix() && fEventsList)
189  {
190  for(Int_t ic=0; ic<GetNCentrBin(); ic++)
191  {
192  for(Int_t iz=0; iz<GetNZvertBin(); iz++)
193  {
194  for(Int_t irp=0; irp<GetNRPBin(); irp++)
195  {
196  Int_t bin = GetEventMixBin(ic,iz,irp);
197  fEventsList[bin]->Delete() ;
198  delete fEventsList[bin] ;
199  }
200  }
201  }
202  delete[] fEventsList;
203  }
204 }
205 
206 //______________________________
209 //______________________________
211 {
212  SetInputAODName("PWG4Particle");
213 
214  AddToHistogramsName("AnaPi0_");
215 
216  fUseAngleEDepCut = kFALSE;
217 
218  fUseAngleCut = kTRUE;
219  fAngleCut = 0.;
220  fAngleMaxCut = DegToRad(80.); // 80 degrees cut, avoid EMCal/DCal combinations
221 
222  fMultiCutAna = kFALSE;
223  fMultiCutAnaAcc = kFALSE;
224  fMultiCutAnaSim = kFALSE;
225 
226  fNPtCuts = 3;
227  fPtCuts[0] = 0.; fPtCuts[1] = 0.3; fPtCuts[2] = 0.5;
228  for(Int_t i = fNPtCuts; i < 10; i++) fPtCuts[i] = 0.;
229  for(Int_t i = 0 ; i < 10; i++) fPtCutsMax[i] = 20.;
230 
231  fNAsymCuts = 2;
232  fAsymCuts[0] = 1.; fAsymCuts[1] = 0.7; //fAsymCuts[2] = 0.6; // fAsymCuts[3] = 0.1;
233  for(Int_t i = fNAsymCuts; i < 10; i++)fAsymCuts[i] = 0.;
234 
235  fNCellNCuts = 3;
236  fCellNCuts[0] = 0; fCellNCuts[1] = 1; fCellNCuts[2] = 2;
237  for(Int_t i = fNCellNCuts; i < 10; i++)fCellNCuts[i] = 0;
238 
239  fNPIDBits = 2;
240  fPIDBits[0] = 0; fPIDBits[1] = 2; // fPIDBits[2] = 4; fPIDBits[3] = 6;// check, no cut, dispersion, neutral, dispersion&&neutral
241  for(Int_t i = fNPIDBits; i < 10; i++)fPIDBits[i] = 0;
242 
243 // fNAngleCutBins = 7;
244 // fAngleCutBinsArray[0] = 0.014; fAngleCutBinsArray[1] = 0.035; fAngleCutBinsArray[2] = 0.07; fAngleCutBinsArray[3] = 0.5;
245 // fAngleCutBinsArray[4] = 0.95 ; fAngleCutBinsArray[5] = 1.03 ; fAngleCutBinsArray[6] = 1.1 ; fAngleCutBinsArray[7] = 1.4 ;
246 // for(Int_t i = fNAngleCutBins+1; i < 11; i++)fAngleCutBinsArray[i] = 1000;
247  fNAngleCutBins = 10;
248  Float_t cellsize = 0.0143;
249  fAngleCutBinsArray[0] = 0; fAngleCutBinsArray[1] = 1*cellsize; fAngleCutBinsArray[2] = 0.02; // 1.5*cellsize
250  fAngleCutBinsArray[3] = 2*cellsize; fAngleCutBinsArray[4] = 3*cellsize; fAngleCutBinsArray[5] = 6*cellsize;
251  fAngleCutBinsArray[6] = 12*cellsize; fAngleCutBinsArray[7] = 24*cellsize; fAngleCutBinsArray[8] = 48*cellsize;
252  fAngleCutBinsArray[9] = 96*cellsize; fAngleCutBinsArray[10]= 128*cellsize;
253 }
254 
255 //_______________________________________
257 //_______________________________________
259 {
260  TString parList ; //this will be list of parameters used for this analysis.
261  const Int_t buffersize = 255;
262  char onePar[buffersize] ;
263  snprintf(onePar,buffersize,"--- AliAnaPi0 ---:") ;
264  parList+=onePar ;
265  snprintf(onePar,buffersize,"Number of bins in Centrality: %d;",GetNCentrBin()) ;
266  parList+=onePar ;
267  snprintf(onePar,buffersize,"Number of bins in Z vert. pos: %d;",GetNZvertBin()) ;
268  parList+=onePar ;
269  snprintf(onePar,buffersize,"Number of bins in Reac. Plain: %d;",GetNRPBin()) ;
270  parList+=onePar ;
271  snprintf(onePar,buffersize,"Depth of event buffer: %d;",GetNMaxEvMix()) ;
272  parList+=onePar ;
273  snprintf(onePar,buffersize,"Select pairs with their angle: %d, edep %d, min angle %2.3f, max angle %2.3f;",fUseAngleCut, fUseAngleEDepCut,fAngleCut,fAngleMaxCut) ;
274  parList+=onePar ;
275  snprintf(onePar,buffersize," Asymmetry cuts: n = %d, asymmetry < ",fNAsymCuts) ;
276  for(Int_t i = 0; i < fNAsymCuts; i++) snprintf(onePar,buffersize,"%s %2.2f;",onePar,fAsymCuts[i]);
277  parList+=onePar ;
278  snprintf(onePar,buffersize," PID selection bits: n = %d, PID bit =",fNPIDBits) ;
279  for(Int_t i = 0; i < fNPIDBits; i++) snprintf(onePar,buffersize,"%s %d;",onePar,fPIDBits[i]);
280  parList+=onePar ;
281  snprintf(onePar,buffersize,"Cuts:") ;
282  parList+=onePar ;
283  snprintf(onePar,buffersize,"Z vertex position: -%f < z < %f;",GetZvertexCut(),GetZvertexCut()) ;
284  parList+=onePar ;
285  snprintf(onePar,buffersize,"Calorimeter: %s;",GetCalorimeterString().Data()) ;
286  parList+=onePar ;
287  snprintf(onePar,buffersize,"Number of modules: %d:",fNModules) ;
288  parList+=onePar ;
290  {
291  snprintf(onePar, buffersize," pT cuts: n = %d, pt > ",fNPtCuts) ;
292  for(Int_t i = 0; i < fNPtCuts; i++) snprintf(onePar,buffersize,"%s %2.2f;",onePar,fPtCuts[i]);
293  parList+=onePar ;
294  snprintf(onePar,buffersize, " N cell in cluster cuts: n = %d, nCell > ",fNCellNCuts) ;
295  for(Int_t i = 0; i < fNCellNCuts; i++) snprintf(onePar,buffersize,"%s %d;",onePar,fCellNCuts[i]);
296  parList+=onePar ;
297  }
298 
299  return new TObjString(parList) ;
300 }
301 
302 //_________________________________________
305 //_________________________________________
307 {
308  TList * outputContainer = new TList() ;
309  outputContainer->SetName(GetName());
310 
311  // Set sizes and ranges
312  const Int_t buffersize = 255;
313  char key[buffersize] ;
314  char title[buffersize] ;
315 
317  Int_t nphibins = GetHistogramRanges()->GetHistoPhiBins();
318  Int_t netabins = GetHistogramRanges()->GetHistoEtaBins();
325 
326  Int_t nmassbins = GetHistogramRanges()->GetHistoMassBins();
332 // Int_t ntrmbins = GetHistogramRanges()->GetHistoTrackMultiplicityBins();
333 // Int_t ntrmmax = GetHistogramRanges()->GetHistoTrackMultiplicityMax();
334 // Int_t ntrmmin = GetHistogramRanges()->GetHistoTrackMultiplicityMin();
338  Int_t ntimebins = GetHistogramRanges()->GetHistoTimeBins();
341 
345 
346  // Start with pure MC kinematics histograms
347  // In case other tasks just need this info like AliAnaPi0EbE
348  if(IsDataMC())
349  {
350  // Pi0
351 
352  fhPrimPi0E = new TH1F("hPrimPi0E","Primary #pi^{0} E, |#it{Y}|<1",
353  nptbins,ptmin,ptmax) ;
354  fhPrimPi0E ->SetXTitle("#it{E} (GeV)");
355  outputContainer->Add(fhPrimPi0E) ;
356 
357  fhPrimPi0Pt = new TH1F("hPrimPi0Pt","Primary #pi^{0} #it{p}_{T} , |#it{Y}|<1",
358  nptbins,ptmin,ptmax) ;
359  fhPrimPi0Pt ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
360  outputContainer->Add(fhPrimPi0Pt) ;
361 
362  Int_t netabinsopen = TMath::Nint(netabins*4/(etamax-etamin));
363  fhPrimPi0Y = new TH2F("hPrimPi0Rapidity","Rapidity of primary #pi^{0}",
364  nptbins,ptmin,ptmax,netabinsopen,-2, 2) ;
365  fhPrimPi0Y ->SetYTitle("#it{Rapidity}");
366  fhPrimPi0Y ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
367  outputContainer->Add(fhPrimPi0Y) ;
368 
369  fhPrimPi0Yeta = new TH2F("hPrimPi0PseudoRapidity","PseudoRapidity of primary #pi^{0}",
370  nptbins,ptmin,ptmax,netabinsopen,-2, 2) ;
371  fhPrimPi0Yeta ->SetYTitle("#eta");
372  fhPrimPi0Yeta ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
373  outputContainer->Add(fhPrimPi0Yeta) ;
374 
375  fhPrimPi0YetaYcut = new TH2F("hPrimPi0PseudoRapidityYcut","PseudoRapidity of primary #pi^{0}, |#it{Y}|<1",
376  nptbins,ptmin,ptmax,netabinsopen,-2, 2) ;
377  fhPrimPi0YetaYcut ->SetYTitle("#eta");
378  fhPrimPi0YetaYcut ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
379  outputContainer->Add(fhPrimPi0YetaYcut) ;
380 
381  Int_t nphibinsopen = TMath::Nint(nphibins*TMath::TwoPi()/(phimax-phimin));
382  fhPrimPi0Phi = new TH2F("hPrimPi0Phi","#phi of primary #pi^{0}, |#it{Y}|<1",
383  nptbins,ptmin,ptmax,nphibinsopen,0,360) ;
384  fhPrimPi0Phi->SetYTitle("#phi (deg)");
385  fhPrimPi0Phi->SetXTitle("#it{p}_{T} (GeV/#it{c})");
386  outputContainer->Add(fhPrimPi0Phi) ;
387 
389  {
390  fhPrimPi0AccE = new TH1F("hPrimPi0AccE","Primary #pi^{0} #it{E} with both photons in acceptance",
391  nptbins,ptmin,ptmax) ;
392  fhPrimPi0AccE->SetXTitle("#it{E} (GeV)");
393  outputContainer->Add(fhPrimPi0AccE) ;
394 
395  fhPrimPi0AccPt = new TH1F("hPrimPi0AccPt","Primary #pi^{0} #it{p}_{T} with both photons in acceptance",
396  nptbins,ptmin,ptmax) ;
397  fhPrimPi0AccPt->SetXTitle("#it{p}_{T} (GeV/#it{c})");
398  outputContainer->Add(fhPrimPi0AccPt) ;
399 
400  fhPrimPi0AccPtPhotonCuts = new TH1F("hPrimPi0AccPtPhotonCuts","Primary #pi^{0} #it{p}_{T} with both photons in acceptance",
401  nptbins,ptmin,ptmax) ;
402  fhPrimPi0AccPtPhotonCuts->SetXTitle("#it{p}_{T} (GeV/#it{c})");
403  outputContainer->Add(fhPrimPi0AccPtPhotonCuts) ;
404 
405  fhPrimPi0AccY = new TH2F("hPrimPi0AccRapidity","Rapidity of primary #pi^{0} with accepted daughters",
406  nptbins,ptmin,ptmax,netabins,etamin,etamax) ;
407  fhPrimPi0AccY->SetYTitle("Rapidity");
408  fhPrimPi0AccY->SetXTitle("#it{p}_{T} (GeV/#it{c})");
409  outputContainer->Add(fhPrimPi0AccY) ;
410 
411  fhPrimPi0AccYeta = new TH2F("hPrimPi0AccPseudoRapidity","PseudoRapidity of primary #pi^{0} with accepted daughters",
412  nptbins,ptmin,ptmax,netabins,etamin,etamax) ;
413  fhPrimPi0AccYeta ->SetYTitle("#eta");
414  fhPrimPi0AccYeta ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
415  outputContainer->Add(fhPrimPi0AccYeta) ;
416 
417  fhPrimPi0AccPhi = new TH2F("hPrimPi0AccPhi","#phi of primary #pi^{0} with accepted daughters",
418  nptbins,ptmin,ptmax,
419  nphibins,phimin*TMath::RadToDeg(),phimax*TMath::RadToDeg()) ;
420  fhPrimPi0AccPhi->SetYTitle("#phi (deg)");
421  fhPrimPi0AccPhi->SetXTitle("#it{p}_{T} (GeV/#it{c})");
422  outputContainer->Add(fhPrimPi0AccPhi) ;
423  }
424 
425  // Eta
426 
427  fhPrimEtaE = new TH1F("hPrimEtaE","Primary eta E",
428  nptbins,ptmin,ptmax) ;
429  fhPrimEtaE ->SetXTitle("#it{E} (GeV)");
430  outputContainer->Add(fhPrimEtaAccE) ;
431 
432  fhPrimEtaPt = new TH1F("hPrimEtaPt","Primary #eta #it{p}_{T}",
433  nptbins,ptmin,ptmax) ;
434  fhPrimEtaPt ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
435  outputContainer->Add(fhPrimEtaPt) ;
436 
437  fhPrimEtaY = new TH2F("hPrimEtaRapidity","Rapidity of primary #eta",
438  nptbins,ptmin,ptmax,netabinsopen,-2, 2) ;
439  fhPrimEtaY->SetYTitle("#it{Rapidity}");
440  fhPrimEtaY->SetXTitle("#it{p}_{T} (GeV/#it{c})");
441  outputContainer->Add(fhPrimEtaY) ;
442 
443  fhPrimEtaYeta = new TH2F("hPrimEtaPseudoRapidityEta","PseudoRapidity of primary #eta",
444  nptbins,ptmin,ptmax,netabinsopen,-2, 2) ;
445  fhPrimEtaYeta->SetYTitle("#it{Rapidity}");
446  fhPrimEtaYeta->SetXTitle("#it{p}_{T} (GeV/#it{c})");
447  outputContainer->Add(fhPrimEtaYeta) ;
448 
449  fhPrimEtaYetaYcut = new TH2F("hPrimEtaPseudoRapidityEtaYcut","PseudoRapidity of primary #eta, |#it{Y}|<1",
450  nptbins,ptmin,ptmax,netabinsopen,-2, 2) ;
451  fhPrimEtaYetaYcut->SetYTitle("#it{Pseudorapidity}");
452  fhPrimEtaYetaYcut->SetXTitle("#it{p}_{T} (GeV/#it{c})");
453  outputContainer->Add(fhPrimEtaYetaYcut) ;
454 
455  fhPrimEtaPhi = new TH2F("hPrimEtaPhi","Azimuthal of primary #eta",
456  nptbins,ptmin,ptmax, nphibinsopen,0,360) ;
457  fhPrimEtaPhi->SetYTitle("#phi (deg)");
458  fhPrimEtaPhi->SetXTitle("#it{p}_{T} (GeV/#it{c})");
459  outputContainer->Add(fhPrimEtaPhi) ;
460 
462  {
463  fhPrimEtaAccE = new TH1F("hPrimEtaAccE","Primary #eta #it{E} with both photons in acceptance",
464  nptbins,ptmin,ptmax) ;
465  fhPrimEtaAccE->SetXTitle("#it{E} (GeV)");
466  outputContainer->Add(fhPrimEtaE) ;
467 
468  fhPrimEtaAccPt = new TH1F("hPrimEtaAccPt","Primary eta #it{p}_{T} with both photons in acceptance",
469  nptbins,ptmin,ptmax) ;
470  fhPrimEtaAccPt->SetXTitle("#it{p}_{T} (GeV/#it{c})");
471  outputContainer->Add(fhPrimEtaAccPt) ;
472 
473  fhPrimEtaAccPtPhotonCuts = new TH1F("hPrimEtaAccPtPhotonCuts","Primary eta #it{p}_{T} with both photons in acceptance",
474  nptbins,ptmin,ptmax) ;
475  fhPrimEtaAccPtPhotonCuts->SetXTitle("#it{p}_{T} (GeV/#it{c})");
476  outputContainer->Add(fhPrimEtaAccPtPhotonCuts) ;
477 
478  fhPrimEtaAccPhi = new TH2F("hPrimEtaAccPhi","Azimuthal of primary #eta with accepted daughters",
479  nptbins,ptmin,ptmax, nphibins,phimin*TMath::RadToDeg(),phimax*TMath::RadToDeg()) ;
480  fhPrimEtaAccPhi->SetYTitle("#phi (deg)");
481  fhPrimEtaAccPhi->SetXTitle("#it{p}_{T} (GeV/#it{c})");
482  outputContainer->Add(fhPrimEtaAccPhi) ;
483 
484  fhPrimEtaAccY = new TH2F("hPrimEtaAccRapidity","Rapidity of primary #eta",
485  nptbins,ptmin,ptmax, netabins,etamin,etamax) ;
486  fhPrimEtaAccY->SetYTitle("#it{Rapidity}");
487  fhPrimEtaAccY->SetXTitle("#it{p}_{T} (GeV/#it{c})");
488  outputContainer->Add(fhPrimEtaAccY) ;
489 
490  fhPrimEtaAccYeta = new TH2F("hPrimEtaAccPseudoRapidity","PseudoRapidity of primary #eta",
491  nptbins,ptmin,ptmax, netabins,etamin,etamax) ;
492  fhPrimEtaAccYeta->SetYTitle("#it{Pseudorapidity}");
493  fhPrimEtaAccYeta->SetXTitle("#it{p}_{T} (GeV/#it{c})");
494  outputContainer->Add(fhPrimEtaAccYeta) ;
495  }
496 
497  // Create histograms only for PbPb or high multiplicity analysis analysis
499  {
500  fhPrimPi0PtCentrality = new TH2F("hPrimPi0PtCentrality","Primary #pi^{0} #it{p}_{T} vs reco centrality, |#it{Y}|<1",
501  nptbins,ptmin,ptmax, 100, 0, 100) ;
502  fhPrimPi0PtCentrality ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
503  fhPrimPi0PtCentrality ->SetYTitle("Centrality");
504  outputContainer->Add(fhPrimPi0PtCentrality) ;
505 
506  fhPrimEtaPtCentrality = new TH2F("hPrimEtaPtCentrality","Primary #eta #it{p}_{T} vs reco centrality, |#it{Y}|<1",
507  nptbins,ptmin,ptmax, 100, 0, 100) ;
508  fhPrimEtaPtCentrality ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
509  fhPrimEtaPtCentrality ->SetYTitle("Centrality");
510  outputContainer->Add(fhPrimEtaPtCentrality) ;
511 
512 
513  fhPrimPi0PtEventPlane = new TH2F("hPrimPi0PtEventPlane","Primary #pi^{0} #it{p}_{T} vs reco event plane angle, |#it{Y}|<1",
514  nptbins,ptmin,ptmax, 100, 0, TMath::Pi()) ;
515  fhPrimPi0PtEventPlane ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
516  fhPrimPi0PtEventPlane ->SetYTitle("Event Plane Angle (rad)");
517  outputContainer->Add(fhPrimPi0PtEventPlane) ;
518 
519 
520  fhPrimEtaPtEventPlane = new TH2F("hPrimEtaPtEventPlane","Primary #eta #it{p}_{T} vs reco event plane angle, |#it{Y}|<1",
521  nptbins,ptmin,ptmax, 100, 0, TMath::Pi()) ;
522  fhPrimEtaPtEventPlane ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
523  fhPrimEtaPtEventPlane ->SetYTitle("Event Plane Angle (rad)");
524  outputContainer->Add(fhPrimEtaPtEventPlane) ;
525 
527  {
528  fhPrimPi0AccPtCentrality = new TH2F("hPrimPi0AccPtCentrality","Primary #pi^{0} with both photons in acceptance #it{p}_{T} vs reco centrality",
529  nptbins,ptmin,ptmax, 100, 0, 100) ;
530  fhPrimPi0AccPtCentrality->SetXTitle("#it{p}_{T} (GeV/#it{c})");
531  fhPrimPi0AccPtCentrality->SetYTitle("Centrality");
532  outputContainer->Add(fhPrimPi0AccPtCentrality) ;
533 
534  fhPrimEtaAccPtCentrality = new TH2F("hPrimEtaAccPtCentrality","Primary #eta with both photons in acceptance #it{p}_{T} vs reco centrality",
535  nptbins,ptmin,ptmax, 100, 0, 100) ;
536  fhPrimEtaAccPtCentrality->SetXTitle("#it{p}_{T} (GeV/#it{c})");
537  fhPrimEtaAccPtCentrality->SetYTitle("Centrality");
538  outputContainer->Add(fhPrimEtaAccPtCentrality) ;
539 
540  fhPrimPi0AccPtEventPlane = new TH2F("hPrimPi0AccPtEventPlane","Primary #pi^{0} with both photons in acceptance #it{p}_{T} vs reco event plane angle",
541  nptbins,ptmin,ptmax, 100, 0, TMath::Pi()) ;
542  fhPrimPi0AccPtEventPlane->SetXTitle("#it{p}_{T} (GeV/#it{c})");
543  fhPrimPi0AccPtEventPlane->SetYTitle("Event Plane Angle (rad)");
544  outputContainer->Add(fhPrimPi0AccPtEventPlane) ;
545 
546  fhPrimEtaAccPtEventPlane = new TH2F("hPrimEtaAccPtEventPlane","Primary #eta with both #gamma_{decay} in acceptance #it{p}_{T} vs reco event plane angle",
547  nptbins,ptmin,ptmax, 100, 0, TMath::Pi()) ;
548  fhPrimEtaAccPtEventPlane->SetXTitle("#it{p}_{T} (GeV/#it{c})");
549  fhPrimEtaAccPtEventPlane->SetYTitle("Event Plane Angle (rad)");
550  outputContainer->Add(fhPrimEtaAccPtEventPlane) ;
551  }
552  }
553 
555  {
557  ("hPrimPi0OpeningAngle","Angle between all primary #gamma pair vs E_{#pi^{0}}, in acceptance",
558  nptbins,ptmin,ptmax,nopanbins,opanmin,opanmax);
559  fhPrimPi0OpeningAngle->SetYTitle("#theta(rad)");
560  fhPrimPi0OpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
561  outputContainer->Add(fhPrimPi0OpeningAngle) ;
562 
564  ("hPrimPi0OpeningAnglePhotonCuts","Angle between all primary #gamma pair vs E_{#pi^{0}} in acceptance",
565  nptbins,ptmin,ptmax,nopanbins,opanmin,opanmax);
566  fhPrimPi0OpeningAnglePhotonCuts->SetYTitle("#theta(rad)");
567  fhPrimPi0OpeningAnglePhotonCuts->SetXTitle("E_{ #pi^{0}} (GeV)");
568  outputContainer->Add(fhPrimPi0OpeningAnglePhotonCuts) ;
569 
571  ("hPrimPi0OpeningAngleAsym","Angle between all primary #gamma pair vs #it{Asymmetry}, in acceptance, #it{p}_{T}>5 GeV/#it{c}",
572  100,0,1,nopanbins,opanmin,opanmax);
573  fhPrimPi0OpeningAngleAsym->SetXTitle("|A|=| (E_{1}-E_{2}) / (E_{1}+E_{2}) |");
574  fhPrimPi0OpeningAngleAsym->SetYTitle("#theta(rad)");
575  outputContainer->Add(fhPrimPi0OpeningAngleAsym) ;
576 
578  ("hPrimPi0CosOpeningAngle","Cosinus of angle between all primary #gamma pair vs E_{#pi^{0}}, in acceptance",
579  nptbins,ptmin,ptmax,100,-1,1);
580  fhPrimPi0CosOpeningAngle->SetYTitle("cos (#theta) ");
581  fhPrimPi0CosOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
582  outputContainer->Add(fhPrimPi0CosOpeningAngle) ;
583 
585  ("hPrimEtaOpeningAngle","Angle between all primary #gamma pair vs E_{#eta}, in acceptance",
586  nptbins,ptmin,ptmax,nopanbins,opanmin,opanmax);
587  fhPrimEtaOpeningAngle->SetYTitle("#theta(rad)");
588  fhPrimEtaOpeningAngle->SetXTitle("E_{#eta} (GeV)");
589  outputContainer->Add(fhPrimEtaOpeningAngle) ;
590 
592  ("hPrimEtaOpeningAnglePhotonCuts","Angle between all primary #gamma pair vs E_{#eta}, in acceptance",
593  nptbins,ptmin,ptmax,nopanbins,opanmin,opanmax);
594  fhPrimEtaOpeningAnglePhotonCuts->SetYTitle("#theta(rad)");
595  fhPrimEtaOpeningAnglePhotonCuts->SetXTitle("E_{#eta} (GeV)");
596  outputContainer->Add(fhPrimEtaOpeningAnglePhotonCuts) ;
597 
599  ("hPrimEtaOpeningAngleAsym","Angle between all primary #gamma pair vs #it{Asymmetry}, #it{p}_{T}>5 GeV/#it{c}, in acceptance",
600  100,0,1,nopanbins,opanmin,opanmax);
601  fhPrimEtaOpeningAngleAsym->SetXTitle("|#it{A}|=| (#it{E}_{1}-#it{E}_{2}) / (#it{E}_{1}+#it{E}_{2}) |");
602  fhPrimEtaOpeningAngleAsym->SetYTitle("#theta(rad)");
603  outputContainer->Add(fhPrimEtaOpeningAngleAsym) ;
604 
606  ("hPrimEtaCosOpeningAngle","Cosinus of angle between all primary #gamma pair vs E_{#eta}, in acceptance",
607  nptbins,ptmin,ptmax,100,-1,1);
608  fhPrimEtaCosOpeningAngle->SetYTitle("cos (#theta) ");
609  fhPrimEtaCosOpeningAngle->SetXTitle("#it{E}_{ #eta} (GeV)");
610  outputContainer->Add(fhPrimEtaCosOpeningAngle) ;
611  }
612 
613  // Primary origin
614  if(fFillOriginHisto)
615  {
616  // Pi0
617  fhPrimPi0PtOrigin = new TH2F("hPrimPi0PtOrigin","Primary #pi^{0} #it{p}_{T} vs origin",nptbins,ptmin,ptmax,11,0,11) ;
618  fhPrimPi0PtOrigin->SetXTitle("#it{p}_{T} (GeV/#it{c})");
619  fhPrimPi0PtOrigin->SetYTitle("Origin");
620  fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(1 ,"Status 21");
621  fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(2 ,"Quark");
622  fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(3 ,"qq Resonances ");
623  fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(4 ,"Resonances");
624  fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(5 ,"#rho");
625  fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(6 ,"#omega");
626  fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(7 ,"K");
627  fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(8 ,"Other");
628  fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(9 ,"#eta");
629  fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(10 ,"#eta prime");
630  outputContainer->Add(fhPrimPi0PtOrigin) ;
631 
632  fhPrimNotResonancePi0PtOrigin = new TH2F("hPrimNotResonancePi0PtOrigin","Primary #pi^{0} #it{p}_{T} vs origin",nptbins,ptmin,ptmax,11,0,11) ;
633  fhPrimNotResonancePi0PtOrigin->SetXTitle("#it{p}_{T} (GeV/#it{c})");
634  fhPrimNotResonancePi0PtOrigin->SetYTitle("Origin");
635  fhPrimNotResonancePi0PtOrigin->GetYaxis()->SetBinLabel(1 ,"Status 21");
636  fhPrimNotResonancePi0PtOrigin->GetYaxis()->SetBinLabel(2 ,"Quark");
637  fhPrimNotResonancePi0PtOrigin->GetYaxis()->SetBinLabel(3 ,"qq Resonances");
638  fhPrimNotResonancePi0PtOrigin->GetYaxis()->SetBinLabel(4 ,"Resonances");
639  fhPrimNotResonancePi0PtOrigin->GetYaxis()->SetBinLabel(5 ,"#rho");
640  fhPrimNotResonancePi0PtOrigin->GetYaxis()->SetBinLabel(6 ,"#omega");
641  fhPrimNotResonancePi0PtOrigin->GetYaxis()->SetBinLabel(7 ,"K");
642  fhPrimNotResonancePi0PtOrigin->GetYaxis()->SetBinLabel(8 ,"Other");
643  fhPrimNotResonancePi0PtOrigin->GetYaxis()->SetBinLabel(9 ,"#eta");
644  fhPrimNotResonancePi0PtOrigin->GetYaxis()->SetBinLabel(10 ,"#eta prime");
645  outputContainer->Add(fhPrimNotResonancePi0PtOrigin) ;
646 
647  fhPrimPi0PtStatus = new TH2F("hPrimPi0PtStatus","Primary #pi^{0} #it{p}_{T} vs status",nptbins,ptmin,ptmax,101,-50,50) ;
648  fhPrimPi0PtStatus->SetXTitle("#it{p}_{T} (GeV/#it{c})");
649  fhPrimPi0PtStatus->SetYTitle("Status");
650  outputContainer->Add(fhPrimPi0PtStatus) ;
651 
652  // Eta
653  fhPrimEtaPtOrigin = new TH2F("hPrimEtaPtOrigin","Primary #pi^{0} #it{p}_{T} vs origin",nptbins,ptmin,ptmax,7,0,7) ;
654  fhPrimEtaPtOrigin->SetXTitle("#it{p}_{T} (GeV/#it{c})");
655  fhPrimEtaPtOrigin->SetYTitle("Origin");
656  fhPrimEtaPtOrigin->GetYaxis()->SetBinLabel(1 ,"Status 21");
657  fhPrimEtaPtOrigin->GetYaxis()->SetBinLabel(2 ,"Quark");
658  fhPrimEtaPtOrigin->GetYaxis()->SetBinLabel(3 ,"qq Resonances");
659  fhPrimEtaPtOrigin->GetYaxis()->SetBinLabel(4 ,"Resonances");
660  fhPrimEtaPtOrigin->GetYaxis()->SetBinLabel(5 ,"Other");
661  fhPrimEtaPtOrigin->GetYaxis()->SetBinLabel(6 ,"#eta prime ");
662  outputContainer->Add(fhPrimEtaPtOrigin) ;
663 
664  // Production vertex
665  fhPrimPi0ProdVertex = new TH2F("hPrimPi0ProdVertex","generated #pi^{0} #it{p}_{T} vs production vertex",
666  200,0.,20.,5000,0,500) ;
667  fhPrimPi0ProdVertex->SetXTitle("#it{p}_{T} (GeV/#it{c})");
668  fhPrimPi0ProdVertex->SetYTitle("#it{R} (cm)");
669  outputContainer->Add(fhPrimPi0ProdVertex) ;
670 
671  fhPrimEtaProdVertex = new TH2F("hPrimEtaProdVertex","generated #eta #it{p}_{T} vs production vertex",
672  200,0.,20.,5000,0,500) ;
673  fhPrimEtaProdVertex->SetXTitle("#it{p}_{T} (GeV/#it{c})");
674  fhPrimEtaProdVertex->SetYTitle("#it{R} (cm)");
675  outputContainer->Add(fhPrimEtaProdVertex) ;
676  }
677 
679  {
680  TString ebin[] = {"8 < E < 12 GeV","12 < E < 16 GeV", "16 < E < 20 GeV", "E > 20 GeV" };
681  Int_t narmbins = 400;
682  Float_t armmin = 0;
683  Float_t armmax = 0.4;
684 
685  for(Int_t i = 0; i < 4; i++)
686  {
687  fhArmPrimPi0[i] = new TH2F(Form("hArmenterosPrimPi0EBin%d",i),
688  Form("Armenteros of primary #pi^{0}, %s",ebin[i].Data()),
689  200, -1, 1, narmbins,armmin,armmax);
690  fhArmPrimPi0[i]->SetYTitle("#it{p}_{T}^{Arm}");
691  fhArmPrimPi0[i]->SetXTitle("#alpha^{Arm}");
692  outputContainer->Add(fhArmPrimPi0[i]) ;
693 
694  fhArmPrimEta[i] = new TH2F(Form("hArmenterosPrimEtaEBin%d",i),
695  Form("Armenteros of primary #eta, %s",ebin[i].Data()),
696  200, -1, 1, narmbins,armmin,armmax);
697  fhArmPrimEta[i]->SetYTitle("#it{p}_{T}^{Arm}");
698  fhArmPrimEta[i]->SetXTitle("#alpha^{Arm}");
699  outputContainer->Add(fhArmPrimEta[i]) ;
700  }
701 
702  // Same as asymmetry ...
704  ("hCosThStarPrimPi0","cos(#theta *) for primary #pi^{0}",nptbins,ptmin,ptmax,200,-1,1);
705  fhCosThStarPrimPi0->SetYTitle("cos(#theta *)");
706  fhCosThStarPrimPi0->SetXTitle("E_{ #pi^{0}} (GeV)");
707  outputContainer->Add(fhCosThStarPrimPi0) ;
708 
710  ("hCosThStarPrimEta","cos(#theta *) for primary #eta",nptbins,ptmin,ptmax,200,-1,1);
711  fhCosThStarPrimEta->SetYTitle("cos(#theta *)");
712  fhCosThStarPrimEta->SetXTitle("E_{ #eta} (GeV)");
713  outputContainer->Add(fhCosThStarPrimEta) ;
714  }
715 
716  if(fFillOnlyMCAcceptanceHisto) return outputContainer;
717  }
718 
719  //
720  // Create mixed event containers
721  //
723 
724  for(Int_t ic=0; ic<GetNCentrBin(); ic++)
725  {
726  for(Int_t iz=0; iz<GetNZvertBin(); iz++)
727  {
728  for(Int_t irp=0; irp<GetNRPBin(); irp++)
729  {
730  Int_t bin = GetEventMixBin(ic,iz,irp);
731  fEventsList[bin] = new TList() ;
732  fEventsList[bin]->SetOwner(kFALSE);
733  }
734  }
735  }
736 
737  // Init the number of modules, set in the class AliCalorimeterUtils
739  if(GetCalorimeter()==kPHOS && fNModules > 4) fNModules = 4;
740 
741  fhReMod = new TH2F*[fNModules] ;
742  fhMiMod = new TH2F*[fNModules] ;
743 
745  {
746  if(GetCalorimeter() == kPHOS)
747  {
748  fhReDiffPHOSMod = new TH2F*[fNModules] ;
749  fhMiDiffPHOSMod = new TH2F*[fNModules] ;
750  }
751  else
752  {
757  }
758  }
759  else
760  {
761  fhReSameSectorDCALPHOSMod = new TH2F*[6] ;
762  fhReDiffSectorDCALPHOSMod = new TH2F*[8] ;
763  fhMiSameSectorDCALPHOSMod = new TH2F*[6] ;
764  fhMiDiffSectorDCALPHOSMod = new TH2F*[8] ;
765  }
766 
770  {
775  }
776  if(fMakeInvPtPlots)
777  {
780  if(fFillBadDistHisto){
785  }
786  }
787 
789  {
790  fhReSecondaryCellInTimeWindow = new TH2F("hReSecondaryCellInTimeWindow","Real Pair, it{t}_{cell} < 50 ns, w_{cell} > 0",
791  nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
792  fhReSecondaryCellInTimeWindow->SetXTitle("#it{p}_{T} (GeV/#it{c})");
793  fhReSecondaryCellInTimeWindow->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
794  outputContainer->Add(fhReSecondaryCellInTimeWindow) ;
795 
796  fhReSecondaryCellOutTimeWindow = new TH2F("hReSecondaryCellOutTimeWindow","Real Pair, it{t}_{cell} < 50 ns, w_{cell} > 0",
797  nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
798  fhReSecondaryCellOutTimeWindow->SetXTitle("#it{p}_{T} (GeV/#it{c})");
799  fhReSecondaryCellOutTimeWindow->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
800  outputContainer->Add(fhReSecondaryCellOutTimeWindow) ;
801 
802  if ( DoOwnMix() )
803  {
804  fhMiSecondaryCellInTimeWindow = new TH2F("hMiSecondaryCellInTimeWindow","Real Pair, it{t}_{cell} < 50 ns, w_{cell} > 0",
805  nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
806  fhMiSecondaryCellInTimeWindow->SetXTitle("#it{p}_{T} (GeV/#it{c})");
807  fhMiSecondaryCellInTimeWindow->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
808  outputContainer->Add(fhMiSecondaryCellInTimeWindow) ;
809 
810  fhMiSecondaryCellOutTimeWindow = new TH2F("hMiSecondaryCellOutTimeWindow","Real Pair, it{t}_{cell} < 50 ns, w_{cell} > 0",
811  nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
812  fhMiSecondaryCellOutTimeWindow->SetXTitle("#it{p}_{T} (GeV/#it{c})");
813  fhMiSecondaryCellOutTimeWindow->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
814  outputContainer->Add(fhMiSecondaryCellOutTimeWindow) ;
815 
816  }
817  }
818 
819  if ( fCheckConversion )
820  {
821  fhReConv = new TH2F("hReConv","Real Pair with one recombined conversion ",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
822  fhReConv->SetXTitle("#it{p}_{T} (GeV/#it{c})");
823  fhReConv->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
824  outputContainer->Add(fhReConv) ;
825 
826  fhReConv2 = new TH2F("hReConv2","Real Pair with 2 recombined conversion ",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
827  fhReConv2->SetXTitle("#it{p}_{T} (GeV/#it{c})");
828  fhReConv2->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
829  outputContainer->Add(fhReConv2) ;
830 
831  if ( DoOwnMix() )
832  {
833  fhMiConv = new TH2F("hMiConv","Mixed Pair with one recombined conversion ",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
834  fhMiConv->SetXTitle("#it{p}_{T} (GeV/#it{c})");
835  fhMiConv->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
836  outputContainer->Add(fhMiConv) ;
837 
838  fhMiConv2 = new TH2F("hMiConv2","Mixed Pair with 2 recombined conversion ",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
839  fhMiConv2->SetXTitle("#it{p}_{T} (GeV/#it{c})");
840  fhMiConv2->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
841  outputContainer->Add(fhMiConv2) ;
842  }
843  }
844 
845  for(Int_t ic=0; ic<GetNCentrBin(); ic++)
846  {
847  for(Int_t ipid=0; ipid<fNPIDBits; ipid++)
848  {
849  for(Int_t iasym=0; iasym<fNAsymCuts; iasym++)
850  {
851  Int_t index = ((ic*fNPIDBits)+ipid)*fNAsymCuts + iasym;
852  //printf("cen %d, pid %d, asy %d, Index %d\n",ic,ipid,iasym,index);
853  //Distance to bad module 1
854  snprintf(key, buffersize,"hRe_cen%d_pidbit%d_asy%d_dist1",ic,ipid,iasym) ;
855  snprintf(title, buffersize,"Real #it{M}_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 1",
856  ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
857  fhRe1[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
858  fhRe1[index]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
859  fhRe1[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
860  //printf("name: %s\n ",fhRe1[index]->GetName());
861  outputContainer->Add(fhRe1[index]) ;
862 
864  {
865  //Distance to bad module 2
866  snprintf(key, buffersize,"hRe_cen%d_pidbit%d_asy%d_dist2",ic,ipid,iasym) ;
867  snprintf(title, buffersize,"Real #it{M}_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 2",
868  ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
869  fhRe2[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
870  fhRe2[index]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
871  fhRe2[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
872  outputContainer->Add(fhRe2[index]) ;
873 
874  //Distance to bad module 3
875  snprintf(key, buffersize,"hRe_cen%d_pidbit%d_asy%d_dist3",ic,ipid,iasym) ;
876  snprintf(title, buffersize,"Real #it{M}_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 3",
877  ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
878  fhRe3[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
879  fhRe3[index]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
880  fhRe3[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
881  outputContainer->Add(fhRe3[index]) ;
882  }
883 
884  //Inverse pT
885  if(fMakeInvPtPlots)
886  {
887  //Distance to bad module 1
888  snprintf(key, buffersize,"hReInvPt_cen%d_pidbit%d_asy%d_dist1",ic,ipid,iasym) ;
889  snprintf(title, buffersize,"Real #it{M}_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 1",
890  ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
891  fhReInvPt1[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
892  fhReInvPt1[index]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
893  fhReInvPt1[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
894  outputContainer->Add(fhReInvPt1[index]) ;
895 
896  if(fFillBadDistHisto){
897  //Distance to bad module 2
898  snprintf(key, buffersize,"hReInvPt_cen%d_pidbit%d_asy%d_dist2",ic,ipid,iasym) ;
899  snprintf(title, buffersize,"Real #it{M}_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 2",
900  ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
901  fhReInvPt2[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
902  fhReInvPt2[index]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
903  fhReInvPt2[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
904  outputContainer->Add(fhReInvPt2[index]) ;
905 
906  //Distance to bad module 3
907  snprintf(key, buffersize,"hReInvPt_cen%d_pidbit%d_asy%d_dist3",ic,ipid,iasym) ;
908  snprintf(title, buffersize,"Real #it{M}_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 3",
909  ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
910  fhReInvPt3[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
911  fhReInvPt3[index]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
912  fhReInvPt3[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
913  outputContainer->Add(fhReInvPt3[index]) ;
914  }
915  }
916 
917  if(DoOwnMix())
918  {
919  //Distance to bad module 1
920  snprintf(key, buffersize,"hMi_cen%d_pidbit%d_asy%d_dist1",ic,ipid,iasym) ;
921  snprintf(title, buffersize,"Mixed #it{M}_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 1",
922  ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
923  fhMi1[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
924  fhMi1[index]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
925  fhMi1[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
926  outputContainer->Add(fhMi1[index]) ;
927  if(fFillBadDistHisto){
928  //Distance to bad module 2
929  snprintf(key, buffersize,"hMi_cen%d_pidbit%d_asy%d_dist2",ic,ipid,iasym) ;
930  snprintf(title, buffersize,"Mixed #it{M}_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 2",
931  ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
932  fhMi2[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
933  fhMi2[index]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
934  fhMi2[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
935  outputContainer->Add(fhMi2[index]) ;
936 
937  //Distance to bad module 3
938  snprintf(key, buffersize,"hMi_cen%d_pidbit%d_asy%d_dist3",ic,ipid,iasym) ;
939  snprintf(title, buffersize,"Mixed #it{M}_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 3",
940  ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
941  fhMi3[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
942  fhMi3[index]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
943  fhMi3[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
944  outputContainer->Add(fhMi3[index]) ;
945  }
946 
947  // Inverse pT
948  if(fMakeInvPtPlots)
949  {
950  // Distance to bad module 1
951  snprintf(key, buffersize,"hMiInvPt_cen%d_pidbit%d_asy%d_dist1",ic,ipid,iasym) ;
952  snprintf(title, buffersize,"Mixed #it{M}_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 1",
953  ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
954  fhMiInvPt1[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
955  fhMiInvPt1[index]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
956  fhMiInvPt1[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
957  outputContainer->Add(fhMiInvPt1[index]) ;
958  if(fFillBadDistHisto){
959  // Distance to bad module 2
960  snprintf(key, buffersize,"hMiInvPt_cen%d_pidbit%d_asy%d_dist2",ic,ipid,iasym) ;
961  snprintf(title, buffersize,"Mixed #it{M}_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 2",
962  ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
963  fhMiInvPt2[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
964  fhMiInvPt2[index]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
965  fhMiInvPt2[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
966  outputContainer->Add(fhMiInvPt2[index]) ;
967 
968  // Distance to bad module 3
969  snprintf(key, buffersize,"hMiInvPt_cen%d_pidbit%d_asy%d_dist3",ic,ipid,iasym) ;
970  snprintf(title, buffersize,"Mixed #it{M}_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f,dist bad 3",
971  ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
972  fhMiInvPt3[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
973  fhMiInvPt3[index]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
974  fhMiInvPt3[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
975  outputContainer->Add(fhMiInvPt3[index]) ;
976  }
977  }
978  }
979  }
980  }
981  }
982 
983  fhEPairDiffTime = new TH2F("hEPairDiffTime","cluster pair time difference vs #it{p}_{T}",nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
984  fhEPairDiffTime->SetXTitle("#it{p}_{T,pair} (GeV/#it{c})");
985  fhEPairDiffTime->SetYTitle("#Delta t (ns)");
986  outputContainer->Add(fhEPairDiffTime);
987 
989  {
990  fhRePtAsym = new TH2F("hRePtAsym","#it{Asymmetry} vs #it{p}_{T} , for pairs",nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax) ;
991  fhRePtAsym->SetXTitle("#it{p}_{T} (GeV/#it{c})");
992  fhRePtAsym->SetYTitle("#it{Asymmetry}");
993  outputContainer->Add(fhRePtAsym);
994 
995  fhRePtAsymPi0 = new TH2F("hRePtAsymPi0","#it{Asymmetry} vs #it{p}_{T} , for pairs close to #pi^{0} mass",nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax) ;
996  fhRePtAsymPi0->SetXTitle("#it{p}_{T} (GeV/#it{c})");
997  fhRePtAsymPi0->SetYTitle("Asymmetry");
998  outputContainer->Add(fhRePtAsymPi0);
999 
1000  fhRePtAsymEta = new TH2F("hRePtAsymEta","#it{Asymmetry} vs #it{p}_{T} , for pairs close to #eta mass",nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax) ;
1001  fhRePtAsymEta->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1002  fhRePtAsymEta->SetYTitle("#it{Asymmetry}");
1003  outputContainer->Add(fhRePtAsymEta);
1004 
1005  if(DoOwnMix())
1006  {
1007  fhMiPtAsym = new TH2F("hMiPtAsym","#it{Asymmetry} vs #it{p}_{T} , for mixed pairs",nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax) ;
1008  fhMiPtAsym->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1009  fhMiPtAsym->SetYTitle("#it{Asymmetry}");
1010  outputContainer->Add(fhMiPtAsym);
1011 
1012  fhMiPtAsymPi0 = new TH2F("hMiPtAsymPi0","#it{Asymmetry} vs #it{p}_{T} , for mixed pairs close to #pi^{0} mass",nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax) ;
1013  fhMiPtAsymPi0->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1014  fhMiPtAsymPi0->SetYTitle("Asymmetry");
1015  outputContainer->Add(fhMiPtAsymPi0);
1016 
1017  fhMiPtAsymEta = new TH2F("hMiPtAsymEta","#it{Asymmetry} vs #it{p}_{T} , for mixed pairs close to #eta mass",nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax) ;
1018  fhMiPtAsymEta->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1019  fhMiPtAsymEta->SetYTitle("#it{Asymmetry}");
1020  outputContainer->Add(fhMiPtAsymEta);
1021  }
1022  }
1023 
1024  if(fMultiCutAnaAcc)
1025  {
1026  for(Int_t ipt=0; ipt<fNPtCuts; ipt++)
1027  {
1028  fhPtBinClusterEtaPhi[ipt] = new TH2F
1029  (Form("hPtBin%d_Cluster_EtaPhi",ipt),
1030  Form("#eta vs #phi, %2.2f<#it{p}_{T}<%2.2f GeV/#it{c}",fPtCuts[ipt],fPtCuts[ipt+1]),
1031  netabins,etamin,etamax,nphibins,phimin,phimax);
1032  fhPtBinClusterEtaPhi[ipt]->SetYTitle("#phi (rad)");
1033  fhPtBinClusterEtaPhi[ipt]->SetXTitle("#eta");
1034  outputContainer->Add(fhPtBinClusterEtaPhi[ipt]) ;
1035 
1036  fhPtBinClusterColRow[ipt] = new TH2F
1037  (Form("hPtBin%d_Cluster_ColRow",ipt),
1038  Form("column vs row, %2.2f<#it{p}_{T}<%2.2f GeV/#it{c}",fPtCuts[ipt],fPtCuts[ipt+1]),
1039  96+2,-1.5,96+0.5,(8*24+2*8)+2,-1.5,(8*24+2*8)+0.5);
1040  fhPtBinClusterColRow[ipt]->SetYTitle("row");
1041  fhPtBinClusterColRow[ipt]->SetXTitle("column");
1042  outputContainer->Add(fhPtBinClusterColRow[ipt]) ;
1043  }
1044  }
1045 
1046  if(fMultiCutAna)
1047  {
1048 // fhRePIDBits = new TH2F*[fNPIDBits];
1049 // for(Int_t ipid=0; ipid<fNPIDBits; ipid++)
1050 // {
1051 // snprintf(key, buffersize,"hRe_pidbit%d",ipid) ;
1052 // snprintf(title, buffersize,"Real #it{M}_{#gamma#gamma} distr. for PIDBit=%d",fPIDBits[ipid]) ;
1053 // fhRePIDBits[ipid] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1054 // fhRePIDBits[ipid]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1055 // fhRePIDBits[ipid]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1056 // outputContainer->Add(fhRePIDBits[ipid]) ;
1057 // }// pid bit loop
1058 
1061 
1062  if(fFillAngleHisto)
1063  {
1066  }
1067 
1069  {
1070  for(Int_t iSM = 0; iSM < fNModules; iSM++)
1071  {
1074  }
1075  }
1076 
1077  for(Int_t ipt=0; ipt<fNPtCuts; ipt++)
1078  {
1079  for(Int_t icell=0; icell<fNCellNCuts; icell++)
1080  {
1081  for(Int_t iasym=0; iasym<fNAsymCuts; iasym++)
1082  {
1083  snprintf(key, buffersize,"hRe_pt%d_cell%d_asym%d",ipt,icell,iasym) ;
1084  snprintf(title, buffersize,"Real #it{M}_{#gamma#gamma} distr. for %1.1f< #it{p}_{T} < %1.1f, ncell>%d and asym<%1.2f ",
1085  fPtCuts[ipt],fPtCutsMax[ipt],fCellNCuts[icell], fAsymCuts[iasym]) ;
1086 
1087  Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
1088  //printf("ipt %d, icell %d, iassym %d, index %d\n",ipt, icell, iasym, index);
1089 
1090  fhRePtNCellAsymCuts[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1091  fhRePtNCellAsymCuts[index]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1092  fhRePtNCellAsymCuts[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1093  outputContainer->Add(fhRePtNCellAsymCuts[index]) ;
1094 
1095  if(DoOwnMix())
1096  {
1097  snprintf(key, buffersize,"hMi_pt%d_cell%d_asym%d",ipt,icell,iasym) ;
1098  snprintf(title, buffersize,"Mixed #it{M}_{#gamma#gamma} distr. for %1.1f< #it{p}_{T} < %1.1f, ncell>%d and asym<%1.2f",
1099  fPtCuts[ipt],fPtCutsMax[ipt],fCellNCuts[icell], fAsymCuts[iasym]) ;
1100  fhMiPtNCellAsymCuts[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1101  fhMiPtNCellAsymCuts[index]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1102  fhMiPtNCellAsymCuts[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1103  outputContainer->Add(fhMiPtNCellAsymCuts[index]) ;
1104  }
1105 
1107  {
1108  for(Int_t iSM = 0; iSM < fNModules; iSM++)
1109  {
1110  //printf("\t sm %d\n",iSM);
1111 
1112  snprintf(key, buffersize,"hRe_pt%d_cell%d_asym%d_SM%d",ipt,icell,iasym,iSM) ;
1113  snprintf(title, buffersize,"Real #it{M}_{#gamma#gamma} distr. for %1.1f< #it{p}_{T} < %1.1f, ncell>%d and asym<%1.2f, SM %d ",
1114  fPtCuts[ipt],fPtCutsMax[ipt],fCellNCuts[icell], fAsymCuts[iasym],iSM) ;
1115  fhRePtNCellAsymCutsSM[iSM][index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1116  fhRePtNCellAsymCutsSM[iSM][index]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1117  fhRePtNCellAsymCutsSM[iSM][index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1118  outputContainer->Add(fhRePtNCellAsymCutsSM[iSM][index]) ;
1119  }
1120 
1121  if(fFillAngleHisto)
1122  {
1123  snprintf(key, buffersize,"hReOpAngle_pt%d_cell%d_asym%d",ipt,icell,iasym) ;
1124  snprintf(title, buffersize,"Real #theta_{#gamma#gamma} distr. for %1.1f< #it{p}_{T} < %1.1f, ncell>%d and asym<%1.2f ",
1125  fPtCuts[ipt],fPtCutsMax[ipt],fCellNCuts[icell], fAsymCuts[iasym]) ;
1126 
1127  Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
1128  //printf("Angle ipt %d, icell %d, iassym %d, index %d, %s, %s\n",ipt, icell, iasym, index, key, title);
1129 
1130  fhRePtNCellAsymCutsOpAngle[index] = new TH2F(key,title,nptbins,ptmin,ptmax,280,0,1.4) ;
1131  fhRePtNCellAsymCutsOpAngle[index]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1132  fhRePtNCellAsymCutsOpAngle[index]->SetYTitle("#theta_{#gamma,#gamma} (rad)");
1133  outputContainer->Add(fhRePtNCellAsymCutsOpAngle[index]) ;
1134 
1135  if(DoOwnMix())
1136  {
1137  snprintf(key, buffersize,"hMiOpAngle_pt%d_cell%d_asym%d",ipt,icell,iasym) ;
1138  snprintf(title, buffersize,"Mixed #theta_{#gamma#gamma} distr. for %1.1f< #it{p}_{T} < %1.1f, ncell>%d and asym<%1.2f",
1139  fPtCuts[ipt],fPtCutsMax[ipt],fCellNCuts[icell], fAsymCuts[iasym]) ;
1140  fhMiPtNCellAsymCutsOpAngle[index] = new TH2F(key,title,nptbins,ptmin,ptmax,280,0,1.4) ;
1141  fhMiPtNCellAsymCutsOpAngle[index]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1142  fhMiPtNCellAsymCutsOpAngle[index]->SetYTitle("#theta_{#gamma,#gamma} (rad)");
1143  outputContainer->Add(fhMiPtNCellAsymCutsOpAngle[index]) ;
1144  }
1145 
1146  //printf("\t %p, %p\n", fhRePtNCellAsymCutsOpAngle[index], fhMiPtNCellAsymCutsOpAngle[index]);
1148  {
1149  for(Int_t iSM = 0; iSM < fNModules; iSM++)
1150  {
1151  snprintf(key, buffersize,"hReOpAngle_pt%d_cell%d_asym%d_SM%d",ipt,icell,iasym,iSM) ;
1152  snprintf(title, buffersize,"Real #it{M}_{#gamma#gamma} distr. for %1.1f< #it{p}_{T} < %1.1f, ncell>%d and asym<%1.2f, SM %d ",
1153  fPtCuts[ipt],fPtCutsMax[ipt],fCellNCuts[icell], fAsymCuts[iasym],iSM) ;
1154  fhRePtNCellAsymCutsSMOpAngle[iSM][index] = new TH2F(key,title,nptbins,ptmin,ptmax,280,0,1.4) ;
1155  fhRePtNCellAsymCutsSMOpAngle[iSM][index]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1156  fhRePtNCellAsymCutsSMOpAngle[iSM][index]->SetYTitle("#theta_{#gamma,#gamma} (rad)");
1157  outputContainer->Add(fhRePtNCellAsymCutsSMOpAngle[iSM][index]) ;
1158  }
1159  }
1160  }
1161  }
1162  }
1163  }
1164  }
1165 
1166 // if(ntrmbins!=0)
1167 // {
1168 // fhRePtMult = new TH3F*[fNAsymCuts] ;
1169 // for(Int_t iasym = 0; iasym<fNAsymCuts; iasym++)
1170 // {
1171 // fhRePtMult[iasym] = new TH3F(Form("hRePtMult_asym%d",iasym),Form("(#it{p}_{T},C,M)_{#gamma#gamma}, A<%1.2f",fAsymCuts[iasym]),
1172 // nptbins,ptmin,ptmax,ntrmbins,ntrmmin,ntrmmax,nmassbins,massmin,massmax);
1173 // fhRePtMult[iasym]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1174 // fhRePtMult[iasym]->SetYTitle("Track multiplicity");
1175 // fhRePtMult[iasym]->SetZTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1176 // outputContainer->Add(fhRePtMult[iasym]) ;
1177 // }
1178 // }
1179  } // multi cuts analysis
1180 
1182  {
1183  fhReSS[0] = new TH2F("hRe_SS_Tight"," 0.01 < #lambda_{0}^{2} < 0.4",
1184  nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
1185  fhReSS[0]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1186  fhReSS[0]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1187  outputContainer->Add(fhReSS[0]) ;
1188 
1189 
1190  fhReSS[1] = new TH2F("hRe_SS_Loose"," #lambda_{0}^{2} > 0.4",
1191  nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
1192  fhReSS[1]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1193  fhReSS[1]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1194  outputContainer->Add(fhReSS[1]) ;
1195 
1196 
1197  fhReSS[2] = new TH2F("hRe_SS_Both"," cluster_{1} #lambda_{0}^{2} > 0.4; cluster_{2} 0.01 < #lambda_{0}^{2} < 0.4",
1198  nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
1199  fhReSS[2]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1200  fhReSS[2]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1201  outputContainer->Add(fhReSS[2]) ;
1202  }
1203 
1204  if(DoOwnMix())
1205  {
1206  fhEventBin=new TH1I("hEventBin","Number of real pairs per bin(cen,vz,rp)",
1209  fhEventBin->SetXTitle("bin");
1210  outputContainer->Add(fhEventBin) ;
1211 
1212 
1213  fhEventMixBin=new TH1I("hEventMixBin","Number of mixed pairs per bin(cen,vz,rp)",
1216  fhEventMixBin->SetXTitle("bin");
1217  outputContainer->Add(fhEventMixBin) ;
1218  }
1219 
1221  {
1222  fhCentrality=new TH1F("hCentralityBin","Number of events in centrality bin",GetNCentrBin(),0.,1.*GetNCentrBin()) ;
1223  fhCentrality->SetXTitle("Centrality bin");
1224  outputContainer->Add(fhCentrality) ;
1225 
1226  fhCentralityNoPair=new TH1F("hCentralityBinNoPair","Number of events in centrality bin, with no cluster pairs",GetNCentrBin(),0.,1.*GetNCentrBin()) ;
1227  fhCentralityNoPair->SetXTitle("Centrality bin");
1228  outputContainer->Add(fhCentralityNoPair) ;
1229 
1230  fhEventPlaneResolution=new TH2F("hEventPlaneResolution","Event plane resolution",GetNCentrBin(),0,GetNCentrBin(),100,0.,TMath::TwoPi()) ;
1231  fhEventPlaneResolution->SetYTitle("Resolution");
1232  fhEventPlaneResolution->SetXTitle("Centrality Bin");
1233  outputContainer->Add(fhEventPlaneResolution) ;
1234  }
1235 
1236  if(fFillAngleHisto)
1237  {
1238  fhRealOpeningAngle = new TH2F
1239  ("hRealOpeningAngle","Angle between all #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,nopanbins,opanmin,opanmax);
1240  fhRealOpeningAngle->SetYTitle("#theta(rad)");
1241  fhRealOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
1242  outputContainer->Add(fhRealOpeningAngle) ;
1243 
1245  ("hRealCosOpeningAngle","Cosinus of angle between all #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,100,0,1);
1246  fhRealCosOpeningAngle->SetYTitle("cos (#theta) ");
1247  fhRealCosOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
1248  outputContainer->Add(fhRealCosOpeningAngle) ;
1249 
1250  if(DoOwnMix())
1251  {
1253  ("hMixedOpeningAngle","Angle between all #gamma pair vs E_{#pi^{0}}, Mixed pairs",nptbins,ptmin,ptmax,nopanbins,opanmin,opanmax);
1254  fhMixedOpeningAngle->SetYTitle("#theta(rad)");
1255  fhMixedOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
1256  outputContainer->Add(fhMixedOpeningAngle) ;
1257 
1259  ("hMixedCosOpeningAngle","Cosinus of angle between all #gamma pair vs E_{#pi^{0}}, Mixed pairs",nptbins,ptmin,ptmax,100,0,1);
1260  fhMixedCosOpeningAngle->SetYTitle("cos (#theta) ");
1261  fhMixedCosOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
1262  outputContainer->Add(fhMixedCosOpeningAngle) ;
1263  }
1264 
1266  {
1267  for(Int_t ism = 0; ism < 20; ism++)
1268  {
1269  fhRealOpeningAnglePerSM[ism] = new TH2F
1270  (Form("hRealOpeningAngleMod_%d",ism),
1271  Form("Angle between all #gamma pair vs E_{#pi^{0}}, SM %d",ism),
1272  nptbins,ptmin,ptmax,nopanbins,opanmin,opanmax);
1273  fhRealOpeningAnglePerSM[ism]->SetYTitle("#theta(rad)");
1274  fhRealOpeningAnglePerSM[ism]->SetXTitle("E_{ #pi^{0}} (GeV)");
1275  outputContainer->Add(fhRealOpeningAnglePerSM[ism]) ;
1276 
1277  if(DoOwnMix())
1278  {
1279  fhMixedOpeningAnglePerSM[ism] = new TH2F
1280  (Form("hMixedOpeningAngleMod_%d",ism),
1281  Form("Angle between all #gamma pair vs E_{#pi^{0}}, Mixed pairs, SM %d",ism),
1282  nptbins,ptmin,ptmax,nopanbins,opanmin,opanmax);
1283  fhMixedOpeningAnglePerSM[ism]->SetYTitle("#theta(rad)");
1284  fhMixedOpeningAnglePerSM[ism]->SetXTitle("E_{ #pi^{0}} (GeV)");
1285  outputContainer->Add(fhMixedOpeningAnglePerSM[ism]) ;
1286  }
1287  }
1288  }
1289  }
1290 
1291  // Histograms filled only if MC data is requested
1292  if(IsDataMC())
1293  {
1294  fhReMCFromConversion = new TH2F("hReMCFromConversion","Invariant mass of 2 clusters originated in conversions",
1295  nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
1296  fhReMCFromConversion->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1297  fhReMCFromConversion->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1298  outputContainer->Add(fhReMCFromConversion) ;
1299 
1300  fhReMCFromNotConversion = new TH2F("hReMCNotFromConversion","Invariant mass of 2 clusters not originated in conversions",
1301  nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
1302  fhReMCFromNotConversion->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1303  fhReMCFromNotConversion->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1304  outputContainer->Add(fhReMCFromNotConversion) ;
1305 
1306  fhReMCFromMixConversion = new TH2F("hReMCFromMixConversion","Invariant mass of 2 clusters one from conversion and the other not",
1307  nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
1308  fhReMCFromMixConversion->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1309  fhReMCFromMixConversion->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1310  outputContainer->Add(fhReMCFromMixConversion) ;
1311 
1312  if(fFillOriginHisto)
1313  {
1314  fhMCPi0PtOrigin = new TH2F("hMCPi0PtOrigin","Reconstructed pair from generated #pi^{0} #it{p}_{T} vs origin",nptbins,ptmin,ptmax,11,0,11) ;
1315  fhMCPi0PtOrigin->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1316  fhMCPi0PtOrigin->SetYTitle("Origin");
1317  fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(1 ,"Status 21");
1318  fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(2 ,"Quark");
1319  fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(3 ,"qq Resonances");
1320  fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(4 ,"Resonances");
1321  fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(5 ,"#rho");
1322  fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(6 ,"#omega");
1323  fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(7 ,"K");
1324  fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(8 ,"Other");
1325  fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(9 ,"#eta");
1326  fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(10 ,"#eta prime");
1327  outputContainer->Add(fhMCPi0PtOrigin) ;
1328 
1329  fhMCNotResonancePi0PtOrigin = new TH2F("hMCNotResonancePi0PtOrigin","Reconstructed pair from generated #pi^{0} #it{p}_{T} vs origin",nptbins,ptmin,ptmax,11,0,11) ;
1330  fhMCNotResonancePi0PtOrigin->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1331  fhMCNotResonancePi0PtOrigin->SetYTitle("Origin");
1332  fhMCNotResonancePi0PtOrigin->GetYaxis()->SetBinLabel(1 ,"Status 21");
1333  fhMCNotResonancePi0PtOrigin->GetYaxis()->SetBinLabel(2 ,"Quark");
1334  fhMCNotResonancePi0PtOrigin->GetYaxis()->SetBinLabel(3 ,"qq Resonances");
1335  fhMCNotResonancePi0PtOrigin->GetYaxis()->SetBinLabel(4 ,"Resonances");
1336  fhMCNotResonancePi0PtOrigin->GetYaxis()->SetBinLabel(5 ,"#rho");
1337  fhMCNotResonancePi0PtOrigin->GetYaxis()->SetBinLabel(6 ,"#omega");
1338  fhMCNotResonancePi0PtOrigin->GetYaxis()->SetBinLabel(7 ,"K");
1339  fhMCNotResonancePi0PtOrigin->GetYaxis()->SetBinLabel(8 ,"Other");
1340  fhMCNotResonancePi0PtOrigin->GetYaxis()->SetBinLabel(9 ,"#eta");
1341  fhMCNotResonancePi0PtOrigin->GetYaxis()->SetBinLabel(10 ,"#eta prime");
1342  outputContainer->Add(fhMCNotResonancePi0PtOrigin) ;
1343 
1344  fhMCPi0PtStatus = new TH2F("hMCPi0PtStatus","Reconstructed pair from generated #pi^{0} #it{p}_{T} vs status",nptbins,ptmin,ptmax,101,-50,50) ;
1345  fhMCPi0PtStatus->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1346  fhMCPi0PtStatus->SetYTitle("Status");
1347  outputContainer->Add(fhMCPi0PtStatus) ;
1348 
1349  // Eta
1350 
1351  fhMCEtaPtOrigin = new TH2F("hMCEtaPtOrigin","Reconstructed pair from generated #pi^{0} #it{p}_{T} vs origin",nptbins,ptmin,ptmax,7,0,7) ;
1352  fhMCEtaPtOrigin->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1353  fhMCEtaPtOrigin->SetYTitle("Origin");
1354  fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(1 ,"Status 21");
1355  fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(2 ,"Quark");
1356  fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(3 ,"qq Resonances");
1357  fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(4 ,"Resonances");
1358  fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(5 ,"Other");
1359  fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(6 ,"#eta prime");
1360  outputContainer->Add(fhMCEtaPtOrigin) ;
1361 
1362  fhMCPi0ProdVertex = new TH2F("hMCPi0ProdVertex","Selected reco pair from generated #pi^{0} #it{p}_{T} vs production vertex",
1363  200,0.,20.,5000,0,500) ;
1364  fhMCPi0ProdVertex->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1365  fhMCPi0ProdVertex->SetYTitle("#it{R} (cm)");
1366  outputContainer->Add(fhMCPi0ProdVertex) ;
1367 
1368  fhMCEtaProdVertex = new TH2F("hMCEtaProdVertex","Selected reco pair from generated #eta #it{p}_{T} vs production vertex",
1369  200,0.,20.,5000,0,500) ;
1370  fhMCEtaProdVertex->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1371  fhMCEtaProdVertex->SetYTitle("#it{R} (cm)");
1372  outputContainer->Add(fhMCEtaProdVertex) ;
1373 
1374  // Name of found ancestors of the cluster pairs. Check definitions in FillMCVsRecDataHistograms
1375  TString ancestorTitle[] = {"Photon","Electron","Pi0","Eta","AntiProton","AntiNeutron","Muon & converted stable particles",
1376  "Resonances","Strings","Initial state interaction","Final state radiations","Colliding protons","not found"};
1377 
1378  for(Int_t i = 0; i<13; i++)
1379  {
1380  fhMCOrgMass[i] = new TH2F(Form("hMCOrgMass_%d",i),Form("#it{M} vs #it{p}_{T}, ancestor %s",ancestorTitle[i].Data()),
1381  nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1382  fhMCOrgMass[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1383  fhMCOrgMass[i]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1384  outputContainer->Add(fhMCOrgMass[i]) ;
1385 
1386  fhMCOrgAsym[i]= new TH2F(Form("hMCOrgAsym_%d",i),Form("#it{Asymmetry} vs #it{p}_{T}, ancestor %s",ancestorTitle[i].Data()),
1387  nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax) ;
1388  fhMCOrgAsym[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1389  fhMCOrgAsym[i]->SetYTitle("A");
1390  outputContainer->Add(fhMCOrgAsym[i]) ;
1391 
1392  fhMCOrgDeltaEta[i] = new TH2F(Form("hMCOrgDeltaEta_%d",i),Form("#Delta #eta of pair vs #it{p}_{T}, ancestor %s",ancestorTitle[i].Data()),
1393  nptbins,ptmin,ptmax,netabins,-1.4,1.4) ;
1394  fhMCOrgDeltaEta[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1395  fhMCOrgDeltaEta[i]->SetYTitle("#Delta #eta");
1396  outputContainer->Add(fhMCOrgDeltaEta[i]) ;
1397 
1398  fhMCOrgDeltaPhi[i]= new TH2F(Form("hMCOrgDeltaPhi_%d",i),Form("#Delta #phi of pair vs #it{p}_{T}, ancestor %s",ancestorTitle[i].Data()),
1399  nptbins,ptmin,ptmax,nphibins,-0.7,0.7) ;
1400  fhMCOrgDeltaPhi[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1401  fhMCOrgDeltaPhi[i]->SetYTitle("#Delta #phi (rad)");
1402  outputContainer->Add(fhMCOrgDeltaPhi[i]) ;
1403  }
1404 
1405  if(fMultiCutAnaSim)
1406  {
1413 
1414  for(Int_t ipt=0; ipt<fNPtCuts; ipt++)
1415  {
1416  for(Int_t icell=0; icell<fNCellNCuts; icell++)
1417  {
1418  for(Int_t iasym=0; iasym<fNAsymCuts; iasym++)
1419  {
1420  Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
1421 
1422  fhMCPi0MassPtRec[index] = new TH2F(Form("hMCPi0MassPtRec_pt%d_cell%d_asym%d",ipt,icell,iasym),
1423  Form("Reconstructed #it{M} vs reconstructed #it{p}_{T} of true #pi^{0} cluster pairs for %1.1f<#it{p}_{T}<%1.1f, #it{N}^{cluster}_{cell}>%d and |#it{A}|<%1.2f",
1424  fPtCuts[ipt],fPtCutsMax[ipt],fCellNCuts[icell], fAsymCuts[iasym]),
1425  nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1426  fhMCPi0MassPtRec[index]->SetXTitle("#it{p}_{T, reconstructed} (GeV/#it{c})");
1427  fhMCPi0MassPtRec[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1428  outputContainer->Add(fhMCPi0MassPtRec[index]) ;
1429 
1430  fhMCPi0MassPtTrue[index] = new TH2F(Form("hMCPi0MassPtTrue_pt%d_cell%d_asym%d",ipt,icell,iasym),
1431  Form("Reconstructed #it{M} vs generated #it{p}_{T} of true #pi^{0} cluster pairs for %1.1f<#it{p}_{T}<%1.1f, #it{N}^{cluster}_{cell}>%d and |#it{A}|<%1.2f",
1432  fPtCuts[ipt],fPtCutsMax[ipt],fCellNCuts[icell], fAsymCuts[iasym]),
1433  nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1434  fhMCPi0MassPtTrue[index]->SetXTitle("#it{p}_{T, generated} (GeV/#it{c})");
1435  fhMCPi0MassPtTrue[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1436  outputContainer->Add(fhMCPi0MassPtTrue[index]) ;
1437 
1438  fhMCPi0PtTruePtRec[index] = new TH2F(Form("hMCPi0PtTruePtRec_pt%d_cell%d_asym%d",ipt,icell,iasym),
1439  Form("Generated vs reconstructed #it{p}_{T} of true #pi^{0} cluster pairs, 0.01 < rec. mass < 0.17 MeV/#it{c}^{2} for %1.1f<#it{p}_{T}<%1.1f, #it{N}^{cluster}_{cell}>%d and |#it{A}|>%1.2f",
1440  fPtCuts[ipt],fPtCutsMax[ipt],fCellNCuts[icell], fAsymCuts[iasym]),
1441  nptbins,ptmin,ptmax,nptbins,ptmin,ptmax) ;
1442  fhMCPi0PtTruePtRec[index]->SetXTitle("#it{p}_{T, generated} (GeV/#it{c})");
1443  fhMCPi0PtTruePtRec[index]->SetYTitle("#it{p}_{T, reconstructed} (GeV/#it{c})");
1444  outputContainer->Add(fhMCPi0PtTruePtRec[index]) ;
1445 
1446  fhMCEtaMassPtRec[index] = new TH2F(Form("hMCEtaMassPtRec_pt%d_cell%d_asym%d",ipt,icell,iasym),
1447  Form("Reconstructed #it{M} vs reconstructed #it{p}_{T} of true #eta cluster pairs for %1.1f<#it{p}_{T}<%1.1f, #it{N}^{cluster}_{cell}>%d and |#it{A}|<%1.2f",
1448  fPtCuts[ipt],fPtCutsMax[ipt],fCellNCuts[icell], fAsymCuts[iasym]),
1449  nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1450  fhMCEtaMassPtRec[index]->SetXTitle("#it{p}_{T, generated} (GeV/#it{c})");
1451  fhMCEtaMassPtRec[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1452  outputContainer->Add(fhMCEtaMassPtRec[index]) ;
1453 
1454  fhMCEtaMassPtTrue[index] = new TH2F(Form("hMCEtaMassPtTrue_pt%d_cell%d_asym%d",ipt,icell,iasym),
1455  Form("Reconstructed #it{M} vs generated #it{p}_{T} of true #eta cluster pairs for %1.1f<#it{p}_{T}<%1.1f, #it{N}^{cluster}_{cell}>%d and |#it{A}|<%1.2f",
1456  fPtCuts[ipt],fPtCutsMax[ipt],fCellNCuts[icell], fAsymCuts[iasym]),
1457  nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1458  fhMCEtaMassPtTrue[index]->SetXTitle("#it{p}_{T, generated} (GeV/#it{c})");
1459  fhMCEtaMassPtTrue[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1460  outputContainer->Add(fhMCEtaMassPtTrue[index]) ;
1461 
1462  fhMCEtaPtTruePtRec[index] = new TH2F(Form("hMCEtaPtTruePtRec_pt%d_cell%d_asym%d",ipt,icell,iasym),
1463  Form("Generated vs reconstructed #it{p}_{T} of true #eta cluster pairs, 0.01 < rec. mass < 0.17 MeV/#it{c}^{2} for %1.1f<#it{p}_{T}<%1.1f, ncell>%d and asym<%1.2f",
1464  fPtCuts[ipt],fPtCutsMax[ipt],fCellNCuts[icell], fAsymCuts[iasym]),
1465  nptbins,ptmin,ptmax,nptbins,ptmin,ptmax) ;
1466  fhMCEtaPtTruePtRec[index]->SetXTitle("#it{p}_{T, generated} (GeV/#it{c})");
1467  fhMCEtaPtTruePtRec[index]->SetYTitle("#it{p}_{T, reconstructed} (GeV/#it{c})");
1468  outputContainer->Add(fhMCEtaPtTruePtRec[index]) ;
1469  }
1470  }
1471  }
1472  }//multi cut ana
1473  else
1474  {
1475  fhMCPi0MassPtTrue = new TH2F*[1];
1476  fhMCPi0PtTruePtRec = new TH2F*[1];
1477  fhMCEtaMassPtTrue = new TH2F*[1];
1478  fhMCEtaPtTruePtRec = new TH2F*[1];
1479 
1480  fhMCPi0MassPtTrue[0] = new TH2F("hMCPi0MassPtTrue","Reconstructed Mass vs generated p_T of true #pi^{0} cluster pairs",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1481  fhMCPi0MassPtTrue[0]->SetXTitle("#it{p}_{T, generated} (GeV/#it{c})");
1482  fhMCPi0MassPtTrue[0]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1483  outputContainer->Add(fhMCPi0MassPtTrue[0]) ;
1484 
1485  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) ;
1486  fhMCPi0PtTruePtRec[0]->SetXTitle("#it{p}_{T, generated} (GeV/#it{c})");
1487  fhMCPi0PtTruePtRec[0]->SetYTitle("#it{p}_{T, reconstructed} (GeV/#it{c})");
1488  outputContainer->Add(fhMCPi0PtTruePtRec[0]) ;
1489 
1490  fhMCEtaMassPtTrue[0] = new TH2F("hMCEtaMassPtTrue","Reconstructed Mass vs generated p_T of true #eta cluster pairs",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1491  fhMCEtaMassPtTrue[0]->SetXTitle("#it{p}_{T, generated} (GeV/#it{c})");
1492  fhMCEtaMassPtTrue[0]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1493  outputContainer->Add(fhMCEtaMassPtTrue[0]) ;
1494 
1495  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) ;
1496  fhMCEtaPtTruePtRec[0]->SetXTitle("#it{p}_{T, generated} (GeV/#it{c})");
1497  fhMCEtaPtTruePtRec[0]->SetYTitle("#it{p}_{T, reconstructed} (GeV/#it{c})");
1498  outputContainer->Add(fhMCEtaPtTruePtRec[0]) ;
1499  }
1500  }
1501  }
1502 
1504  {
1506  {
1507  TString pairnamePHOS[] = {"(0-1)","(0-2)","(1-2)","(0-3)","(0-4)","(1-3)","(1-4)","(2-3)","(2-4)","(3-4)"};
1508  for(Int_t imod=0; imod<fNModules; imod++)
1509  {
1510  //Module dependent invariant mass
1511  snprintf(key, buffersize,"hReMod_%d",imod) ;
1512  snprintf(title, buffersize,"Real #it{M}_{#gamma#gamma} distr. for Module %d",imod) ;
1513  fhReMod[imod] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1514  fhReMod[imod]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1515  fhReMod[imod]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1516  outputContainer->Add(fhReMod[imod]) ;
1517  if(GetCalorimeter()==kPHOS)
1518  {
1519  snprintf(key, buffersize,"hReDiffPHOSMod_%d",imod) ;
1520  snprintf(title, buffersize,"Real pairs PHOS, clusters in different Modules: %s",(pairnamePHOS[imod]).Data()) ;
1521  fhReDiffPHOSMod[imod] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1522  fhReDiffPHOSMod[imod]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1523  fhReDiffPHOSMod[imod]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1524  outputContainer->Add(fhReDiffPHOSMod[imod]) ;
1525  }
1526  else
1527  {//EMCAL
1528  if(imod<fNModules/2)
1529  {
1530  snprintf(key, buffersize,"hReSameSectorEMCAL_%d",imod) ;
1531  snprintf(title, buffersize,"Real pairs EMCAL, clusters in same sector, SM(%d,%d)",imod*2,imod*2+1) ;
1532  fhReSameSectorEMCALMod[imod] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1533  fhReSameSectorEMCALMod[imod]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1534  fhReSameSectorEMCALMod[imod]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1535  outputContainer->Add(fhReSameSectorEMCALMod[imod]) ;
1536  }
1537  if(imod<fNModules-2)
1538  {
1539  snprintf(key, buffersize,"hReSameSideEMCAL_%d",imod) ;
1540  snprintf(title, buffersize,"Real pairs EMCAL, clusters in same side SM(%d,%d)",imod, imod+2) ;
1541  fhReSameSideEMCALMod[imod] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1542  fhReSameSideEMCALMod[imod]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1543  fhReSameSideEMCALMod[imod]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1544  outputContainer->Add(fhReSameSideEMCALMod[imod]) ;
1545  }
1546  }//EMCAL
1547 
1548  if(DoOwnMix())
1549  {
1550  snprintf(key, buffersize,"hMiMod_%d",imod) ;
1551  snprintf(title, buffersize,"Mixed #it{M}_{#gamma#gamma} distr. for Module %d",imod) ;
1552  fhMiMod[imod] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1553  fhMiMod[imod]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1554  fhMiMod[imod]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1555  outputContainer->Add(fhMiMod[imod]) ;
1556 
1557  if(GetCalorimeter()==kPHOS)
1558  {
1559  snprintf(key, buffersize,"hMiDiffPHOSMod_%d",imod) ;
1560  snprintf(title, buffersize,"Mixed pairs PHOS, clusters in different Modules: %s",(pairnamePHOS[imod]).Data()) ;
1561  fhMiDiffPHOSMod[imod] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1562  fhMiDiffPHOSMod[imod]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1563  fhMiDiffPHOSMod[imod]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1564  outputContainer->Add(fhMiDiffPHOSMod[imod]) ;
1565  }//PHOS
1566  else
1567  {//EMCAL
1568  if(imod<fNModules/2)
1569  {
1570  snprintf(key, buffersize,"hMiSameSectorEMCALMod_%d",imod) ;
1571  snprintf(title, buffersize,"Mixed pairs EMCAL, clusters in same sector, SM(%d,%d)",imod*2,imod*2+1) ;
1572  fhMiSameSectorEMCALMod[imod] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1573  fhMiSameSectorEMCALMod[imod]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1574  fhMiSameSectorEMCALMod[imod]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1575  outputContainer->Add(fhMiSameSectorEMCALMod[imod]) ;
1576  }
1577  if(imod<fNModules-2){
1578 
1579  snprintf(key, buffersize,"hMiSameSideEMCALMod_%d",imod) ;
1580  snprintf(title, buffersize,"Mixed pairs EMCAL, clusters in same side SM(%d,%d)",imod, imod+2) ;
1581  fhMiSameSideEMCALMod[imod] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1582  fhMiSameSideEMCALMod[imod]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1583  fhMiSameSideEMCALMod[imod]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1584  outputContainer->Add(fhMiSameSideEMCALMod[imod]) ;
1585  }
1586  } // EMCAL
1587  } // own mix
1588  } // loop combinations
1589  } // Not pair of detectors
1590  else
1591  {
1592  Int_t dcSameSM[6] = {12,13,14,15,16,17}; // Check eta order
1593  Int_t phSameSM[6] = {3, 3, 2, 2, 1, 1};
1594 
1595  Int_t dcDiffSM[8] = {12,13,14,15,16,17,0,0};
1596  Int_t phDiffSM[8] = {2, 2, 1, 1, 3, 3,0,0};
1597 
1598  for(Int_t icombi = 0; icombi < 8; icombi++)
1599  {
1600  snprintf(key, buffersize,"hReDiffSectorDCALPHOS_%d",icombi) ;
1601  snprintf(title, buffersize,"Real pairs DCAL-PHOS, clusters in different sector, SM(%d,%d)",dcDiffSM[icombi],phDiffSM[icombi]) ;
1602  fhReDiffSectorDCALPHOSMod[icombi] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1603  fhReDiffSectorDCALPHOSMod[icombi]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1604  fhReDiffSectorDCALPHOSMod[icombi]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1605  outputContainer->Add(fhReDiffSectorDCALPHOSMod[icombi]) ;
1606  if(DoOwnMix())
1607  {
1608  snprintf(key, buffersize,"hMiDiffSectorDCALPHOS_%d",icombi) ;
1609  snprintf(title, buffersize,"Mixed pairs DCAL-PHOS, clusters in different sector, SM(%d,%d)",dcDiffSM[icombi],phDiffSM[icombi]) ;
1610  fhMiDiffSectorDCALPHOSMod[icombi] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1611  fhMiDiffSectorDCALPHOSMod[icombi]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1612  fhMiDiffSectorDCALPHOSMod[icombi]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1613  outputContainer->Add(fhMiDiffSectorDCALPHOSMod[icombi]) ;
1614  }
1615 
1616  if(icombi > 5) continue ;
1617 
1618  snprintf(key, buffersize,"hReSameSectorDCALPHOS_%d",icombi) ;
1619  snprintf(title, buffersize,"Real pairs DCAL-PHOS, clusters in same sector, SM(%d,%d)",dcSameSM[icombi],phSameSM[icombi]) ;
1620  fhReSameSectorDCALPHOSMod[icombi] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1621  fhReSameSectorDCALPHOSMod[icombi]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1622  fhReSameSectorDCALPHOSMod[icombi]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1623  outputContainer->Add(fhReSameSectorDCALPHOSMod[icombi]) ;
1624  if(DoOwnMix())
1625  {
1626  snprintf(key, buffersize,"hMiSameSectorDCALPHOS_%d",icombi) ;
1627  snprintf(title, buffersize,"Mixed pairs DCAL-PHOS, clusters in same sector, SM(%d,%d)",dcSameSM[icombi],phSameSM[icombi]) ;
1628  fhMiSameSectorDCALPHOSMod[icombi] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1629  fhMiSameSectorDCALPHOSMod[icombi]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1630  fhMiSameSectorDCALPHOSMod[icombi]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1631  outputContainer->Add(fhMiSameSectorDCALPHOSMod[icombi]) ;
1632  }
1633  }
1634 
1635  }
1636  } // SM combinations
1637 
1638  //
1640  {
1641  for(Int_t icut = 0; icut < fNAngleCutBins; icut++)
1642  {
1644  (Form("hReOpAngleBin%d_ClusterMin_EtaPhi",icut),
1645  Form("#eta vs #phi, cluster pair lower #it{E}, pair %1.4f<#theta<%1.4f",fAngleCutBinsArray[icut],fAngleCutBinsArray[icut+1]),
1646  netabins,etamin,etamax,nphibins,phimin,phimax);
1647  fhReOpAngleBinMinClusterEtaPhi[icut]->SetYTitle("#phi (rad)");
1648  fhReOpAngleBinMinClusterEtaPhi[icut]->SetXTitle("#eta");
1649  outputContainer->Add(fhReOpAngleBinMinClusterEtaPhi[icut]) ;
1650 
1652  (Form("hReOpAngleBin%d_ClusterMax_EtaPhi",icut),
1653  Form("#eta vs #phi, cluster pair higher #it{E}, pair %1.4f<#theta<%1.4f",fAngleCutBinsArray[icut],fAngleCutBinsArray[icut+1]),
1654  netabins,etamin,etamax,nphibins,phimin,phimax);
1655  fhReOpAngleBinMaxClusterEtaPhi[icut]->SetYTitle("#phi (rad)");
1656  fhReOpAngleBinMaxClusterEtaPhi[icut]->SetXTitle("#eta");
1657  outputContainer->Add(fhReOpAngleBinMaxClusterEtaPhi[icut]) ;
1658 
1660  (Form("hReOpAngleBin%d_ClusterMin_ColRow",icut),
1661  Form("highest #it{E} cell, column vs row, cluster pair lower #it{E}, pair %1.4f<#theta<%1.4f",fAngleCutBinsArray[icut],fAngleCutBinsArray[icut+1]),
1662  96+2,-1.5,96+0.5,(8*24+2*8)+2,-1.5,(8*24+2*8)+0.5);
1663  fhReOpAngleBinMinClusterColRow[icut]->SetYTitle("row");
1664  fhReOpAngleBinMinClusterColRow[icut]->SetXTitle("column");
1665  outputContainer->Add(fhReOpAngleBinMinClusterColRow[icut]) ;
1666 
1668  (Form("hReOpAngleBin%d_ClusterMax_ColRow",icut),
1669  Form("highest #it{E} cell, column vs row, cluster pair higher #it{E}, pair %1.4f<#theta<%1.4f",fAngleCutBinsArray[icut],fAngleCutBinsArray[icut+1]),
1670  96+2,-1.5,96+0.5,(8*24+2*8)+2,-1.5,(8*24+2*8)+0.5);
1671  fhReOpAngleBinMaxClusterColRow[icut]->SetYTitle("row");
1672  fhReOpAngleBinMaxClusterColRow[icut]->SetXTitle("column");
1673  outputContainer->Add(fhReOpAngleBinMaxClusterColRow[icut]) ;
1674 
1676  (Form("hReOpAngleBin%d_ClusterMin_EPerSM",icut),
1677  Form("cluster pair lower #it{E}, pair %1.4f<#theta<%1.4f",fAngleCutBinsArray[icut],fAngleCutBinsArray[icut+1]),
1678  nptbins,ptmin,ptmax,20,0,20);
1679  fhReOpAngleBinMinClusterEPerSM[icut]->SetYTitle("SM");
1680  fhReOpAngleBinMinClusterEPerSM[icut]->SetXTitle("#it{E} (GeV/#it{c})");
1681  outputContainer->Add(fhReOpAngleBinMinClusterEPerSM[icut]) ;
1682 
1684  (Form("hReOpAngleBin%d_ClusterMax_EPerSM",icut),
1685  Form("cluster pair higher #it{E}, pair %1.4f<#theta<%1.4f",fAngleCutBinsArray[icut],fAngleCutBinsArray[icut+1]),
1686  nptbins,ptmin,ptmax,20,0,20);
1687  fhReOpAngleBinMaxClusterEPerSM[icut]->SetYTitle("SM");
1688  fhReOpAngleBinMaxClusterEPerSM[icut]->SetXTitle("#it{E} (GeV/#it{c})");
1689  outputContainer->Add(fhReOpAngleBinMaxClusterEPerSM[icut]) ;
1690 
1692  (Form("hReOpAngleBin%d_ClusterMin_TimePerSM",icut),
1693  Form("cluster pair lower #it{E}, pair %1.4f<#theta<%1.4f",fAngleCutBinsArray[icut],fAngleCutBinsArray[icut+1]),
1694  ntimebins,timemin,timemax,20,0,20);
1695  fhReOpAngleBinMinClusterTimePerSM[icut]->SetYTitle("SM");
1696  fhReOpAngleBinMinClusterTimePerSM[icut]->SetXTitle("#it{t} (ns)");
1697  outputContainer->Add(fhReOpAngleBinMinClusterTimePerSM[icut]) ;
1698 
1700  (Form("hReOpAngleBin%d_ClusterMax_TimePerSM",icut),
1701  Form("cluster pair higher #it{E}, pair %1.4f<#theta<%1.4f",fAngleCutBinsArray[icut],fAngleCutBinsArray[icut+1]),
1702  ntimebins,timemin,timemax,20,0,20);
1703  fhReOpAngleBinMaxClusterTimePerSM[icut]->SetYTitle("SM");
1704  fhReOpAngleBinMaxClusterTimePerSM[icut]->SetXTitle("#it{t} (ns)");
1705  outputContainer->Add(fhReOpAngleBinMaxClusterTimePerSM[icut]) ;
1706 
1708  (Form("hReOpAngleBin%d_ClusterMin_NCellPerSM",icut),
1709  Form("cluster pair lower #it{E}, pair %1.4f<#theta<%1.4f",fAngleCutBinsArray[icut],fAngleCutBinsArray[icut+1]),
1710  30,0,30,20,0,20);
1711  fhReOpAngleBinMinClusterNCellPerSM[icut]->SetYTitle("SM");
1712  fhReOpAngleBinMinClusterNCellPerSM[icut]->SetXTitle("# cells");
1713  outputContainer->Add(fhReOpAngleBinMinClusterNCellPerSM[icut]) ;
1714 
1716  (Form("hReOpAngleBin%d_ClusterMax_NCellPerSM",icut),
1717  Form("cluster pair higher #it{E}, pair %1.4f<#theta<%1.4f",fAngleCutBinsArray[icut],fAngleCutBinsArray[icut+1]),
1718  30,0,30,20,0,20);
1719  fhReOpAngleBinMaxClusterNCellPerSM[icut]->SetYTitle("SM");
1720  fhReOpAngleBinMaxClusterNCellPerSM[icut]->SetXTitle("# cells");
1721  outputContainer->Add(fhReOpAngleBinMaxClusterNCellPerSM[icut]) ;
1722 
1724  (Form("hReOpAngleBin%d_PairCluster_RatioPerSM",icut),
1725  Form("cluster pair #it{E}_{high}/ #it{E}_{low}, pair %1.4f<#theta<%1.4f",fAngleCutBinsArray[icut],fAngleCutBinsArray[icut+1]),
1726  100,0,1,20,0,20);
1727  fhReOpAngleBinPairClusterRatioPerSM[icut]->SetYTitle("SM");
1728  fhReOpAngleBinPairClusterRatioPerSM[icut]->SetXTitle("#it{E}_{low}/ #it{E}_{high}");
1729  outputContainer->Add(fhReOpAngleBinPairClusterRatioPerSM[icut]) ;
1730 
1732  (Form("hReOpAngleBin%d_PairCluster_MassPerSM",icut),
1733  Form("cluster pair #it{M}, pair %1.4f<#theta<%1.4f",fAngleCutBinsArray[icut],fAngleCutBinsArray[icut+1]),
1734  nmassbins,massmin,massmax,20,0,20);
1735  fhReOpAngleBinPairClusterMassPerSM[icut]->SetXTitle("#it{M} (GeV/#it{c}^2)");
1736  fhReOpAngleBinPairClusterMassPerSM[icut]->SetYTitle("SM");
1737  outputContainer->Add(fhReOpAngleBinPairClusterMassPerSM[icut]) ;
1738 
1740  (Form("hReOpAngleBin%d_PairCluster_Mass",icut),
1741  Form("cluster pair #it{M}, pair %1.4f<#theta<%1.4f",fAngleCutBinsArray[icut],fAngleCutBinsArray[icut+1]),
1742  nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
1743  fhReOpAngleBinPairClusterMass[icut]->SetYTitle("#it{M} (GeV/#it{c}^2)");
1744  fhReOpAngleBinPairClusterMass[icut]->SetXTitle("#it{p}_{T} GeV/#it{c}");
1745  outputContainer->Add(fhReOpAngleBinPairClusterMass[icut]) ;
1746 
1747  if(IsDataMC())
1748  {
1749  if(fFillOriginHisto)
1750  {
1752  (Form("hReOpAngleBin%d_PairCluster_MassMCTruePi0",icut),
1753  Form("cluster pair #it{M}, pair %1.4f<#theta<%1.4f, true mc #pi^{0}",fAngleCutBinsArray[icut],fAngleCutBinsArray[icut+1]),
1754  nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
1755  fhReOpAngleBinPairClusterMassMCTruePi0[icut]->SetYTitle("#it{M} (GeV/#it{c}^2)");
1756  fhReOpAngleBinPairClusterMassMCTruePi0[icut]->SetXTitle("#it{p}_{T} GeV/#it{c}");
1757  outputContainer->Add(fhReOpAngleBinPairClusterMassMCTruePi0[icut]) ;
1758 
1760  (Form("hReOpAngleBin%d_PairCluster_MassMCTrueEta",icut),
1761  Form("cluster pair #it{M}, pair %1.4f<#theta<%1.4f, true mc #eta",fAngleCutBinsArray[icut],fAngleCutBinsArray[icut+1]),
1762  nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
1763  fhReOpAngleBinPairClusterMassMCTrueEta[icut]->SetYTitle("#it{M} (GeV/#it{c}^2)");
1764  fhReOpAngleBinPairClusterMassMCTrueEta[icut]->SetXTitle("#it{p}_{T} GeV/#it{c}");
1765  outputContainer->Add(fhReOpAngleBinPairClusterMassMCTrueEta[icut]) ;
1766  }
1767 
1768  fhPrimPi0AccPtOpAngCuts[icut] = new TH1F
1769  (Form("hPrimPi0AccPt_OpAngleBin%d",icut),
1770  Form("Primary #pi^{0} #it{p}_{T} with both photons in acceptance, pair %1.4f<#theta<%1.4f",fAngleCutBinsArray[icut],fAngleCutBinsArray[icut+1]),
1771  nptbins,ptmin,ptmax) ;
1772  fhPrimPi0AccPtOpAngCuts[icut]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1773  outputContainer->Add(fhPrimPi0AccPtOpAngCuts[icut]) ;
1774 
1775  fhPrimEtaAccPtOpAngCuts[icut] = new TH1F
1776  (Form("hPrimEtaAccPt_OpAngleBin%d",icut),
1777  Form("Primary #eta #it{p}_{T} with both photons in acceptance, pair %1.4f<#theta<%1.4f",fAngleCutBinsArray[icut],fAngleCutBinsArray[icut+1]),
1778  nptbins,ptmin,ptmax) ;
1779  fhPrimEtaAccPtOpAngCuts[icut]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1780  outputContainer->Add(fhPrimEtaAccPtOpAngCuts[icut]) ;
1781  }
1782 
1783 // fhReOpAngleBinPairClusterAbsIdMaxCell[icut] = new TH2F
1784 // (Form("hReOpAngleBin%d_PairCluster_AbsIdCell",icut),
1785 // Form("cluster pair Abs Cell ID low #it{E} vs high #it{E}, pair %1.4f<#theta<%1.4f",fAngleCutBinsArray[icut],fAngleCutBinsArray[icut+1]),
1786 // //17664,0,17664,17664,0,17664);
1787 // 1689,0,16896,1689,0,16896);
1788 // fhReOpAngleBinPairClusterAbsIdMaxCell[icut]->SetYTitle("AbsId-higher");
1789 // fhReOpAngleBinPairClusterAbsIdMaxCell[icut]->SetXTitle("AbsId-lower");
1790 // outputContainer->Add(fhReOpAngleBinPairClusterAbsIdMaxCell[icut]) ;
1791 
1792  if( DoOwnMix() )
1793  {
1795  (Form("hMiOpAngleBin%d_ClusterMin_EtaPhi",icut),
1796  Form("#eta vs #phi, mixed cluster pair lower #it{E}, pair %1.4f<#theta<%1.4f",fAngleCutBinsArray[icut],fAngleCutBinsArray[icut+1]),
1797  netabins,etamin,etamax,nphibins,phimin,phimax);
1798  fhMiOpAngleBinMinClusterEtaPhi[icut]->SetYTitle("#phi (rad)");
1799  fhMiOpAngleBinMinClusterEtaPhi[icut]->SetXTitle("#eta");
1800  outputContainer->Add(fhMiOpAngleBinMinClusterEtaPhi[icut]) ;
1801 
1803  (Form("hMiOpAngleBin%d_ClusterMax_EtaPhi",icut),
1804  Form("#eta vs #phi, mixed cluster pair higher #it{E}, pair %1.4f<#theta<%1.4f",fAngleCutBinsArray[icut],fAngleCutBinsArray[icut+1]),
1805  netabins,etamin,etamax,nphibins,phimin,phimax);
1806  fhMiOpAngleBinMaxClusterEtaPhi[icut]->SetYTitle("#phi (rad)");
1807  fhMiOpAngleBinMaxClusterEtaPhi[icut]->SetXTitle("#eta");
1808  outputContainer->Add(fhMiOpAngleBinMaxClusterEtaPhi[icut]) ;
1809 
1810 // fhMiOpAngleBinMinClusterColRow[icut] = new TH2F
1811 // (Form("hMiOpAngleBin%d_ClusterMin_ColRow",icut),
1812 // Form("highest #it{E} cell, column vs row, mixed cluster pair lower #it{E}, pair %1.4f<#theta<%1.4f",fAngleCutBinsArray[icut],fAngleCutBinsArray[icut+1]),
1813 // 96+2,-1.5,96+0.5,(8*24+2*8)+2,-1.5,(8*24+2*8)+0.5);
1814 // fhMiOpAngleBinMinClusterColRow[icut]->SetYTitle("row");
1815 // fhMiOpAngleBinMinClusterColRow[icut]->SetXTitle("column");
1816 // outputContainer->Add(fhMiOpAngleBinMinClusterColRow[icut]) ;
1817 //
1818 // fhMiOpAngleBinMaxClusterColRow[icut] = new TH2F
1819 // (Form("hMiOpAngleBin%d_ClusterMax_ColRow",icut),
1820 // Form("highest #it{E} cell, column vs row, mixed cluster pair higher #it{E}, pair %1.4f<#theta<%1.4f",fAngleCutBinsArray[icut],fAngleCutBinsArray[icut+1]),
1821 // 96+2,-1.5,96+0.5,(8*24+2*8)+2,-1.5,(8*24+2*8)+0.5);
1822 // fhMiOpAngleBinMaxClusterColRow[icut]->SetYTitle("row");
1823 // fhMiOpAngleBinMaxClusterColRow[icut]->SetXTitle("column");
1824 // outputContainer->Add(fhMiOpAngleBinMaxClusterColRow[icut]) ;
1825 
1827  (Form("hMiOpAngleBin%d_ClusterMin_EPerSM",icut),
1828  Form("mixed cluster pair lower #it{E}, pair %1.4f<#theta<%1.4f",fAngleCutBinsArray[icut],fAngleCutBinsArray[icut+1]),
1829  nptbins,ptmin,ptmax,20,0,20);
1830  fhMiOpAngleBinMinClusterEPerSM[icut]->SetYTitle("SM");
1831  fhMiOpAngleBinMinClusterEPerSM[icut]->SetXTitle("#it{E} (GeV/#it{c})");
1832  outputContainer->Add(fhMiOpAngleBinMinClusterEPerSM[icut]) ;
1833 
1835  (Form("hMiOpAngleBin%d_ClusterMax_EPerSM",icut),
1836  Form("mixed cluster pair higher #it{E}, pair %1.4f<#theta<%1.4f",fAngleCutBinsArray[icut],fAngleCutBinsArray[icut+1]),
1837  nptbins,ptmin,ptmax,20,0,20);
1838  fhMiOpAngleBinMaxClusterEPerSM[icut]->SetYTitle("SM");
1839  fhMiOpAngleBinMaxClusterEPerSM[icut]->SetXTitle("#it{E} (GeV/#it{c})");
1840  outputContainer->Add(fhMiOpAngleBinMaxClusterEPerSM[icut]) ;
1841 
1843  (Form("hMiOpAngleBin%d_ClusterMin_TimePerSM",icut),
1844  Form("mixed cluster pair lower #it{E}, pair %1.4f<#theta<%1.4f",fAngleCutBinsArray[icut],fAngleCutBinsArray[icut+1]),
1845  ntimebins,timemin,timemax,20,0,20);
1846  fhMiOpAngleBinMinClusterTimePerSM[icut]->SetYTitle("SM");
1847  fhMiOpAngleBinMinClusterTimePerSM[icut]->SetXTitle("#it{t} (ns)");
1848  outputContainer->Add(fhMiOpAngleBinMinClusterTimePerSM[icut]) ;
1849 
1851  (Form("hMiOpAngleBin%d_ClusterMax_TimePerSM",icut),
1852  Form("mixed cluster pair higher #it{E}, pair %1.4f<#theta<%1.4f",fAngleCutBinsArray[icut],fAngleCutBinsArray[icut+1]),
1853  ntimebins,timemin,timemax,20,0,20);
1854  fhMiOpAngleBinMaxClusterTimePerSM[icut]->SetYTitle("SM");
1855  fhMiOpAngleBinMaxClusterTimePerSM[icut]->SetXTitle("#it{t} (ns)");
1856  outputContainer->Add(fhMiOpAngleBinMaxClusterTimePerSM[icut]) ;
1857 
1859  (Form("hMiOpAngleBin%d_ClusterMin_NCellPerSM",icut),
1860  Form("mixed cluster pair lower #it{E}, pair %1.4f<#theta<%1.4f",fAngleCutBinsArray[icut],fAngleCutBinsArray[icut+1]),
1861  30,0,30,20,0,20);
1862  fhMiOpAngleBinMinClusterNCellPerSM[icut]->SetYTitle("SM");
1863  fhMiOpAngleBinMinClusterNCellPerSM[icut]->SetXTitle("# cells");
1864  outputContainer->Add(fhMiOpAngleBinMinClusterNCellPerSM[icut]) ;
1865 
1867  (Form("hMiOpAngleBin%d_ClusterMax_NCellPerSM",icut),
1868  Form("mixed cluster pair higher #it{E}, pair %1.4f<#theta<%1.4f",fAngleCutBinsArray[icut],fAngleCutBinsArray[icut+1]),
1869  30,0,30,20,0,20);
1870  fhMiOpAngleBinMaxClusterNCellPerSM[icut]->SetYTitle("SM");
1871  fhMiOpAngleBinMaxClusterNCellPerSM[icut]->SetXTitle("# cells");
1872  outputContainer->Add(fhMiOpAngleBinMaxClusterNCellPerSM[icut]) ;
1873 
1875  (Form("hMiOpAngleBin%d_PairCluster_RatioPerSM",icut),
1876  Form("mixed cluster pair #it{E}_{high}/ #it{E}_{low}, pair %1.4f<#theta<%1.4f",fAngleCutBinsArray[icut],fAngleCutBinsArray[icut+1]),
1877  100,0,1,20,0,20);
1878  fhMiOpAngleBinPairClusterRatioPerSM[icut]->SetYTitle("SM");
1879  fhMiOpAngleBinPairClusterRatioPerSM[icut]->SetXTitle("#it{E}_{low}/ #it{E}_{high}");
1880  outputContainer->Add(fhMiOpAngleBinPairClusterRatioPerSM[icut]) ;
1881 
1883  (Form("hMiOpAngleBin%d_PairCluster_MassPerSM",icut),
1884  Form("cluster pair #it{M}, pair %1.4f<#theta<%1.4f",fAngleCutBinsArray[icut],fAngleCutBinsArray[icut+1]),
1885  nmassbins,massmin,massmax,20,0,20);
1886  fhMiOpAngleBinPairClusterMassPerSM[icut]->SetXTitle("#it{M} (GeV/#it{c}^2)");
1887  fhMiOpAngleBinPairClusterMassPerSM[icut]->SetYTitle("SM");
1888  outputContainer->Add(fhMiOpAngleBinPairClusterMassPerSM[icut]) ;
1889 
1891  (Form("hMiOpAngleBin%d_PairCluster_Mass",icut),
1892  Form("cluster pair #it{M}, pair %1.4f<#theta<%1.4f",fAngleCutBinsArray[icut],fAngleCutBinsArray[icut+1]),
1893  nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
1894  fhMiOpAngleBinPairClusterMass[icut]->SetYTitle("#it{M} (GeV/#it{c}^2)");
1895  fhMiOpAngleBinPairClusterMass[icut]->SetXTitle("#it{p}_{T} GeV/#it{c}");
1896  outputContainer->Add(fhMiOpAngleBinPairClusterMass[icut]) ;
1897 
1898 // fhMiOpAngleBinPairClusterAbsIdMaxCell[icut] = new TH2F
1899 // (Form("hMiOpAngleBin%d_PairCluster_AbsIdCell",icut),
1900 // Form("mixed cluster pair Abs Cell ID low #it{E} vs high #it{E}, pair %1.4f<#theta<%1.4f",fAngleCutBinsArray[icut],fAngleCutBinsArray[icut+1]),
1901 // ntimebins,timemin,timemax,20,0,20);
1902 // fhMiOpAngleBinPairClusterRatioPerSM[icut]->SetYTitle("AbsId-higher");
1903 // fhMiOpAngleBinPairClusterRatioPerSM[icut]->SetXTitle("AbsId-lower");
1904 // outputContainer->Add(fhMiOpAngleBinPairClusterRatioPerSM[icut]) ;
1905  }
1906  } // angle bin
1907  }
1908 
1909 
1910 // for(Int_t i = 0; i < outputContainer->GetEntries() ; i++){
1911 //
1912 // printf("Histogram %d, name: %s\n ",i, outputContainer->At(i)->GetName());
1913 //
1914 // }
1915 
1916  return outputContainer;
1917 }
1918 
1919 //___________________________________________________
1921 //___________________________________________________
1922 void AliAnaPi0::Print(const Option_t * /*opt*/) const
1923 {
1924  printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
1926 
1927  printf("Number of bins in Centrality: %d \n",GetNCentrBin()) ;
1928  printf("Number of bins in Z vert. pos: %d \n",GetNZvertBin()) ;
1929  printf("Number of bins in Reac. Plain: %d \n",GetNRPBin()) ;
1930  printf("Depth of event buffer: %d \n",GetNMaxEvMix()) ;
1931  printf("Pair in same Module: %d \n",fSameSM) ;
1932  printf("Cuts: \n") ;
1933  // printf("Z vertex position: -%2.3f < z < %2.3f \n",GetZvertexCut(),GetZvertexCut()) ; //It crashes here, why?
1934  printf("Number of modules: %d \n",fNModules) ;
1935  printf("Select pairs with their angle: %d, edep %d, min angle %2.3f, max angle %2.3f \n",fUseAngleCut, fUseAngleEDepCut, fAngleCut, fAngleMaxCut) ;
1936  printf("Asymmetry cuts: n = %d, \n",fNAsymCuts) ;
1937  printf("\tasymmetry < ");
1938  for(Int_t i = 0; i < fNAsymCuts; i++) printf("%2.2f ",fAsymCuts[i]);
1939  printf("\n");
1940 
1941  printf("PID selection bits: n = %d, \n",fNPIDBits) ;
1942  printf("\tPID bit = ");
1943  for(Int_t i = 0; i < fNPIDBits; i++) printf("%d ",fPIDBits[i]);
1944  printf("\n");
1945 
1947  {
1948  printf("pT cuts: n = %d, \n",fNPtCuts) ;
1949  printf("\tpT > ");
1950  for(Int_t i = 0; i < fNPtCuts; i++) printf("%2.2f ",fPtCuts[i]);
1951  printf("GeV/c\n");
1952  printf("\tpT < ");
1953  for(Int_t i = 0; i < fNPtCuts; i++) printf("%2.2f ",fPtCutsMax[i]);
1954  printf("GeV/c\n");
1955 
1956  printf("N cell in cluster cuts: n = %d, \n",fNCellNCuts) ;
1957  printf("\tnCell > ");
1958  for(Int_t i = 0; i < fNCellNCuts; i++) printf("%d ",fCellNCuts[i]);
1959  printf("\n");
1960 
1961  }
1962  printf("------------------------------------------------------\n") ;
1963 }
1964 
1965 //________________________________________
1967 //________________________________________
1969 {
1970  Double_t mesonY = -100 ;
1971  Double_t mesonE = -1 ;
1972  Double_t mesonPt = -1 ;
1973  Double_t mesonPhi = 100 ;
1974  Double_t mesonYeta= -1 ;
1975 
1976  Int_t pdg = 0 ;
1977  Int_t nprim = 0 ;
1978  Int_t nDaught = 0 ;
1979  Int_t iphot1 = 0 ;
1980  Int_t iphot2 = 0 ;
1981 
1982  Float_t cen = GetEventCentrality();
1983  Float_t ep = GetEventPlaneAngle();
1984 
1985  TParticle * primStack = 0;
1986  AliAODMCParticle * primAOD = 0;
1987 
1988  // Get the ESD MC particles container
1989  AliStack * stack = 0;
1990  if( GetReader()->ReadStack() )
1991  {
1992  stack = GetMCStack();
1993  if(!stack ) return;
1994  nprim = stack->GetNtrack();
1995  }
1996 
1997  // Get the AOD MC particles container
1998  TClonesArray * mcparticles = 0;
1999  if( GetReader()->ReadAODMCParticles() )
2000  {
2001  mcparticles = GetReader()->GetAODMCParticles();
2002  if( !mcparticles ) return;
2003  nprim = mcparticles->GetEntriesFast();
2004  }
2005 
2006  for(Int_t i=0 ; i < nprim; i++)
2007  {
2008  if ( !GetReader()->AcceptParticleMCLabel( i ) ) continue ;
2009 
2010  if(GetReader()->ReadStack())
2011  {
2012  primStack = stack->Particle(i) ;
2013  if(!primStack)
2014  {
2015  AliWarning("ESD primaries pointer not available!!");
2016  continue;
2017  }
2018 
2019  // If too small skip
2020  if( primStack->Energy() < 0.4 ) continue;
2021 
2022  pdg = primStack->GetPdgCode();
2023  nDaught = primStack->GetNDaughters();
2024  iphot1 = primStack->GetDaughter(0) ;
2025  iphot2 = primStack->GetDaughter(1) ;
2026 
2027  // Protection against floating point exception
2028  if ( primStack->Energy() == TMath::Abs(primStack->Pz()) ||
2029  (primStack->Energy() - primStack->Pz()) < 1e-3 ||
2030  (primStack->Energy() + primStack->Pz()) < 0 ) continue ;
2031 
2032  //printf("i %d, %s %d %s %d \n",i, stack->Particle(i)->GetName(), stack->Particle(i)->GetPdgCode(),
2033  // prim->GetName(), prim->GetPdgCode());
2034 
2035  //Photon kinematics
2036  primStack->Momentum(fPi0Mom);
2037 
2038  mesonY = 0.5*TMath::Log((primStack->Energy()+primStack->Pz())/(primStack->Energy()-primStack->Pz())) ;
2039  }
2040  else
2041  {
2042  primAOD = (AliAODMCParticle *) mcparticles->At(i);
2043  if(!primAOD)
2044  {
2045  AliWarning("AOD primaries pointer not available!!");
2046  continue;
2047  }
2048 
2049  // If too small skip
2050  if( primAOD->E() < 0.4 ) continue;
2051 
2052  pdg = primAOD->GetPdgCode();
2053  nDaught = primAOD->GetNDaughters();
2054  iphot1 = primAOD->GetFirstDaughter() ;
2055  iphot2 = primAOD->GetLastDaughter() ;
2056 
2057  // Protection against floating point exception
2058  if ( primAOD->E() == TMath::Abs(primAOD->Pz()) ||
2059  (primAOD->E() - primAOD->Pz()) < 1e-3 ||
2060  (primAOD->E() + primAOD->Pz()) < 0 ) continue ;
2061 
2062  // Photon kinematics
2063  fPi0Mom.SetPxPyPzE(primAOD->Px(),primAOD->Py(),primAOD->Pz(),primAOD->E());
2064 
2065  mesonY = 0.5*TMath::Log((primAOD->E()+primAOD->Pz())/(primAOD->E()-primAOD->Pz())) ;
2066  }
2067 
2068  // Select only pi0 or eta
2069  if( pdg != 111 && pdg != 221) continue ;
2070 
2071  mesonPt = fPi0Mom.Pt () ;
2072  mesonE = fPi0Mom.E () ;
2073  mesonYeta= fPi0Mom.Eta() ;
2074  mesonPhi = fPi0Mom.Phi() ;
2075  if( mesonPhi < 0 ) mesonPhi+=TMath::TwoPi();
2076  mesonPhi *= TMath::RadToDeg();
2077 
2078  if(pdg == 111)
2079  {
2080  if(TMath::Abs(mesonY) < 1.0)
2081  {
2082  fhPrimPi0E ->Fill(mesonE , GetEventWeight()) ;
2083  fhPrimPi0Pt ->Fill(mesonPt, GetEventWeight()) ;
2084  fhPrimPi0Phi->Fill(mesonPt, mesonPhi, GetEventWeight()) ;
2085 
2086  fhPrimPi0YetaYcut->Fill(mesonPt, mesonYeta, GetEventWeight()) ;
2087 
2089  {
2090  fhPrimPi0PtCentrality->Fill(mesonPt, cen, GetEventWeight()) ;
2091  fhPrimPi0PtEventPlane->Fill(mesonPt, ep , GetEventWeight()) ;
2092  }
2093  }
2094 
2095  fhPrimPi0Y ->Fill(mesonPt, mesonY , GetEventWeight()) ;
2096  fhPrimPi0Yeta->Fill(mesonPt, mesonYeta, GetEventWeight()) ;
2097  }
2098  else if(pdg == 221)
2099  {
2100  if(TMath::Abs(mesonY) < 1.0)
2101  {
2102  fhPrimEtaE ->Fill(mesonE , GetEventWeight()) ;
2103  fhPrimEtaPt ->Fill(mesonPt, GetEventWeight()) ;
2104  fhPrimEtaPhi->Fill(mesonPt, mesonPhi, GetEventWeight()) ;
2105 
2106  fhPrimEtaYetaYcut->Fill(mesonPt, mesonYeta, GetEventWeight()) ;
2107 
2109  {
2110  fhPrimEtaPtCentrality->Fill(mesonPt, cen, GetEventWeight()) ;
2111  fhPrimEtaPtEventPlane->Fill(mesonPt, ep , GetEventWeight()) ;
2112  }
2113  }
2114 
2115  fhPrimEtaY ->Fill(mesonPt, mesonY , GetEventWeight()) ;
2116  fhPrimEtaYeta->Fill(mesonPt, mesonYeta, GetEventWeight()) ;
2117  }
2118 
2119  // Origin of meson
2120  if(fFillOriginHisto && TMath::Abs(mesonY) < 0.7)
2121  {
2122  Int_t momindex = -1;
2123  Int_t mompdg = -1;
2124  Int_t momstatus = -1;
2125  Int_t status = -1;
2126  Float_t momR = 0;
2127  if(GetReader()->ReadStack()) momindex = primStack->GetFirstMother();
2128  if(GetReader()->ReadAODMCParticles()) momindex = primAOD ->GetMother();
2129 
2130  if(momindex >= 0 && momindex < nprim)
2131  {
2132  if(GetReader()->ReadStack())
2133  {
2134  status = primStack->GetStatusCode();
2135  TParticle* mother = stack->Particle(momindex);
2136  mompdg = TMath::Abs(mother->GetPdgCode());
2137  momstatus = mother->GetStatusCode();
2138  momR = mother->R();
2139  }
2140 
2141  if(GetReader()->ReadAODMCParticles())
2142  {
2143  status = primAOD->GetStatus();
2144  AliAODMCParticle* mother = (AliAODMCParticle*) mcparticles->At(momindex);
2145  mompdg = TMath::Abs(mother->GetPdgCode());
2146  momstatus = mother->GetStatus();
2147  momR = TMath::Sqrt(mother->Xv()*mother->Xv()+mother->Yv()*mother->Yv());
2148  }
2149 
2150  if(pdg == 111)
2151  {
2152  fhPrimPi0ProdVertex->Fill(mesonPt, momR , GetEventWeight());
2153  fhPrimPi0PtStatus ->Fill(mesonPt, status, GetEventWeight());
2154 
2155 
2156  if (momstatus == 21) fhPrimPi0PtOrigin->Fill(mesonPt, 0.5, GetEventWeight());//parton
2157  else if(mompdg < 22 ) fhPrimPi0PtOrigin->Fill(mesonPt, 1.5, GetEventWeight());//quark
2158  else if(mompdg > 2100 && mompdg < 2210)
2159  fhPrimPi0PtOrigin->Fill(mesonPt, 2.5, GetEventWeight());// resonances
2160  else if(mompdg == 221) fhPrimPi0PtOrigin->Fill(mesonPt, 8.5, GetEventWeight());//eta
2161  else if(mompdg == 331) fhPrimPi0PtOrigin->Fill(mesonPt, 9.5, GetEventWeight());//eta prime
2162  else if(mompdg == 213) fhPrimPi0PtOrigin->Fill(mesonPt, 4.5, GetEventWeight());//rho
2163  else if(mompdg == 223) fhPrimPi0PtOrigin->Fill(mesonPt, 5.5, GetEventWeight());//omega
2164  else if(mompdg >= 310 && mompdg <= 323)
2165  fhPrimPi0PtOrigin->Fill(mesonPt, 6.5, GetEventWeight());//k0S, k+-,k*
2166  else if(mompdg == 130) fhPrimPi0PtOrigin->Fill(mesonPt, 6.5, GetEventWeight());//k0L
2167  else if(momstatus == 11 || momstatus == 12 )
2168  fhPrimPi0PtOrigin->Fill(mesonPt, 3.5, GetEventWeight());//resonances
2169  else fhPrimPi0PtOrigin->Fill(mesonPt, 7.5, GetEventWeight());//other?
2170 
2171  //printf("Prim Pi0: index %d, pt %2.2f Prod vertex %3.3f, origin pdg %d, origin status %d, origin UI %d\n",
2172  // momindex, mesonPt,mother->R(),mompdg,momstatus,mother->GetUniqueID());
2173 
2174  if(status!=11)
2175  {
2176  if (momstatus == 21) fhPrimNotResonancePi0PtOrigin->Fill(mesonPt, 0.5, GetEventWeight());//parton
2177  else if(mompdg < 22 ) fhPrimNotResonancePi0PtOrigin->Fill(mesonPt, 1.5, GetEventWeight());//quark
2178  else if(mompdg > 2100 && mompdg < 2210)
2179  fhPrimNotResonancePi0PtOrigin->Fill(mesonPt, 2.5, GetEventWeight());// resonances
2180  else if(mompdg == 221) fhPrimNotResonancePi0PtOrigin->Fill(mesonPt, 8.5, GetEventWeight());//eta
2181  else if(mompdg == 331) fhPrimNotResonancePi0PtOrigin->Fill(mesonPt, 9.5, GetEventWeight());//eta prime
2182  else if(mompdg == 213) fhPrimNotResonancePi0PtOrigin->Fill(mesonPt, 4.5, GetEventWeight());//rho
2183  else if(mompdg == 223) fhPrimNotResonancePi0PtOrigin->Fill(mesonPt, 5.5, GetEventWeight());//omega
2184  else if(mompdg >= 310 && mompdg <= 323)
2185  fhPrimNotResonancePi0PtOrigin->Fill(mesonPt, 6.5, GetEventWeight());//k0S, k+-,k*
2186  else if(mompdg == 130) fhPrimNotResonancePi0PtOrigin->Fill(mesonPt, 6.5, GetEventWeight());//k0L
2187  else if(momstatus == 11 || momstatus == 12 )
2188  fhPrimNotResonancePi0PtOrigin->Fill(mesonPt, 3.5, GetEventWeight());//resonances
2189  else fhPrimNotResonancePi0PtOrigin->Fill(mesonPt, 7.5, GetEventWeight());//other?
2190  }
2191 
2192  }//pi0
2193  else
2194  {
2195  if (momstatus == 21 ) fhPrimEtaPtOrigin->Fill(mesonPt, 0.5, GetEventWeight());//parton
2196  else if(mompdg < 22 ) fhPrimEtaPtOrigin->Fill(mesonPt, 1.5, GetEventWeight());//quark
2197  else if(mompdg > 2100 && mompdg < 2210)
2198  fhPrimEtaPtOrigin->Fill(mesonPt, 2.5, GetEventWeight());//qq resonances
2199  else if(mompdg == 331) fhPrimEtaPtOrigin->Fill(mesonPt, 5.5, GetEventWeight());//eta prime
2200  else if(momstatus == 11 || momstatus == 12 )
2201  fhPrimEtaPtOrigin->Fill(mesonPt, 3.5, GetEventWeight());//resonances
2202  else fhPrimEtaPtOrigin->Fill(mesonPt, 4.5, GetEventWeight());//stable, conversions?
2203  //printf("Other Meson pdg %d, Mother %s, pdg %d, status %d\n",pdg, TDatabasePDG::Instance()->GetParticle(mompdg)->GetName(),mompdg, momstatus );
2204 
2205  fhPrimEtaProdVertex->Fill(mesonPt, momR, GetEventWeight());
2206  }
2207  } // pi0 has mother
2208  }
2209 
2210  // Check if both photons hit calorimeter or a given fiducial region
2211  // only if those settings are specified.
2212  if ( !IsFiducialCutOn() && !IsRealCaloAcceptanceOn() ) continue ;
2213 
2214  if ( nDaught != 2 ) continue; //Only interested in 2 gamma decay
2215 
2216  if ( iphot1 < 0 || iphot1 >= nprim || iphot2 < 0 || iphot2 >= nprim ) continue ;
2217 
2218  Int_t pdg1 = 0;
2219  Int_t pdg2 = 0;
2220  Bool_t inacceptance1 = kTRUE;
2221  Bool_t inacceptance2 = kTRUE;
2222 
2223  if(GetReader()->ReadStack())
2224  {
2225  TParticle * phot1 = stack->Particle(iphot1) ;
2226  TParticle * phot2 = stack->Particle(iphot2) ;
2227 
2228  if(!phot1 || !phot2) continue ;
2229 
2230  pdg1 = phot1->GetPdgCode();
2231  pdg2 = phot2->GetPdgCode();
2232 
2233  phot1->Momentum(fPhotonMom1);
2234  phot2->Momentum(fPhotonMom2);
2235 
2236  // Check if photons hit the Calorimeter acceptance
2238  {
2239  if( !GetCaloUtils()->IsMCParticleInCalorimeterAcceptance( GetCalorimeter(), phot1 )) inacceptance1 = kFALSE ;
2240  if( !GetCaloUtils()->IsMCParticleInCalorimeterAcceptance( GetCalorimeter(), phot2 )) inacceptance2 = kFALSE ;
2241  }
2242  }
2243 
2244  if(GetReader()->ReadAODMCParticles())
2245  {
2246  AliAODMCParticle * phot1 = (AliAODMCParticle *) mcparticles->At(iphot1) ;
2247  AliAODMCParticle * phot2 = (AliAODMCParticle *) mcparticles->At(iphot2) ;
2248 
2249  if(!phot1 || !phot2) continue ;
2250 
2251  pdg1 = phot1->GetPdgCode();
2252  pdg2 = phot2->GetPdgCode();
2253 
2254  fPhotonMom1.SetPxPyPzE(phot1->Px(),phot1->Py(),phot1->Pz(),phot1->E());
2255  fPhotonMom2.SetPxPyPzE(phot2->Px(),phot2->Py(),phot2->Pz(),phot2->E());
2256 
2257  // Check if photons hit the Calorimeter acceptance
2259  {
2260  if( !GetCaloUtils()->IsMCParticleInCalorimeterAcceptance( GetCalorimeter(), phot1 )) inacceptance1 = kFALSE ;
2261  if( !GetCaloUtils()->IsMCParticleInCalorimeterAcceptance( GetCalorimeter(), phot2 )) inacceptance2 = kFALSE ;
2262  }
2263  }
2264 
2265  if( pdg1 != 22 || pdg2 !=22) continue ;
2266 
2267  // Check if photons hit desired acceptance in the fidutial borders fixed in the analysis
2268  if(IsFiducialCutOn())
2269  {
2270  if( inacceptance1 && !GetFiducialCut()->IsInFiducialCut(fPhotonMom1.Eta(), fPhotonMom1.Phi(), GetCalorimeter()) ) inacceptance1 = kFALSE ;
2271  if( inacceptance2 && !GetFiducialCut()->IsInFiducialCut(fPhotonMom2.Eta(), fPhotonMom2.Phi(), GetCalorimeter()) ) inacceptance2 = kFALSE ;
2272  }
2273 
2275 
2276  if(GetCalorimeter()==kEMCAL && inacceptance1 && inacceptance2 && fCheckAccInSector)
2277  {
2278  Int_t absID1=0;
2279  Int_t absID2=0;
2280 
2281  Float_t photonPhi1 = fPhotonMom1.Phi();
2282  Float_t photonPhi2 = fPhotonMom2.Phi();
2283 
2284  if(photonPhi1 < 0) photonPhi1+=TMath::TwoPi();
2285  if(photonPhi2 < 0) photonPhi2+=TMath::TwoPi();
2286 
2287  GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(fPhotonMom1.Eta(),photonPhi1,absID1);
2288  GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(fPhotonMom2.Eta(),photonPhi2,absID2);
2289 
2290  Int_t sm1 = GetEMCALGeometry()->GetSuperModuleNumber(absID1);
2291  Int_t sm2 = GetEMCALGeometry()->GetSuperModuleNumber(absID2);
2292 
2293  Int_t j=0;
2294  Bool_t sameSector = kFALSE;
2295  for(Int_t isector = 0; isector < fNModules/2; isector++)
2296  {
2297  j=2*isector;
2298  if((sm1==j && sm2==j+1) || (sm1==j+1 && sm2==j)) sameSector = kTRUE;
2299  }
2300 
2301  if(sm1!=sm2 && !sameSector)
2302  {
2303  inacceptance1 = kFALSE;
2304  inacceptance2 = kFALSE;
2305  }
2306  //if(sm1!=sm2)printf("sm1 %d, sm2 %d, same sector %d, in acceptance %d\n",sm1,sm2,sameSector,inacceptance);
2307  // if(GetEMCALGeometry()->Impact(phot1) && GetEMCALGeometry()->Impact(phot2))
2308  // inacceptance = kTRUE;
2309  }
2310 
2311  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",
2313  mesonPt,mesonYeta,mesonPhi,
2314  fPhotonMom1.Pt(),fPhotonMom1.Eta(),fPhotonMom1.Phi()*TMath::RadToDeg(),
2315  fPhotonMom2.Pt(),fPhotonMom2.Eta(),fPhotonMom2.Phi()*TMath::RadToDeg(),
2316  inacceptance1, inacceptance2));
2317 
2318  if(inacceptance1 && inacceptance2)
2319  {
2320  Float_t asym = TMath::Abs((fPhotonMom1.E()-fPhotonMom2.E()) / (fPhotonMom1.E()+fPhotonMom2.E()));
2321  Double_t angle = fPhotonMom1.Angle(fPhotonMom2.Vect());
2322 
2323  Bool_t cutAngle = kFALSE;
2324  if(fUseAngleCut && (angle < fAngleCut || angle > fAngleMaxCut)) cutAngle = kTRUE;
2325 
2326  AliDebug(2,Form("\t ACCEPTED pdg %d: pt %2.2f, phi %2.2f, eta %2.2f",pdg,mesonPt,mesonPhi,mesonYeta));
2327 
2328  if(pdg==111)
2329  {
2330  fhPrimPi0AccE ->Fill(mesonE , GetEventWeight()) ;
2331  fhPrimPi0AccPt ->Fill(mesonPt, GetEventWeight()) ;
2332  fhPrimPi0AccPhi ->Fill(mesonPt, mesonPhi , GetEventWeight()) ;
2333  fhPrimPi0AccY ->Fill(mesonPt, mesonY , GetEventWeight()) ;
2334  fhPrimPi0AccYeta->Fill(mesonPt, mesonYeta, GetEventWeight()) ;
2335 
2337  {
2338  fhPrimPi0AccPtCentrality->Fill(mesonPt, cen, GetEventWeight()) ;
2339  fhPrimPi0AccPtEventPlane->Fill(mesonPt, ep , GetEventWeight()) ;
2340  }
2341 
2342  if(fFillAngleHisto)
2343  {
2344  fhPrimPi0OpeningAngle ->Fill(mesonPt, angle, GetEventWeight());
2345  if(mesonPt > 5)fhPrimPi0OpeningAngleAsym->Fill(asym, angle, GetEventWeight());
2346  fhPrimPi0CosOpeningAngle ->Fill(mesonPt, TMath::Cos(angle), GetEventWeight());
2347  }
2348 
2349  if(fPhotonMom1.Pt() > GetMinPt() && fPhotonMom2.Pt() > GetMinPt() && !cutAngle)
2350  {
2351  fhPrimPi0AccPtPhotonCuts->Fill(mesonPt, GetEventWeight()) ;
2352 
2353  if(fFillAngleHisto)
2354  fhPrimPi0OpeningAnglePhotonCuts->Fill(mesonPt, angle, GetEventWeight());
2355 
2356  if(fNAngleCutBins > 0)
2357  {
2358  Int_t angleBin = -1;
2359  for(Int_t ibin = 0; ibin < fNAngleCutBins; ibin++)
2360  {
2361  if(angle >= fAngleCutBinsArray[ibin] &&
2362  angle < fAngleCutBinsArray[ibin+1]) angleBin = ibin;
2363  }
2364 
2365  if( angleBin >= 0 && angleBin < fNAngleCutBins)
2366  fhPrimPi0AccPtOpAngCuts[angleBin]->Fill(mesonPt,GetEventWeight());
2367  }
2368  }
2369  }
2370  else if(pdg==221)
2371  {
2372  fhPrimEtaAccE ->Fill(mesonE , GetEventWeight()) ;
2373  fhPrimEtaAccPt ->Fill(mesonPt, GetEventWeight()) ;
2374  fhPrimEtaAccPhi ->Fill(mesonPt, mesonPhi , GetEventWeight()) ;
2375  fhPrimEtaAccY ->Fill(mesonPt, mesonY , GetEventWeight()) ;
2376  fhPrimEtaAccYeta->Fill(mesonPt, mesonYeta, GetEventWeight()) ;
2377 
2379  {
2380  fhPrimEtaAccPtCentrality->Fill(mesonPt, cen, GetEventWeight()) ;
2381  fhPrimEtaAccPtEventPlane->Fill(mesonPt, ep , GetEventWeight()) ;
2382  }
2383 
2384  if(fFillAngleHisto)
2385  {
2386  fhPrimEtaOpeningAngle ->Fill(mesonPt, angle, GetEventWeight());
2387  if(mesonPt > 5)fhPrimEtaOpeningAngleAsym->Fill(asym, angle, GetEventWeight());
2388  fhPrimEtaCosOpeningAngle ->Fill(mesonPt, TMath::Cos(angle), GetEventWeight());
2389 
2390  if(fPhotonMom1.Pt() > GetMinPt() && fPhotonMom2.Pt() > GetMinPt() && !cutAngle)
2391  {
2392  if(fUseAngleCut && angle > fAngleCut && angle < fAngleMaxCut)
2393  fhPrimEtaAccPtPhotonCuts->Fill(mesonPt, GetEventWeight()) ;
2394 
2395  if(fFillAngleHisto)
2396  fhPrimEtaOpeningAnglePhotonCuts->Fill(mesonPt, angle, GetEventWeight());
2397 
2398  if(fNAngleCutBins > 0)
2399  {
2400  Int_t angleBin = -1;
2401  for(Int_t ibin = 0; ibin < fNAngleCutBins; ibin++)
2402  {
2403  if(angle >= fAngleCutBinsArray[ibin] &&
2404  angle < fAngleCutBinsArray[ibin+1]) angleBin = ibin;
2405  }
2406 
2407  if( angleBin >= 0 && angleBin < fNAngleCutBins)
2408  fhPrimEtaAccPtOpAngCuts[angleBin]->Fill(mesonPt,GetEventWeight());
2409  }
2410  }
2411  }
2412  }
2413  } // Accepted
2414  } // loop on primaries
2415 }
2416 
2417 //________________________________________________
2419 //________________________________________________
2421 {
2422  // Get pTArm and AlphaArm
2423  Float_t momentumSquaredMother = fPi0Mom.P()*fPi0Mom.P();
2424  Float_t momentumDaughter1AlongMother = 0.;
2425  Float_t momentumDaughter2AlongMother = 0.;
2426 
2427  if (momentumSquaredMother > 0.)
2428  {
2429  momentumDaughter1AlongMother = (fPhotonMom1.Px()*fPi0Mom.Px() + fPhotonMom1.Py()*fPi0Mom.Py()+ fPhotonMom1.Pz()*fPi0Mom.Pz()) / sqrt(momentumSquaredMother);
2430  momentumDaughter2AlongMother = (fPhotonMom2.Px()*fPi0Mom.Px() + fPhotonMom2.Py()*fPi0Mom.Py()+ fPhotonMom2.Pz()*fPi0Mom.Pz()) / sqrt(momentumSquaredMother);
2431  }
2432 
2433  Float_t momentumSquaredDaughter1 = fPhotonMom1.P()*fPhotonMom1.P();
2434  Float_t ptArmSquared = momentumSquaredDaughter1 - momentumDaughter1AlongMother*momentumDaughter1AlongMother;
2435 
2436  Float_t pTArm = 0.;
2437  if (ptArmSquared > 0.)
2438  pTArm = sqrt(ptArmSquared);
2439 
2440  Float_t alphaArm = 0.;
2441  if(momentumDaughter1AlongMother + momentumDaughter2AlongMother > 0)
2442  alphaArm = (momentumDaughter1AlongMother -momentumDaughter2AlongMother) / (momentumDaughter1AlongMother + momentumDaughter2AlongMother);
2443 
2445  fPhotonMom1Boost.Boost(-fPi0Mom.BoostVector());
2446  Float_t cosThStar=TMath::Cos(fPhotonMom1Boost.Vect().Angle(fPi0Mom.Vect()));
2447 
2448  Float_t en = fPi0Mom.Energy();
2449  Int_t ebin = -1;
2450  if(en > 8 && en <= 12) ebin = 0;
2451  if(en > 12 && en <= 16) ebin = 1;
2452  if(en > 16 && en <= 20) ebin = 2;
2453  if(en > 20) ebin = 3;
2454  if(ebin < 0 || ebin > 3) return ;
2455 
2456  if(pdg==111)
2457  {
2458  fhCosThStarPrimPi0->Fill(en, cosThStar, GetEventWeight());
2459  fhArmPrimPi0[ebin]->Fill(alphaArm, pTArm , GetEventWeight());
2460  }
2461  else
2462  {
2463  fhCosThStarPrimEta->Fill(en, cosThStar, GetEventWeight());
2464  fhArmPrimEta[ebin]->Fill(alphaArm,pTArm , GetEventWeight());
2465  }
2466 
2467  // if(GetDebug() > 2 )
2468  // {
2469  // Float_t asym = 0;
2470  // if(fPhotonMom1.E()+fPhotonMom2.E() > 0) asym = TMath::Abs(fPhotonMom1.E()-fPhotonMom2.E())/(fPhotonMom1.E()+fPhotonMom2.E());
2471  //
2472  // printf("AliAnaPi0::FillArmenterosThetaStar() - E %f, alphaArm %f, pTArm %f, cos(theta*) %f, asymmetry %1.3f\n",
2473  // en,alphaArm,pTArm,cosThStar,asym);
2474  // }
2475 }
2476 
2477 //_______________________________________________________________________________________
2482 //_______________________________________________________________________________________
2484  Float_t pt1, Float_t pt2,
2485  Int_t ncell1, Int_t ncell2,
2486  Double_t mass, Double_t pt, Double_t asym,
2487  Double_t deta, Double_t dphi, Double_t angle)
2488 {
2489  Int_t ancPDG = 0;
2490  Int_t ancStatus = 0;
2491  Int_t ancLabel = GetMCAnalysisUtils()->CheckCommonAncestor(index1, index2,
2492  GetReader(), ancPDG, ancStatus,fPi0Mom, fProdVertex);
2493 
2494  Int_t momindex = -1;
2495  Int_t mompdg = -1;
2496  Int_t momstatus = -1;
2497  Int_t status = -1;
2498  Float_t prodR = -1;
2499  Int_t mcIndex = -1;
2500 
2501  if(ancLabel > -1)
2502  {
2503  AliDebug(1,Form("Common ancestor label %d, pdg %d, name %s, status %d",
2504  ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus));
2505 
2506  if(ancPDG==22)
2507  {//gamma
2508  mcIndex = 0;
2509  }
2510  else if(TMath::Abs(ancPDG)==11)
2511  {//e
2512  mcIndex = 1;
2513  }
2514  else if(ancPDG==111)
2515  {//Pi0
2516  mcIndex = 2;
2517  if(fMultiCutAnaSim)
2518  {
2519  for(Int_t ipt=0; ipt<fNPtCuts; ipt++)
2520  {
2521  for(Int_t icell=0; icell<fNCellNCuts; icell++)
2522  {
2523  for(Int_t iasym=0; iasym<fNAsymCuts; iasym++)
2524  {
2525  Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
2526  if(pt1 > fPtCuts[ipt] && pt2 > fPtCuts[ipt] &&
2527  pt1 < fPtCutsMax[ipt] && pt2 < fPtCutsMax[ipt] &&
2528  asym < fAsymCuts[iasym] &&
2529  ncell1 >= fCellNCuts[icell] && ncell2 >= fCellNCuts[icell])
2530  {
2531  fhMCPi0MassPtRec [index]->Fill(pt, mass, GetEventWeight());
2532  fhMCPi0MassPtTrue[index]->Fill(fPi0Mom.Pt(),mass, GetEventWeight());
2533  if(mass < 0.17 && mass > 0.1)
2534  fhMCPi0PtTruePtRec[index]->Fill(fPi0Mom.Pt(), pt, GetEventWeight());
2535  }//pass the different cuts
2536  }// pid bit cut loop
2537  }// icell loop
2538  }// pt cut loop
2539  }// Multi cut ana sim
2540  else
2541  {
2542  fhMCPi0MassPtTrue[0]->Fill(fPi0Mom.Pt(), mass, GetEventWeight());
2543 
2544  if(mass < 0.17 && mass > 0.1)
2545  {
2546  fhMCPi0PtTruePtRec[0]->Fill(fPi0Mom.Pt(), pt, GetEventWeight());
2547 
2548  //Int_t uniqueId = -1;
2549  if(GetReader()->ReadStack())
2550  {
2551  TParticle* ancestor = GetMCStack()->Particle(ancLabel);
2552  status = ancestor->GetStatusCode();
2553  momindex = ancestor->GetFirstMother();
2554  if(momindex < 0) return;
2555  TParticle* mother = GetMCStack()->Particle(momindex);
2556  mompdg = TMath::Abs(mother->GetPdgCode());
2557  momstatus = mother->GetStatusCode();
2558  prodR = mother->R();
2559  //uniqueId = mother->GetUniqueID();
2560  }
2561  else
2562  {
2563  TClonesArray * mcparticles = GetReader()->GetAODMCParticles();
2564  AliAODMCParticle* ancestor = (AliAODMCParticle *) mcparticles->At(ancLabel);
2565  status = ancestor->GetStatus();
2566  momindex = ancestor->GetMother();
2567  if(momindex < 0) return;
2568  AliAODMCParticle* mother = (AliAODMCParticle *) mcparticles->At(momindex);
2569  mompdg = TMath::Abs(mother->GetPdgCode());
2570  momstatus = mother->GetStatus();
2571  prodR = TMath::Sqrt(mother->Xv()*mother->Xv()+mother->Yv()*mother->Yv());
2572  //uniqueId = mother->GetUniqueID();
2573  }
2574 
2575  // printf("Reco Pi0: pt %2.2f Prod vertex %3.3f, origin pdg %d, origin status %d, origin UI %d\n",
2576  // pt,prodR,mompdg,momstatus,uniqueId);
2577 
2578  fhMCPi0ProdVertex->Fill(pt, prodR , GetEventWeight());
2579  fhMCPi0PtStatus ->Fill(pt, status, GetEventWeight());
2580 
2581  if (momstatus == 21) fhMCPi0PtOrigin->Fill(pt, 0.5, GetEventWeight());//parton
2582  else if(mompdg < 22 ) fhMCPi0PtOrigin->Fill(pt, 1.5, GetEventWeight());//quark
2583  else if(mompdg > 2100 && mompdg < 2210)
2584  fhMCPi0PtOrigin->Fill(pt, 2.5, GetEventWeight());// resonances
2585  else if(mompdg == 221) fhMCPi0PtOrigin->Fill(pt, 8.5, GetEventWeight());//eta
2586  else if(mompdg == 331) fhMCPi0PtOrigin->Fill(pt, 9.5, GetEventWeight());//eta prime
2587  else if(mompdg == 213) fhMCPi0PtOrigin->Fill(pt, 4.5, GetEventWeight());//rho
2588  else if(mompdg == 223) fhMCPi0PtOrigin->Fill(pt, 5.5, GetEventWeight());//omega
2589  else if(mompdg >= 310 && mompdg <= 323)
2590  fhMCPi0PtOrigin->Fill(pt, 6.5, GetEventWeight());//k0S, k+-,k*
2591  else if(mompdg == 130) fhMCPi0PtOrigin->Fill(pt, 6.5, GetEventWeight());//k0L
2592  else if(momstatus == 11 || momstatus == 12 )
2593  fhMCPi0PtOrigin->Fill(pt, 3.5, GetEventWeight());//resonances
2594  else fhMCPi0PtOrigin->Fill(pt, 7.5, GetEventWeight());//other?
2595 
2596  if(status!=11)
2597  {
2598  if (momstatus == 21) fhMCNotResonancePi0PtOrigin->Fill(pt, 0.5, GetEventWeight());//parton
2599  else if(mompdg < 22 ) fhMCNotResonancePi0PtOrigin->Fill(pt, 1.5, GetEventWeight());//quark
2600  else if(mompdg > 2100 && mompdg < 2210)
2601  fhMCNotResonancePi0PtOrigin->Fill(pt, 2.5, GetEventWeight());// resonances
2602  else if(mompdg == 221) fhMCNotResonancePi0PtOrigin->Fill(pt, 8.5, GetEventWeight());//eta
2603  else if(mompdg == 331) fhMCNotResonancePi0PtOrigin->Fill(pt, 9.5, GetEventWeight());//eta prime
2604  else if(mompdg == 213) fhMCNotResonancePi0PtOrigin->Fill(pt, 4.5, GetEventWeight());//rho
2605  else if(mompdg == 223) fhMCNotResonancePi0PtOrigin->Fill(pt, 5.5, GetEventWeight());//omega
2606  else if(mompdg >= 310 && mompdg <= 323)
2607  fhMCNotResonancePi0PtOrigin->Fill(pt, 6.5, GetEventWeight());//k0S, k+-,k*
2608  else if(mompdg == 130) fhMCNotResonancePi0PtOrigin->Fill(pt, 6.5, GetEventWeight());//k0L
2609  else if(momstatus == 11 || momstatus == 12 )
2610  fhMCNotResonancePi0PtOrigin->Fill(pt, 3.5, GetEventWeight());//resonances
2611  else fhMCNotResonancePi0PtOrigin->Fill(pt, 7.5, GetEventWeight());//other?
2612  }
2613  }//pi0 mass region
2614  }
2615 
2616  if(fNAngleCutBins > 0)
2617  {
2618  Int_t angleBin = -1;
2619  for(Int_t ibin = 0; ibin < fNAngleCutBins; ibin++)
2620  {
2621  if(angle >= fAngleCutBinsArray[ibin] &&
2622  angle < fAngleCutBinsArray[ibin+1]) angleBin = ibin;
2623  }
2624 
2625  if( angleBin >= 0 && angleBin < fNAngleCutBins)
2627  }
2628 
2629  }
2630  else if(ancPDG==221)
2631  {//Eta
2632  mcIndex = 3;
2633  if(fMultiCutAnaSim)
2634  {
2635  for(Int_t ipt=0; ipt<fNPtCuts; ipt++)
2636  {
2637  for(Int_t icell=0; icell<fNCellNCuts; icell++)
2638  {
2639  for(Int_t iasym=0; iasym<fNAsymCuts; iasym++)
2640  {
2641  Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
2642  if(pt1 > fPtCuts[ipt] && pt2 > fPtCuts[ipt] &&
2643  pt1 < fPtCutsMax[ipt] && pt2 < fPtCutsMax[ipt] &&
2644  asym < fAsymCuts[iasym] &&
2645  ncell1 >= fCellNCuts[icell] && ncell2 >= fCellNCuts[icell])
2646  {
2647  fhMCEtaMassPtRec [index]->Fill(pt , mass, GetEventWeight());
2648  fhMCEtaMassPtTrue[index]->Fill(fPi0Mom.Pt(), mass, GetEventWeight());
2649  if(mass < 0.65 && mass > 0.45)
2650  fhMCEtaPtTruePtRec[index]->Fill(fPi0Mom.Pt(), pt, GetEventWeight());
2651  }//pass the different cuts
2652  }// pid bit cut loop
2653  }// icell loop
2654  }// pt cut loop
2655  } //Multi cut ana sim
2656  else
2657  {
2658  fhMCEtaMassPtTrue[0]->Fill(fPi0Mom.Pt(), mass, GetEventWeight());
2659  if(mass < 0.65 && mass > 0.45) fhMCEtaPtTruePtRec[0]->Fill(fPi0Mom.Pt(), pt, GetEventWeight());
2660 
2661  if(GetReader()->ReadStack())
2662  {
2663  TParticle* ancestor = GetMCStack()->Particle(ancLabel);
2664  momindex = ancestor->GetFirstMother();
2665  if(momindex < 0) return;
2666  TParticle* mother = GetMCStack()->Particle(momindex);
2667  mompdg = TMath::Abs(mother->GetPdgCode());
2668  momstatus = mother->GetStatusCode();
2669  prodR = mother->R();
2670  }
2671  else
2672  {
2673  TClonesArray * mcparticles = GetReader()->GetAODMCParticles();
2674  AliAODMCParticle* ancestor = (AliAODMCParticle *) mcparticles->At(ancLabel);
2675  momindex = ancestor->GetMother();
2676  if(momindex < 0) return;
2677  AliAODMCParticle* mother = (AliAODMCParticle *) mcparticles->At(momindex);
2678  mompdg = TMath::Abs(mother->GetPdgCode());
2679  momstatus = mother->GetStatus();
2680  prodR = TMath::Sqrt(mother->Xv()*mother->Xv()+mother->Yv()*mother->Yv());
2681  }
2682 
2683  fhMCEtaProdVertex->Fill(pt, prodR, GetEventWeight());
2684 
2685  if (momstatus == 21 ) fhMCEtaPtOrigin->Fill(pt, 0.5, GetEventWeight());//parton
2686  else if(mompdg < 22 ) fhMCEtaPtOrigin->Fill(pt, 1.5, GetEventWeight());//quark
2687  else if(mompdg > 2100 && mompdg < 2210)
2688  fhMCEtaPtOrigin->Fill(pt, 2.5, GetEventWeight());//qq resonances
2689  else if(mompdg == 331) fhMCEtaPtOrigin->Fill(pt, 5.5, GetEventWeight());//eta prime
2690  else if(momstatus == 11 || momstatus == 12 )
2691  fhMCEtaPtOrigin->Fill(pt, 3.5, GetEventWeight());//resonances
2692  else fhMCEtaPtOrigin->Fill(pt, 4.5, GetEventWeight());//stable, conversions?
2693  //printf("Other Meson pdg %d, Mother %s, pdg %d, status %d\n",pdg, TDatabasePDG::Instance()->GetParticle(mompdg)->GetName(),mompdg, momstatus );
2694 
2695  }// eta mass region
2696 
2697  if(fNAngleCutBins > 0)
2698  {
2699  Int_t angleBin = -1;
2700  for(Int_t ibin = 0; ibin < fNAngleCutBins; ibin++)
2701  {
2702  if(angle >= fAngleCutBinsArray[ibin] &&
2703  angle < fAngleCutBinsArray[ibin+1]) angleBin = ibin;
2704  }
2705 
2706  if( angleBin >= 0 && angleBin < fNAngleCutBins)
2708  }
2709  }
2710  else if(ancPDG==-2212)
2711  {//AProton
2712  mcIndex = 4;
2713  }
2714  else if(ancPDG==-2112)
2715  {//ANeutron
2716  mcIndex = 5;
2717  }
2718  else if(TMath::Abs(ancPDG)==13)
2719  {//muons
2720  mcIndex = 6;
2721  }
2722  else if (TMath::Abs(ancPDG) > 100 && ancLabel > 7)
2723  {
2724  if(ancStatus==1)
2725  {//Stable particles, converted? not decayed resonances
2726  mcIndex = 6;
2727  }
2728  else
2729  {//resonances and other decays, more hadron conversions?
2730  mcIndex = 7;
2731  }
2732  }
2733  else
2734  {//Partons, colliding protons, strings, intermediate corrections
2735  if(ancStatus==11 || ancStatus==12)
2736  {//String fragmentation
2737  mcIndex = 8;
2738  }
2739  else if (ancStatus==21){
2740  if(ancLabel < 2)
2741  {//Colliding protons
2742  mcIndex = 11;
2743  }//colliding protons
2744  else if(ancLabel < 6)
2745  {//partonic initial states interactions
2746  mcIndex = 9;
2747  }
2748  else if(ancLabel < 8)
2749  {//Final state partons radiations?
2750  mcIndex = 10;
2751  }
2752  // else {
2753  // printf("AliAnaPi0::FillMCVersusRecDataHistograms() - Check ** Common ancestor label %d, pdg %d, name %s, status %d; \n",
2754  // ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus);
2755  // }
2756  }//status 21
2757  //else {
2758  // printf("AliAnaPi0::FillMCVersusRecDataHistograms() - Check *** Common ancestor label %d, pdg %d, name %s, status %d; \n",
2759  // ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus);
2760  // }
2761  }
2762  }//ancLabel > -1
2763  else
2764  { //ancLabel <= -1
2765  //printf("Not related at all label = %d\n",ancLabel);
2766  AliDebug(1,"Common ancestor not found");
2767 
2768  mcIndex = 12;
2769  }
2770 
2771  if(mcIndex >= 0 && mcIndex < 13)
2772  {
2773  fhMCOrgMass [mcIndex]->Fill(pt, mass, GetEventWeight());
2774  fhMCOrgAsym [mcIndex]->Fill(pt, asym, GetEventWeight());
2775  fhMCOrgDeltaEta[mcIndex]->Fill(pt, deta, GetEventWeight());
2776  fhMCOrgDeltaPhi[mcIndex]->Fill(pt, dphi, GetEventWeight());
2777  }
2778 }
2779 
2780 //__________________________________________
2784 //__________________________________________
2786 {
2787  // In case of simulated data, fill acceptance histograms
2788  if(IsDataMC())
2789  {
2791 
2792  if(fFillOnlyMCAcceptanceHisto) return;
2793  }
2794 
2795 //if (GetReader()->GetEventNumber()%10000 == 0)
2796 // printf("--- Event %d ---\n",GetReader()->GetEventNumber());
2797 
2798  if(!GetInputAODBranch())
2799  {
2800  AliFatal(Form("No input aod photons in AOD with name branch < %s >, STOP",GetInputAODName().Data()));
2801  return;
2802  }
2803 
2804  //
2805  // Init some variables
2806  //
2807 
2808  // Analysis within the same detector:
2809  TClonesArray* secondLoopInputData = GetInputAODBranch();
2810 
2811  Int_t nPhot = GetInputAODBranch()->GetEntriesFast() ;
2812  Int_t nPhot2 = nPhot; // second loop
2813  Int_t minEntries = 2;
2814  Int_t last = 1; // last entry in first loop to be removed
2815 
2816  // Combine 2 detectors:
2818  {
2819  // Input from CaloTrackCorr tasks
2820  secondLoopInputData = (TClonesArray *) GetReader()->GetAODBranchList()->FindObject(fOtherDetectorInputName);
2821 
2822  // In case input from PCM or other external source
2823  if(!secondLoopInputData) secondLoopInputData = (TClonesArray *) GetReader()->GetOutputEvent()->FindObject(fOtherDetectorInputName);
2824  if(!secondLoopInputData) secondLoopInputData = (TClonesArray *) GetReader()->GetInputEvent ()->FindObject(fOtherDetectorInputName);
2825 
2826  if(!secondLoopInputData)
2827  {
2828  AliFatal(Form("Input for other detector not found in branch %s!",fOtherDetectorInputName.Data()));
2829  return ; // coverity
2830  }
2831 
2832  nPhot2 = secondLoopInputData->GetEntriesFast() ; // add one since we remove one for the normal case
2833  minEntries = 1;
2834  last = 0;
2835  }
2836 
2837  AliDebug(1,Form("Photon entries %d", nPhot));
2838 
2839  // If less than photon 2 (1) entries in the list, skip this event
2840  if ( nPhot < minEntries )
2841  {
2842  AliDebug(1,Form("nPhotons %d, cent bin %d continue to next event",nPhot, GetEventCentralityBin()));
2843 
2844  if ( GetNCentrBin() > 1 && IsHighMultiplicityAnalysisOn() )
2846 
2847  return ;
2848  }
2849  else if(fPairWithOtherDetector && nPhot2 < minEntries)
2850  {
2851  AliDebug(1,Form("nPhotons %d, cent bin %d continue to next event",nPhot, GetEventCentralityBin()));
2852 
2853  if ( GetNCentrBin() > 1 && IsHighMultiplicityAnalysisOn() )
2855 
2856  return ;
2857  }
2858 
2859  Int_t ncentr = GetNCentrBin();
2860  if(GetNCentrBin()==0) ncentr = 1;
2861 
2862  //Init variables
2863  Int_t module1 = -1;
2864  Int_t module2 = -1;
2865  Double_t vert[] = {0.0, 0.0, 0.0} ; //vertex
2866  Int_t evtIndex1 = 0 ;
2867  Int_t currentEvtIndex = -1;
2868  Int_t curCentrBin = GetEventCentralityBin();
2869  //Int_t curVzBin = GetEventVzBin();
2870  //Int_t curRPBin = GetEventRPBin();
2871  Int_t eventbin = GetEventMixBin();
2872 
2873  if(eventbin > GetNCentrBin()*GetNZvertBin()*GetNRPBin())
2874  {
2875  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",
2877  return;
2878  }
2879 
2880  // Get shower shape information of clusters
2881 // TObjArray *clusters = 0;
2882 // if (GetCalorimeter()==kEMCAL) clusters = GetEMCALClusters();
2883 // else if(GetCalorimeter()==kPHOS ) clusters = GetPHOSClusters() ;
2884 
2885  //---------------------------------
2886  // First loop on photons/clusters
2887  //---------------------------------
2888  for(Int_t i1 = 0; i1 < nPhot-last; i1++)
2889  {
2890  AliAODPWG4Particle * p1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i1)) ;
2891 
2892  // Select photons within a pT range
2893  if ( p1->Pt() < GetMinPt() || p1->Pt() > GetMaxPt() ) continue ;
2894 
2895  //printf("AliAnaPi0::MakeAnalysisFillHistograms() : cluster1 id %d/%d\n",i1,nPhot-1);
2896 
2897  // get the event index in the mixed buffer where the photon comes from
2898  // in case of mixing with analysis frame, not own mixing
2899  evtIndex1 = GetEventIndex(p1, vert) ;
2900  if ( evtIndex1 == -1 )
2901  return ;
2902  if ( evtIndex1 == -2 )
2903  continue ;
2904 
2905  // Only effective in case of mixed event frame
2906  if(TMath::Abs(vert[2]) > GetZvertexCut()) continue ; //vertex cut
2907 
2908  if (evtIndex1 != currentEvtIndex)
2909  {
2910  // Fill event bin info
2911  if( DoOwnMix() ) fhEventBin->Fill(eventbin, GetEventWeight()) ;
2912 
2914  {
2915  fhCentrality->Fill(curCentrBin, GetEventWeight());
2916 
2917  if( GetEventPlane() )
2918  fhEventPlaneResolution->Fill(curCentrBin, TMath::Cos(2.*GetEventPlane()->GetQsubRes()), GetEventWeight());
2919  }
2920 
2921  currentEvtIndex = evtIndex1 ;
2922  }
2923 
2924  //printf("AliAnaPi0::MakeAnalysisFillHistograms(): Photon 1 Evt %d Vertex : %f,%f,%f\n",evtIndex1, GetVertex(evtIndex1)[0] ,GetVertex(evtIndex1)[1],GetVertex(evtIndex1)[2]);
2925 
2926  //Get the momentum of this cluster
2927  fPhotonMom1.SetPxPyPzE(p1->Px(),p1->Py(),p1->Pz(),p1->E());
2928 
2929  //Get (Super)Module number of this cluster
2930  module1 = p1->GetSModNumber();// GetModuleNumber(p1);
2931 
2932  //------------------------------------------
2933  // Recover original cluster
2934  // Declare variables for absid col-row identification of pair
2935  // Also fill col-row, eta-phi histograms depending on pt bin
2936 
2937  Int_t iclus1 = -1, iclus2 = -1 ;
2938  Float_t maxCellFraction1 = 0, maxCellFraction2 = 0;
2939  Int_t absIdMax1 = -1, absIdMax2 = -1;
2940  Int_t icol1 = -1, icol2 = -1, icolAbs1 = -1, icolAbs2 = -1;
2941  Int_t irow1 = -1, irow2 = -1, irowAbs1 = -1, irowAbs2 = -1;
2942  Int_t iRCU1 = -1, iRCU2 = -1;
2943 
2945  {
2946  AliVCluster * cluster1 = FindCluster(GetEMCALClusters(),p1->GetCaloLabel(0),iclus1);
2947  if(!cluster1) AliWarning("Cluster1 not found!");
2948 
2949  absIdMax1 = GetCaloUtils()->GetMaxEnergyCell(GetEMCALCells(),cluster1,maxCellFraction1);
2950 
2951  GetModuleNumberCellIndexesAbsCaloMap(absIdMax1,GetCalorimeter(), icol1, irow1, iRCU1, icolAbs1, irowAbs1);
2952 
2953  if(fMultiCutAnaAcc)
2954  {
2955  for(Int_t ipt = 0; ipt < fNPtCuts; ipt++)
2956  {
2957  if( p1->Pt() > fPtCuts[ipt] && p1->Pt() < fPtCuts[ipt+1] )
2958  {
2959  fhPtBinClusterEtaPhi[ipt]->Fill(p1->Eta(),GetPhi(p1->Phi()),GetEventWeight()) ;
2960 
2961  fhPtBinClusterColRow[ipt]->Fill(icolAbs1,irowAbs1,GetEventWeight()) ;
2962  }
2963  }
2964  }
2965  }
2966 
2967  //---------------------------------
2968  // Second loop on photons/clusters
2969  //---------------------------------
2970  Int_t first = i1+1;
2971  if(fPairWithOtherDetector) first = 0;
2972 
2973  for(Int_t i2 = first; i2 < nPhot2; i2++)
2974  {
2975  //AliAODPWG4Particle * p2 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i2)) ;
2976  AliAODPWG4Particle * p2 = (AliAODPWG4Particle*) (secondLoopInputData->At(i2)) ;
2977 
2978  // Select photons within a pT range
2979  if ( p2->Pt() < GetMinPt() || p2->Pt() > GetMaxPt() ) continue ;
2980 
2981  //printf("AliAnaPi0::MakeAnalysisFillHistograms() : cluster2 i %d/%d\n",i2,nPhot);
2982 
2983  //In case of mixing frame, check we are not in the same event as the first cluster
2984  Int_t evtIndex2 = GetEventIndex(p2, vert) ;
2985  if ( evtIndex2 == -1 )
2986  return ;
2987  if ( evtIndex2 == -2 )
2988  continue ;
2989  if (GetMixedEvent() && (evtIndex1 == evtIndex2))
2990  continue ;
2991 
2992 // //------------------------------------------
2993 // // Recover original cluster
2994 // Int_t iclus2 = -1;
2995 // AliVCluster * cluster2 = FindCluster(clusters,p2->GetCaloLabel(0),iclus2,iclus1+1);
2996 // // start new loop from iclus1+1 to gain some time
2997 // if(!cluster2) AliWarning("Cluster2 not found!");
2998 //
2999 // // Get the TOF,l0 and ncells from the clusters
3000 // Float_t tof1 = -1;
3001 // Float_t l01 = -1;
3002 // Int_t ncell1 = 0;
3003 // if(cluster1)
3004 // {
3005 // tof1 = cluster1->GetTOF()*1e9;
3006 // l01 = cluster1->GetM02();
3007 // ncell1 = cluster1->GetNCells();
3008 // //printf("cluster1: E %2.2f (%2.2f), l0 %2.2f, tof %2.2f\n",cluster1->E(),p1->E(),l01,tof1);
3009 // }
3010 // //else printf("cluster1 not available: calo label %d / %d, cluster ID %d\n",
3011 // // p1->GetCaloLabel(0),(GetReader()->GetInputEvent())->GetNumberOfCaloClusters()-1,cluster1->GetID());
3012 //
3013 // Float_t tof2 = -1;
3014 // Float_t l02 = -1;
3015 // Int_t ncell2 = 0;
3016 // if(cluster2)
3017 // {
3018 // tof2 = cluster2->GetTOF()*1e9;
3019 // l02 = cluster2->GetM02();
3020 // ncell2 = cluster2->GetNCells();
3021 // //printf("cluster2: E %2.2f (%2.2f), l0 %2.2f, tof %2.2f\n",cluster2->E(),p2->E(),l02,tof2);
3022 // }
3023 // //else printf("cluster2 not available: calo label %d / %d, cluster ID %d\n",
3024 // // p2->GetCaloLabel(0),(GetReader()->GetInputEvent())->GetNumberOfCaloClusters()-1,cluster2->GetID());
3025 //
3026 // if(cluster1 && cluster2)
3027 // {
3028 // Double_t t12diff = tof1-tof2;
3029 // if(TMath::Abs(t12diff) > GetPairTimeCut()) continue;
3030 // }
3031 
3032  Float_t tof1 = p1->GetTime();
3033  Float_t l01 = p1->GetM02();
3034  Int_t ncell1 = p1->GetNCells();
3035  //printf("cluster1: E %2.2f, l0 %2.2f, tof %2.2f\n",p1->E(),l01,tof1);
3036 
3037  Float_t tof2 = p2->GetTime();
3038  Float_t l02 = p2->GetM02();
3039  Int_t ncell2 = p2->GetNCells();
3040  //printf("cluster2: E %2.2f, l0 %2.2f, tof %2.2f\n",p2->E(),l02,tof2);
3041 
3042  Double_t t12diff = tof1-tof2;
3043  fhEPairDiffTime->Fill((fPhotonMom1 + fPhotonMom2).Pt(), t12diff, GetEventWeight());
3044  if(TMath::Abs(t12diff) > GetPairTimeCut()) continue;
3045 
3046  //------------------------------------------
3047 
3048  //printf("AliAnaPi0::MakeAnalysisFillHistograms(): Photon 2 Evt %d Vertex : %f,%f,%f\n",evtIndex2, GetVertex(evtIndex2)[0] ,GetVertex(evtIndex2)[1],GetVertex(evtIndex2)[2]);
3049 
3050  // Get the momentum of this cluster
3051  fPhotonMom2.SetPxPyPzE(p2->Px(),p2->Py(),p2->Pz(),p2->E());
3052 
3053  // Get module number
3054  module2 = p2->GetSModNumber(); //GetModuleNumber(p2);
3055 
3056  //---------------------------------
3057  // Get pair kinematics
3058  //---------------------------------
3059  Double_t m = (fPhotonMom1 + fPhotonMom2).M() ;
3060  Double_t pt = (fPhotonMom1 + fPhotonMom2).Pt();
3061  Double_t deta = fPhotonMom1.Eta() - fPhotonMom2.Eta();
3062  Double_t dphi = fPhotonMom1.Phi() - fPhotonMom2.Phi();
3063  Double_t a = TMath::Abs(p1->E()-p2->E())/(p1->E()+p2->E()) ;
3064 
3065  AliDebug(2,Form("E: fPhotonMom1 %f, fPhotonMom2 %f; Pair: pT %f, mass %f, a %f", p1->E(), p2->E(), (fPhotonMom1 + fPhotonMom2).E(),m,a));
3066 
3067  //--------------------------------
3068  // Opening angle selection
3069  //--------------------------------
3070  // Check if opening angle is too large or too small compared to what is expected
3071  Double_t angle = fPhotonMom1.Angle(fPhotonMom2.Vect());
3073  {
3074  AliDebug(2,Form("Real pair angle %f (deg) not in E %f window",RadToDeg(angle), (fPhotonMom1+fPhotonMom2).E()));
3075  continue;
3076  }
3077 
3078  if(fUseAngleCut && (angle < fAngleCut || angle > fAngleMaxCut))
3079  {
3080  AliDebug(2,Form("Real pair cut %f < angle %f < cut %f (deg)",RadToDeg(fAngleCut), RadToDeg(angle), RadToDeg(fAngleMaxCut)));
3081  continue;
3082  }
3083 
3084  //-------------------------------------------------------------------------------------------------
3085  // Fill module dependent histograms, put a cut on assymmetry on the first available cut in the array
3086  //-------------------------------------------------------------------------------------------------
3087  if(a < fAsymCuts[0] && fFillSMCombinations)
3088  {
3090  {
3091  if(module1==module2 && module1 >=0 && module1<fNModules)
3092  {
3093  fhReMod[module1]->Fill(pt, m, GetEventWeight()) ;
3094  if(fFillAngleHisto) fhRealOpeningAnglePerSM[module1]->Fill(pt, angle, GetEventWeight());
3095  }
3096 
3097  if (GetCalorimeter() == kEMCAL )
3098  {
3099  // Same sector
3100  Int_t j=0;
3101  for(Int_t i = 0; i < fNModules/2; i++)
3102  {
3103  j=2*i;
3104  if((module1==j && module2==j+1) || (module1==j+1 && module2==j)) fhReSameSectorEMCALMod[i]->Fill(pt, m, GetEventWeight()) ;
3105  }
3106 
3107  // Same side
3108  for(Int_t i = 0; i < fNModules-2; i++){
3109  if((module1==i && module2==i+2) || (module1==i+2 && module2==i)) fhReSameSideEMCALMod[i]->Fill(pt, m, GetEventWeight());
3110  }
3111  } // EMCAL
3112  else
3113  { // PHOS
3114  if((module1==0 && module2==1) || (module1==1 && module2==0)) fhReDiffPHOSMod[0]->Fill(pt, m, GetEventWeight()) ;
3115  if((module1==0 && module2==2) || (module1==2 && module2==0)) fhReDiffPHOSMod[1]->Fill(pt, m, GetEventWeight()) ;
3116  if((module1==1 && module2==2) || (module1==2 && module2==1)) fhReDiffPHOSMod[2]->Fill(pt, m, GetEventWeight()) ;
3117  } // PHOS
3118  }
3119  else
3120  {
3121  Float_t phi1 = GetPhi(fPhotonMom1.Phi());
3122  Float_t phi2 = GetPhi(fPhotonMom2.Phi());
3123  Bool_t etaside = 0;
3124  if( (p1->GetDetectorTag()==kEMCAL && fPhotonMom1.Eta() < 0)
3125  || (p2->GetDetectorTag()==kEMCAL && fPhotonMom2.Eta() < 0)) etaside = 1;
3126 
3127  if ( phi1 > DegToRad(260) && phi2 > DegToRad(260) && phi1 < DegToRad(280) && phi2 < DegToRad(280)) fhReSameSectorDCALPHOSMod[0+etaside]->Fill(pt, m, GetEventWeight());
3128  else if ( phi1 > DegToRad(280) && phi2 > DegToRad(280) && phi1 < DegToRad(300) && phi2 < DegToRad(300)) fhReSameSectorDCALPHOSMod[2+etaside]->Fill(pt, m, GetEventWeight());
3129  else if ( phi1 > DegToRad(300) && phi2 > DegToRad(300) && phi1 < DegToRad(320) && phi2 < DegToRad(320)) fhReSameSectorDCALPHOSMod[4+etaside]->Fill(pt, m, GetEventWeight());
3130  else if ( (phi1 > DegToRad(260) && phi2 > DegToRad(280) && phi1 < DegToRad(280) && phi2 < DegToRad(300))
3131  || (phi1 > DegToRad(280) && phi2 > DegToRad(260) && phi1 < DegToRad(300) && phi2 < DegToRad(280))) fhReDiffSectorDCALPHOSMod[0+etaside]->Fill(pt, m, GetEventWeight());
3132  else if ( (phi1 > DegToRad(280) && phi2 > DegToRad(300) && phi1 < DegToRad(300) && phi2 < DegToRad(320))
3133  || (phi1 > DegToRad(300) && phi2 > DegToRad(280) && phi1 < DegToRad(320) && phi2 < DegToRad(300))) fhReDiffSectorDCALPHOSMod[2+etaside]->Fill(pt, m, GetEventWeight());
3134  else if ( (phi1 > DegToRad(260) && phi2 > DegToRad(300) && phi1 < DegToRad(280) && phi2 < DegToRad(320))
3135  || (phi1 > DegToRad(300) && phi2 > DegToRad(260) && phi1 < DegToRad(320) && phi2 < DegToRad(280))) fhReDiffSectorDCALPHOSMod[4+etaside]->Fill(pt, m, GetEventWeight());
3136  else fhReDiffSectorDCALPHOSMod[6+etaside]->Fill(pt, m, GetEventWeight());
3137  }
3138  } // Fill SM combinations
3139 
3140  // In case we want only pairs in same (super) module, check their origin.
3141  Bool_t ok = kTRUE;
3142 
3143  if(fSameSM)
3144  {
3146  {
3147  if(module1!=module2) ok=kFALSE;
3148  }
3149  else // PHOS and DCal in same sector
3150  {
3151  Float_t phi1 = GetPhi(fPhotonMom1.Phi());
3152  Float_t phi2 = GetPhi(fPhotonMom2.Phi());
3153  ok=kFALSE;
3154  if ( phi1 > DegToRad(260) && phi2 > DegToRad(260) && phi1 < DegToRad(280) && phi2 < DegToRad(280)) ok = kTRUE;
3155  else if ( phi1 > DegToRad(280) && phi2 > DegToRad(280) && phi1 < DegToRad(300) && phi2 < DegToRad(300)) ok = kTRUE;
3156  else if ( phi1 > DegToRad(300) && phi2 > DegToRad(300) && phi1 < DegToRad(320) && phi2 < DegToRad(320)) ok = kTRUE;
3157  }
3158  } // Pair only in same SM
3159 
3160  if(!ok) continue;
3161 
3162  //
3163  // Fill histograms with selected cluster pairs
3164  //
3165 
3166  // Check if one of the clusters comes from a conversion
3167  if(fCheckConversion)
3168  {
3169  if (p1->IsTagged() && p2->IsTagged()) fhReConv2->Fill(pt, m, GetEventWeight());
3170  else if(p1->IsTagged() || p2->IsTagged()) fhReConv ->Fill(pt, m, GetEventWeight());
3171  }
3172 
3173  // Fill shower shape cut histograms
3175  {
3176  if ( l01 > 0.01 && l01 < 0.4 &&
3177  l02 > 0.01 && l02 < 0.4 ) fhReSS[0]->Fill(pt, m, GetEventWeight()); // Tight
3178  else if( l01 > 0.4 && l02 > 0.4 ) fhReSS[1]->Fill(pt, m, GetEventWeight()); // Loose
3179  else if( l01 > 0.01 && l01 < 0.4 && l02 > 0.4 ) fhReSS[2]->Fill(pt, m, GetEventWeight()); // Both
3180  else if( l02 > 0.01 && l02 < 0.4 && l01 > 0.4 ) fhReSS[2]->Fill(pt, m, GetEventWeight()); // Both
3181  }
3182 
3183  // Fill histograms for different bad channel distance, centrality, assymmetry cut and pid bit
3184  for(Int_t ipid=0; ipid<fNPIDBits; ipid++)
3185  {
3186  if((p1->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton)) && (p2->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton)))
3187  {
3188  for(Int_t iasym=0; iasym < fNAsymCuts; iasym++)
3189  {
3190  if(a < fAsymCuts[iasym])
3191  {
3192  Int_t index = ((curCentrBin*fNPIDBits)+ipid)*fNAsymCuts + iasym;
3193  //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);
3194 
3195  if(index < 0 || index >= ncentr*fNPIDBits*fNAsymCuts) continue ;
3196 
3197  fhRe1 [index]->Fill(pt, m, GetEventWeight());
3198  if(fMakeInvPtPlots)fhReInvPt1[index]->Fill(pt, m, 1./pt * GetEventWeight()) ;
3199 
3200  if(fFillBadDistHisto)
3201  {
3202  if(p1->DistToBad()>0 && p2->DistToBad()>0)
3203  {
3204  fhRe2 [index]->Fill(pt, m, GetEventWeight()) ;
3205  if(fMakeInvPtPlots)fhReInvPt2[index]->Fill(pt, m, 1./pt * GetEventWeight()) ;
3206 
3207  if(p1->DistToBad()>1 && p2->DistToBad()>1)
3208  {
3209  fhRe3 [index]->Fill(pt, m, GetEventWeight()) ;
3210  if(fMakeInvPtPlots)fhReInvPt3[index]->Fill(pt, m, 1./pt * GetEventWeight()) ;
3211  }// bad 3
3212  }// bad2
3213  }// Fill bad dist histos
3214  }//assymetry cut
3215  }// asymmetry cut loop
3216  }// bad 1
3217  }// pid bit loop
3218 
3219  //
3220  // Fill histograms with opening angle
3221  if(fFillAngleHisto)
3222  {
3223  fhRealOpeningAngle ->Fill(pt, angle, GetEventWeight());
3224  fhRealCosOpeningAngle->Fill(pt, TMath::Cos(angle), GetEventWeight());
3225  }
3226 
3227  //
3228  // Fill histograms for different opening angle bins
3230  {
3231  Int_t angleBin = -1;
3232  for(Int_t ibin = 0; ibin < fNAngleCutBins; ibin++)
3233  {
3234  if(angle >= fAngleCutBinsArray[ibin] &&
3235  angle < fAngleCutBinsArray[ibin+1]) angleBin = ibin;
3236  }
3237 
3238  if( angleBin >= 0 && angleBin < fNAngleCutBins)
3239  {
3240  Float_t e1 = fPhotonMom1.E();
3241  Float_t e2 = fPhotonMom2.E();
3242 
3243  Float_t t1 = tof1;
3244  Float_t t2 = tof2;
3245 
3246  Int_t nc1 = ncell1;
3247  Int_t nc2 = ncell2;
3248 
3249  Float_t eta1 = fPhotonMom1.Eta();
3250  Float_t eta2 = fPhotonMom2.Eta();
3251 
3252  Float_t phi1 = GetPhi(fPhotonMom1.Phi());
3253  Float_t phi2 = GetPhi(fPhotonMom2.Phi());
3254 
3255  Int_t mod1 = module1;
3256  Int_t mod2 = module2;
3257 
3258  // Recover original cluster
3259  AliVCluster * cluster2 = FindCluster(GetEMCALClusters(),p2->GetCaloLabel(0),iclus2);
3260  if(!cluster2) AliWarning("Cluster2 not found!");
3261 
3262  absIdMax2 = GetCaloUtils()->GetMaxEnergyCell(GetEMCALCells(),cluster2,maxCellFraction2);
3263 
3264  if(e2 > e1)
3265  {
3266  e1 = fPhotonMom2.E();
3267  e2 = fPhotonMom1.E();
3268 
3269  t1 = tof2;
3270  t2 = tof1;
3271 
3272  nc1 = ncell2;
3273  nc2 = ncell1;
3274 
3275  eta1 = fPhotonMom2.Eta();
3276  eta2 = fPhotonMom1.Eta();
3277 
3278  phi1 = GetPhi(fPhotonMom2.Phi());
3279  phi2 = GetPhi(fPhotonMom1.Phi());
3280 
3281  mod1 = module2;
3282  mod2 = module1;
3283 
3284  Int_t tmp = absIdMax2;
3285  absIdMax2 = absIdMax1;
3286  absIdMax1 = tmp;
3287  }
3288 
3289  fhReOpAngleBinMinClusterEPerSM[angleBin]->Fill(e2,mod2,GetEventWeight()) ;
3290  fhReOpAngleBinMaxClusterEPerSM[angleBin]->Fill(e1,mod1,GetEventWeight()) ;
3291 
3292  fhReOpAngleBinMinClusterTimePerSM[angleBin]->Fill(t2,mod2,GetEventWeight()) ;
3293  fhReOpAngleBinMaxClusterTimePerSM[angleBin]->Fill(t1,mod1,GetEventWeight()) ;
3294 
3295  fhReOpAngleBinMinClusterNCellPerSM[angleBin]->Fill(nc2,mod2,GetEventWeight()) ;
3296  fhReOpAngleBinMaxClusterNCellPerSM[angleBin]->Fill(nc1,mod1,GetEventWeight()) ;
3297 
3298  fhReOpAngleBinPairClusterMass[angleBin]->Fill(pt,m,GetEventWeight()) ;
3299  if(mod2 == mod1) fhReOpAngleBinPairClusterMassPerSM[angleBin]->Fill(m,mod1,GetEventWeight()) ;
3300 
3301  if(e1 > 0.01) fhReOpAngleBinPairClusterRatioPerSM[angleBin]->Fill(e2/e1,mod1,GetEventWeight()) ;
3302 
3303  fhReOpAngleBinMinClusterEtaPhi[angleBin]->Fill(eta2,phi2,GetEventWeight()) ;
3304  fhReOpAngleBinMaxClusterEtaPhi[angleBin]->Fill(eta1,phi1,GetEventWeight()) ;
3305 
3306  GetModuleNumberCellIndexesAbsCaloMap(absIdMax2,GetCalorimeter(), icol2, irow2, iRCU2, icolAbs2, irowAbs2);
3307 
3308  //fhReOpAngleBinPairClusterAbsIdMaxCell[angleBin]->Fill(absIdMax1,absIdMax2,GetEventWeight());
3309 
3310  fhReOpAngleBinMinClusterColRow[angleBin]->Fill(icolAbs2,irowAbs2,GetEventWeight()) ;
3311  fhReOpAngleBinMaxClusterColRow[angleBin]->Fill(icolAbs1,irowAbs1,GetEventWeight()) ;
3312  }
3313  } // fFillOpAngleHisto
3314 
3315 
3316  // Fill histograms with pair assymmetry
3318  {
3319  fhRePtAsym->Fill(pt, a, GetEventWeight());
3320  if(m > 0.10 && m < 0.17) fhRePtAsymPi0->Fill(pt, a, GetEventWeight());
3321  if(m > 0.45 && m < 0.65) fhRePtAsymEta->Fill(pt, a, GetEventWeight());
3322  }
3323 
3324  // Check cell time content in cluster
3326  {
3327  if ( p1->GetFiducialArea() == 0 && p2->GetFiducialArea() == 0 )
3329 
3330  else if ( p1->GetFiducialArea() != 0 && p2->GetFiducialArea() != 0 )
3332  }
3333 
3334  //---------
3335  // MC data
3336  //---------
3337  // Do some MC checks on the origin of the pair, is there any common ancestor and if there is one, who?
3338  if(IsDataMC())
3339  {
3340  if(GetMCAnalysisUtils()->CheckTagBit(p1->GetTag(),AliMCAnalysisUtils::kMCConversion) &&
3342  {
3343  fhReMCFromConversion->Fill(pt, m, GetEventWeight());
3344  }
3345  else if(!GetMCAnalysisUtils()->CheckTagBit(p1->GetTag(),AliMCAnalysisUtils::kMCConversion) &&
3347  {
3348  fhReMCFromNotConversion->Fill(pt, m, GetEventWeight());
3349  }
3350  else
3351  {
3352  fhReMCFromMixConversion->Fill(pt, m, GetEventWeight());
3353  }
3354 
3355  if(fFillOriginHisto)
3356  FillMCVersusRecDataHistograms(p1->GetLabel(), p2->GetLabel(),p1->Pt(), p2->Pt(),
3357  ncell1, ncell2, m, pt, a,deta, dphi, angle);
3358  }
3359 
3360  //-----------------------
3361  // Multi cuts analysis
3362  //-----------------------
3363  if(fMultiCutAna)
3364  {
3365 // // Histograms for different PID bits selection
3366 // for(Int_t ipid=0; ipid<fNPIDBits; ipid++)
3367 // {
3368 // if(p1->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton) &&
3369 // p2->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton)) fhRePIDBits[ipid]->Fill(pt, m, GetEventWeight()) ;
3370 //
3371 // //printf("RE ipid%d, name %s\n",ipid, fhRePIDBits[ipid]->GetName());
3372 // } // pid bit cut loop
3373 
3374  // Several pt,ncell and asymmetry cuts
3375  for(Int_t ipt = 0; ipt < fNPtCuts; ipt++)
3376  {
3377  for(Int_t icell = 0; icell < fNCellNCuts; icell++)
3378  {
3379  for(Int_t iasym = 0; iasym < fNAsymCuts; iasym++)
3380  {
3381  Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
3382  if(p1->Pt() > fPtCuts[ipt] && p2->Pt() > fPtCuts[ipt] &&
3383  p1->Pt() < fPtCutsMax[ipt] && p2->Pt() < fPtCutsMax[ipt] &&
3384  a < fAsymCuts[iasym] &&
3385  ncell1 >= fCellNCuts[icell] && ncell2 >= fCellNCuts[icell])
3386  {
3387  fhRePtNCellAsymCuts[index]->Fill(pt, m, GetEventWeight()) ;
3388  if(fFillAngleHisto) fhRePtNCellAsymCutsOpAngle[index]->Fill(pt, angle, GetEventWeight()) ;
3389 
3390  if(fFillSMCombinations && module1==module2)
3391  {
3392  fhRePtNCellAsymCutsSM[module1][index]->Fill(pt, m, GetEventWeight()) ;
3393  if(fFillAngleHisto) fhRePtNCellAsymCutsSMOpAngle[module1][index]->Fill(pt, angle, GetEventWeight()) ;
3394  }
3395  }
3396  }// pid bit cut loop
3397  }// icell loop
3398  }// pt cut loop
3399 //
3400 // if(GetHistogramRanges()->GetHistoTrackMultiplicityBins())
3401 // {
3402 // for(Int_t iasym = 0; iasym < fNAsymCuts; iasym++)
3403 // {
3404 // if(a < fAsymCuts[iasym]) fhRePtMult[iasym]->Fill(pt, GetTrackMultiplicity(), m, GetEventWeight()) ;
3405 // }
3406 // }
3407  }// multiple cuts analysis
3408 
3409  }// second same event particle
3410  }// first cluster
3411 
3412  //-------------------------------------------------------------
3413  // Mixing
3414  //-------------------------------------------------------------
3415  if(DoOwnMix())
3416  {
3417  // Recover events in with same characteristics as the current event
3418 
3419  // Check that the bin exists, if not (bad determination of RP, centrality or vz bin) do nothing
3420  if(eventbin < 0) return ;
3421 
3422  TList * evMixList=fEventsList[eventbin] ;
3423 
3424  if(!evMixList)
3425  {
3426  AliWarning(Form("Mix event list not available, bin %d",eventbin));
3427  return;
3428  }
3429 
3430  Int_t nMixed = evMixList->GetSize() ;
3431  for(Int_t ii=0; ii<nMixed; ii++)
3432  {
3433  TClonesArray* ev2= (TClonesArray*) (evMixList->At(ii));
3434  Int_t nPhot2=ev2->GetEntriesFast() ;
3435  Double_t m = -999;
3436  AliDebug(1,Form("Mixed event %d photon entries %d, centrality bin %d",ii, nPhot2, GetEventCentralityBin()));
3437 
3438  fhEventMixBin->Fill(eventbin, GetEventWeight()) ;
3439 
3440  //---------------------------------
3441  // First loop on photons/clusters
3442  //---------------------------------
3443  for(Int_t i1 = 0; i1 < nPhot; i1++)
3444  {
3445  AliAODPWG4Particle * p1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i1)) ;
3446 
3447  // Select photons within a pT range
3448  if ( p1->Pt() < GetMinPt() || p1->Pt() > GetMaxPt() ) continue ;
3449 
3450  // Not sure why this line is here
3451  //if(fSameSM && GetModuleNumber(p1)!=module1) continue;
3452 
3453  //Get kinematics of cluster and (super) module of this cluster
3454  fPhotonMom1.SetPxPyPzE(p1->Px(),p1->Py(),p1->Pz(),p1->E());
3455  module1 = GetModuleNumber(p1);
3456 
3457  //---------------------------------
3458  // Second loop on other mixed event photons/clusters
3459  //---------------------------------
3460  for(Int_t i2 = 0; i2 < nPhot2; i2++)
3461  {
3462  AliAODPWG4Particle * p2 = (AliAODPWG4Particle*) (ev2->At(i2)) ;
3463 
3464  // Select photons within a pT range
3465  if ( p2->Pt() < GetMinPt() || p2->Pt() > GetMaxPt() ) continue ;
3466 
3467  // Get kinematics of second cluster and calculate those of the pair
3468  fPhotonMom2.SetPxPyPzE(p2->Px(),p2->Py(),p2->Pz(),p2->E());
3469  m = (fPhotonMom1+fPhotonMom2).M() ;
3470  Double_t pt = (fPhotonMom1 + fPhotonMom2).Pt();
3471  Double_t a = TMath::Abs(p1->E()-p2->E())/(p1->E()+p2->E()) ;
3472 
3473  // Check if opening angle is too large or too small compared to what is expected
3474  Double_t angle = fPhotonMom1.Angle(fPhotonMom2.Vect());
3476  {
3477  AliDebug(2,Form("Mix pair angle %f (deg) not in E %f window",RadToDeg(angle), (fPhotonMom1+fPhotonMom2).E()));
3478  continue;
3479  }
3480 
3481  if(fUseAngleCut && (angle < fAngleCut || angle > fAngleMaxCut))
3482  {
3483  AliDebug(2,Form("Mix pair cut %f < angle %f < cut %f (deg)",RadToDeg(fAngleCut),RadToDeg(angle),RadToDeg(fAngleMaxCut)));
3484  continue;
3485  }
3486 
3487  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));
3488 
3489  // In case we want only pairs in same (super) module, check their origin.
3490  module2 = GetModuleNumber(p2);
3491 
3492  //-------------------------------------------------------------------------------------------------
3493  // Fill module dependent histograms, put a cut on assymmetry on the first available cut in the array
3494  //-------------------------------------------------------------------------------------------------
3495  if(a < fAsymCuts[0] && fFillSMCombinations)
3496  {
3498  {
3499  if(module1==module2 && module1 >=0 && module1<fNModules)
3500  {
3501  fhMiMod[module1]->Fill(pt, m, GetEventWeight()) ;
3502  if(fFillAngleHisto) fhMixedOpeningAnglePerSM[module1]->Fill(pt, angle, GetEventWeight());
3503  }
3504 
3505  if(GetCalorimeter()==kEMCAL)
3506  {
3507  // Same sector
3508  Int_t j=0;
3509  for(Int_t i = 0; i < fNModules/2; i++)
3510  {
3511  j=2*i;
3512  if((module1==j && module2==j+1) || (module1==j+1 && module2==j)) fhMiSameSectorEMCALMod[i]->Fill(pt, m, GetEventWeight()) ;
3513  }
3514 
3515  // Same side
3516  for(Int_t i = 0; i < fNModules-2; i++)
3517  {
3518  if((module1==i && module2==i+2) || (module1==i+2 && module2==i)) fhMiSameSideEMCALMod[i]->Fill(pt, m, GetEventWeight());
3519  }
3520  } // EMCAL
3521  else
3522  { // PHOS
3523  if((module1==0 && module2==1) || (module1==1 && module2==0)) fhMiDiffPHOSMod[0]->Fill(pt, m, GetEventWeight()) ;
3524  if((module1==0 && module2==2) || (module1==2 && module2==0)) fhMiDiffPHOSMod[1]->Fill(pt, m, GetEventWeight()) ;
3525  if((module1==1 && module2==2) || (module1==2 && module2==1)) fhMiDiffPHOSMod[2]->Fill(pt, m, GetEventWeight()) ;
3526  } // PHOS
3527  }
3528  else
3529  {
3530  Float_t phi1 = GetPhi(fPhotonMom1.Phi());
3531  Float_t phi2 = GetPhi(fPhotonMom2.Phi());
3532  Bool_t etaside = 0;
3533  if( (p1->GetDetectorTag()==kEMCAL && fPhotonMom1.Eta() < 0)
3534  || (p2->GetDetectorTag()==kEMCAL && fPhotonMom2.Eta() < 0)) etaside = 1;
3535 
3536  if ( phi1 > DegToRad(260) && phi2 > DegToRad(260) && phi1 < DegToRad(280) && phi2 < DegToRad(280)) fhMiSameSectorDCALPHOSMod[0+etaside]->Fill(pt, m, GetEventWeight());
3537  else if ( phi1 > DegToRad(280) && phi2 > DegToRad(280) && phi1 < DegToRad(300) && phi2 < DegToRad(300)) fhMiSameSectorDCALPHOSMod[2+etaside]->Fill(pt, m, GetEventWeight());
3538  else if ( phi1 > DegToRad(300) && phi2 > DegToRad(300) && phi1 < DegToRad(320) && phi2 < DegToRad(320)) fhMiSameSectorDCALPHOSMod[4+etaside]->Fill(pt, m, GetEventWeight());
3539  else if ( (phi1 > DegToRad(260) && phi2 > DegToRad(280) && phi1 < DegToRad(280) && phi2 < DegToRad(300))
3540  || (phi1 > DegToRad(280) && phi2 > DegToRad(260) && phi1 < DegToRad(300) && phi2 < DegToRad(280))) fhMiDiffSectorDCALPHOSMod[0+etaside]->Fill(pt, m, GetEventWeight());
3541  else if ( (phi1 > DegToRad(280) && phi2 > DegToRad(300) && phi1 < DegToRad(300) && phi2 < DegToRad(320))
3542  || (phi1 > DegToRad(300) && phi2 > DegToRad(280) && phi1 < DegToRad(320) && phi2 < DegToRad(300))) fhMiDiffSectorDCALPHOSMod[2+etaside]->Fill(pt, m, GetEventWeight());
3543  else if ( (phi1 > DegToRad(260) && phi2 > DegToRad(300) && phi1 < DegToRad(280) && phi2 < DegToRad(320))
3544  || (phi1 > DegToRad(300) && phi2 > DegToRad(260) && phi1 < DegToRad(320) && phi2 < DegToRad(280))) fhMiDiffSectorDCALPHOSMod[4+etaside]->Fill(pt, m, GetEventWeight());
3545  else fhMiDiffSectorDCALPHOSMod[6+etaside]->Fill(pt, m, GetEventWeight());
3546  }
3547  }
3548 
3549  Bool_t ok = kTRUE;
3550  if(fSameSM)
3551  {
3553  {
3554  if(module1!=module2) ok=kFALSE;
3555  }
3556  else // PHOS and DCal in same sector
3557  {
3558  Float_t phi1 = GetPhi(fPhotonMom1.Phi());
3559  Float_t phi2 = GetPhi(fPhotonMom2.Phi());
3560  ok=kFALSE;
3561  if ( phi1 > DegToRad(260) && phi2 > DegToRad(260) && phi1 < DegToRad(280) && phi2 < DegToRad(280)) ok = kTRUE;
3562  else if ( phi1 > DegToRad(280) && phi2 > DegToRad(280) && phi1 < DegToRad(300) && phi2 < DegToRad(300)) ok = kTRUE;
3563  else if ( phi1 > DegToRad(300) && phi2 > DegToRad(300) && phi1 < DegToRad(320) && phi2 < DegToRad(320)) ok = kTRUE;
3564  }
3565  } // Pair only in same SM
3566 
3567  if(!ok) continue ;
3568 
3569  //
3570  // Do the final histograms with the selected clusters
3571  //
3572 
3573  // Check if one of the clusters comes from a conversion
3574  if(fCheckConversion)
3575  {
3576  if (p1->IsTagged() && p2->IsTagged()) fhMiConv2->Fill(pt, m, GetEventWeight());
3577  else if(p1->IsTagged() || p2->IsTagged()) fhMiConv ->Fill(pt, m, GetEventWeight());
3578  }
3579 
3580  // Fill histograms for different bad channel distance, centrality, assymmetry cut and pid bit
3581  for(Int_t ipid=0; ipid<fNPIDBits; ipid++)
3582  {
3583  if((p1->IsPIDOK(ipid,AliCaloPID::kPhoton)) && (p2->IsPIDOK(ipid,AliCaloPID::kPhoton)))
3584  {
3585  for(Int_t iasym=0; iasym < fNAsymCuts; iasym++)
3586  {
3587  if(a < fAsymCuts[iasym])
3588  {
3589  Int_t index = ((curCentrBin*fNPIDBits)+ipid)*fNAsymCuts + iasym;
3590 
3591  if(index < 0 || index >= ncentr*fNPIDBits*fNAsymCuts) continue ;
3592 
3593  fhMi1[index]->Fill(pt, m, GetEventWeight()) ;
3594  if(fMakeInvPtPlots)fhMiInvPt1[index]->Fill(pt, m, 1./pt * GetEventWeight()) ;
3595 
3596  if(fFillBadDistHisto)
3597  {
3598  if(p1->DistToBad()>0 && p2->DistToBad()>0)
3599  {
3600  fhMi2[index]->Fill(pt, m, GetEventWeight()) ;
3601  if(fMakeInvPtPlots)fhMiInvPt2[index]->Fill(pt, m, 1./pt * GetEventWeight()) ;
3602 
3603  if(p1->DistToBad()>1 && p2->DistToBad()>1)
3604  {
3605  fhMi3[index]->Fill(pt, m, GetEventWeight()) ;
3606  if(fMakeInvPtPlots)fhMiInvPt3[index]->Fill(pt, m, 1./pt * GetEventWeight()) ;
3607  }
3608  }
3609  }// Fill bad dist histo
3610 
3611  //-----------------------
3612  // Multi cuts analysis
3613  //-----------------------
3614  Int_t ncell1 = p1->GetNCells();
3615  Int_t ncell2 = p1->GetNCells();
3616 
3617  if(fMultiCutAna)
3618  {
3619  // Several pt,ncell and asymmetry cuts
3620  for(Int_t ipt=0; ipt<fNPtCuts; ipt++)
3621  {
3622  for(Int_t icell=0; icell<fNCellNCuts; icell++)
3623  {
3624  for(Int_t iasym=0; iasym<fNAsymCuts; iasym++)
3625  {
3626  Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
3627 
3628  if(p1->Pt() > fPtCuts[ipt] && p2->Pt() > fPtCuts[ipt] &&
3629  p1->Pt() < fPtCutsMax[ipt] && p2->Pt() < fPtCutsMax[ipt] &&
3630  a < fAsymCuts[iasym] &&
3631  ncell1 >= fCellNCuts[icell] && ncell2 >= fCellNCuts[icell]
3632  )
3633  {
3634  //printf("MI ipt %d, iasym%d, icell %d, index %d \n",ipt, iasym, icell, index);
3635  //printf("\t %p, %p\n",fhMiPtNCellAsymCuts[index],fhMiPtNCellAsymCutsOpAngle[index]);
3636 
3637  fhMiPtNCellAsymCuts[index]->Fill(pt, m, GetEventWeight()) ;
3638  if(fFillAngleHisto) fhMiPtNCellAsymCutsOpAngle[index]->Fill(pt, angle, GetEventWeight()) ;
3639 
3640  //printf("ipt %d, icell%d, iasym %d, name %s\n",ipt, icell, iasym, fhRePtNCellAsymCuts[((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym]->GetName());
3641  }
3642  }// pid bit cut loop
3643  }// icell loop
3644  }// pt cut loop
3645  } // Multi cut ana
3646 
3647  //
3648  // Fill histograms with opening angle
3649  if(fFillAngleHisto)
3650  {
3651  fhMixedOpeningAngle ->Fill(pt, angle, GetEventWeight());
3652  fhMixedCosOpeningAngle->Fill(pt, TMath::Cos(angle), GetEventWeight());
3653  }
3654 
3655  //
3656  // Fill histograms for different opening angle bins
3658  {
3659  Int_t angleBin = -1;
3660  for(Int_t ibin = 0; ibin < fNAngleCutBins; ibin++)
3661  {
3662  if(angle >= fAngleCutBinsArray[ibin] &&
3663  angle < fAngleCutBinsArray[ibin+1]) angleBin = ibin;
3664  }
3665 
3666  if( angleBin >= 0 && angleBin < fNAngleCutBins)
3667  {
3668  Float_t e1 = fPhotonMom1.E();
3669  Float_t e2 = fPhotonMom2.E();
3670 
3671  Float_t t1 = p1->GetTime();
3672  Float_t t2 = p2->GetTime();
3673 
3674  Int_t nc1 = ncell1;
3675  Int_t nc2 = ncell2;
3676 
3677  Float_t eta1 = fPhotonMom1.Eta();
3678  Float_t eta2 = fPhotonMom2.Eta();
3679 
3680  Float_t phi1 = GetPhi(fPhotonMom1.Phi());
3681  Float_t phi2 = GetPhi(fPhotonMom2.Phi());
3682 
3683  Int_t mod1 = module1;
3684  Int_t mod2 = module2;
3685 
3686  // // Recover original cluster
3687  // Int_t iclus1 = -1, iclus2 = -1 ;
3688  // AliVCluster * cluster1 = FindCluster(GetEMCALClusters(),p1->GetCaloLabel(0),iclus1);
3689  // AliVCluster * cluster2 = FindCluster(GetEMCALClusters(),p2->GetCaloLabel(0),iclus2);
3690  //
3691  // Float_t maxCellFraction1 = 0, maxCellFraction2 = 0;
3692  // Int_t absIdMax1 = GetCaloUtils()->GetMaxEnergyCell(GetEMCALCells(),cluster1,maxCellFraction1);
3693  // Int_t absIdMax2 = GetCaloUtils()->GetMaxEnergyCell(GetEMCALCells(),cluster2,maxCellFraction2);
3694 
3695  if(e2 > e1)
3696  {
3697  e1 = fPhotonMom2.E();
3698  e2 = fPhotonMom1.E();
3699 
3700  t1 = p2->GetTime();
3701  t2 = p1->GetTime();
3702 
3703  nc1 = ncell2;
3704  nc2 = ncell1;
3705 
3706  eta1 = fPhotonMom2.Eta();
3707  eta2 = fPhotonMom1.Eta();
3708 
3709  phi1 = GetPhi(fPhotonMom2.Phi());
3710  phi2 = GetPhi(fPhotonMom1.Phi());
3711 
3712  mod1 = module2;
3713  mod2 = module1;
3714 
3715  // Int_t tmp = absIdMax2;
3716  // absIdMax2 = absIdMax1;
3717  // absIdMax1 = tmp;
3718  }
3719 
3720  fhMiOpAngleBinMinClusterEPerSM[angleBin]->Fill(e2,mod2,GetEventWeight()) ;
3721  fhMiOpAngleBinMaxClusterEPerSM[angleBin]->Fill(e1,mod1,GetEventWeight()) ;
3722 
3723  fhMiOpAngleBinMinClusterTimePerSM[angleBin]->Fill(t2,mod2,GetEventWeight()) ;
3724  fhMiOpAngleBinMaxClusterTimePerSM[angleBin]->Fill(t1,mod1,GetEventWeight()) ;
3725 
3726  fhMiOpAngleBinMinClusterNCellPerSM[angleBin]->Fill(nc2,mod2,GetEventWeight()) ;
3727  fhMiOpAngleBinMaxClusterNCellPerSM[angleBin]->Fill(nc1,mod1,GetEventWeight()) ;
3728 
3729  fhMiOpAngleBinPairClusterMass[angleBin]->Fill(pt,m,GetEventWeight()) ;
3730  if(mod2 == mod1) fhMiOpAngleBinPairClusterMassPerSM[angleBin]->Fill(m,mod1,GetEventWeight()) ;
3731 
3732  if(e1 > 0.01) fhMiOpAngleBinPairClusterRatioPerSM[angleBin]->Fill(e2/e1,mod1,GetEventWeight()) ;
3733 
3734  fhMiOpAngleBinMinClusterEtaPhi[angleBin]->Fill(eta2,phi2,GetEventWeight()) ;
3735  fhMiOpAngleBinMaxClusterEtaPhi[angleBin]->Fill(eta1,phi1,GetEventWeight()) ;
3736 
3737  // Int_t icol1 = -1, icol2 = -1, icolAbs1 = -1, icolAbs2 = -1;
3738  // Int_t irow1 = -1, irow2 = -1, irowAbs1 = -1, irowAbs2 = -1;
3739  // Int_t iRCU1 = -1, iRCU2 = -1;
3740  // GetModuleNumberCellIndexesAbsCaloMap(absIdMax1,GetCalorimeter(), icol1, irow1, iRCU1, icolAbs1, irowAbs1);
3741  // GetModuleNumberCellIndexesAbsCaloMap(absIdMax2,GetCalorimeter(), icol2, irow2, iRCU1, icolAbs2, irowAbs2);
3742  //
3743  // fhMiOpAngleBinPairClusterAbsIdMaxCell[angleBin]->Fill(absIdMax1,absIdMax2,GetEventWeight());
3744  //
3745  // fhMiColRowClusterMinOpAngleBin[angleBin]->Fill(icolAbs2,irowAbs2,GetEventWeight()) ;
3746  // fhMiOpAngleBinMaxClusterColRow[angleBin]->Fill(icolAbs1,irowAbs1,GetEventWeight()) ;
3747  }
3748  }
3749 
3750  // Fill histograms with pair assymmetry
3752  {
3753  fhMiPtAsym->Fill(pt, a, GetEventWeight());
3754  if(m > 0.10 && m < 0.17) fhMiPtAsymPi0->Fill(pt, a, GetEventWeight());
3755  if(m > 0.45 && m < 0.65) fhMiPtAsymEta->Fill(pt, a, GetEventWeight());
3756  }
3757 
3758  // Check cell time content in cluster
3760  {
3761  if ( p1->GetFiducialArea() == 0 && p2->GetFiducialArea() == 0 )
3763 
3764  else if ( p1->GetFiducialArea() != 0 && p2->GetFiducialArea() != 0 )
3766  }
3767 
3768  }//Asymmetry cut
3769  }// Asymmetry loop
3770  }//PID cut
3771  }//loop for histograms
3772 
3773  }// second cluster loop
3774  }//first cluster loop
3775  }//loop on mixed events
3776 
3777  //--------------------------------------------------------
3778  // Add the current event to the list of events for mixing
3779  //--------------------------------------------------------
3780 
3781  //TClonesArray *currentEvent = new TClonesArray(*GetInputAODBranch());
3782  TClonesArray *currentEvent = new TClonesArray(*secondLoopInputData);
3783 
3784  // Add current event to buffer and Remove redundant events
3785  if( currentEvent->GetEntriesFast() > 0 )
3786  {
3787  evMixList->AddFirst(currentEvent) ;
3788  currentEvent=0 ; //Now list of particles belongs to buffer and it will be deleted with buffer
3789  if( evMixList->GetSize() >= GetNMaxEvMix() )
3790  {
3791  TClonesArray * tmp = (TClonesArray*) (evMixList->Last()) ;
3792  evMixList->RemoveLast() ;
3793  delete tmp ;
3794  }
3795  }
3796  else
3797  { // empty event
3798  delete currentEvent ;
3799  currentEvent=0 ;
3800  }
3801  }// DoOwnMix
3802 
3803  AliDebug(1,"End fill histograms");
3804 }
3805 
3806 //________________________________________________________________________
3810 //________________________________________________________________________
3811 Int_t AliAnaPi0::GetEventIndex(AliAODPWG4Particle * part, Double_t * vert)
3812 {
3813  Int_t evtIndex = -1 ;
3814 
3815  if(GetReader()->GetDataType()!=AliCaloTrackReader::kMC)
3816  {
3817  if (GetMixedEvent())
3818  {
3819  evtIndex = GetMixedEvent()->EventIndexForCaloCluster(part->GetCaloLabel(0)) ;
3820  GetVertex(vert,evtIndex);
3821 
3822  if(TMath::Abs(vert[2])> GetZvertexCut())
3823  evtIndex = -2 ; //Event can not be used (vertex, centrality,... cuts not fulfilled)
3824  }
3825  else
3826  {
3827  // Single event
3828  GetVertex(vert);
3829 
3830  if(TMath::Abs(vert[2])> GetZvertexCut())
3831  evtIndex = -1 ; //Event can not be used (vertex, centrality,... cuts not fulfilled)
3832  else
3833  evtIndex = 0 ;
3834  }
3835  } // No MC reader
3836  else
3837  {
3838  evtIndex = 0;
3839  vert[0] = 0. ;
3840  vert[1] = 0. ;
3841  vert[2] = 0. ;
3842  }
3843 
3844  return evtIndex ;
3845 }
3846 
TH2F * fhPrimEtaY
! Rapidity distribution of primary particles vs pT
Definition: AliAnaPi0.h:423
Float_t GetHistoPtMax() const
TH2F ** fhRe3
[GetNCentrBin()*fNPIDBits*fNAsymCuts]
Definition: AliAnaPi0.h:303
TH2F * fhPrimEtaAccPtCentrality
! primary eta with accepted daughters reconstructed centrality vs pT
Definition: AliAnaPi0.h:436
TH2F * fhMiOpAngleBinMinClusterEPerSM[10]
! energy of lowest energy cluster in pair, depending on opening angle cut, y axis is SM number...
Definition: AliAnaPi0.h:525
TH2F * fhReOpAngleBinMaxClusterEPerSM[10]
! energy of highest energy cluster in pair, depending on opening angle cut, y axis is SM number ...
Definition: AliAnaPi0.h:508
Int_t pdg
TH2F * fhPrimPi0Y
! Rapidity distribution of primary particles vs pT
Definition: AliAnaPi0.h:399
TLorentzVector fPhotonMom1Boost
! Photon cluster momentum, temporary array
Definition: AliAnaPi0.h:226
TH2F * fhMiOpAngleBinPairClusterRatioPerSM[10]
! lowest/highest energy cluster in pair, depending on opening angle cut, y axis is SM number...
Definition: AliAnaPi0.h:531
Int_t fNModules
[GetNCentrBin()*GetNZvertBin()*GetNRPBin()]
Definition: AliAnaPi0.h:181
TH2F * fhPrimEtaYetaYcut
! PseudoRapidity distribution of primary particles vs pT, Y<1
Definition: AliAnaPi0.h:426
TH1F * fhPrimPi0AccPt
! Spectrum of primary with accepted daughters
Definition: AliAnaPi0.h:396
Float_t GetHistoPtMin() const
TH2F * fhPrimEtaAccPhi
! Azimutal distribution of primary with accepted daughters vs pT
Definition: AliAnaPi0.h:429
TLorentzVector fPi0Mom
! Pi0 cluster momentum, temporary array
Definition: AliAnaPi0.h:228
double Double_t
Definition: External.C:58
TH2F * fhReMCFromMixConversion
! Invariant mass of 2 clusters one from conversion and the other not
Definition: AliAnaPi0.h:487
Int_t fPIDBits[10]
Array with different PID bits.
Definition: AliAnaPi0.h:200
TH2F * fhMiPtAsymPi0
! Mix two-photon pt vs asymmetry, close to pi0 mass
Definition: AliAnaPi0.h:368
virtual void AddToHistogramsName(TString add)
virtual AliVCaloCells * GetEMCALCells() const
TH2F * fhReConv
[8]
Definition: AliAnaPi0.h:282
TH2F * fhReOpAngleBinPairClusterMassMCTrueEta[10]
! cluster pair mass vs pT, depending on opening angle cut, true eta decay pairs from MC ...
Definition: AliAnaPi0.h:519
Definition: External.C:236
Bool_t fFillArmenterosThetaStar
Fill armenteros histograms.
Definition: AliAnaPi0.h:215
Int_t fNPIDBits
Number of possible PID bit combinations.
Definition: AliAnaPi0.h:199
TLorentzVector fPhotonMom2
! Photon cluster momentum, temporary array
Definition: AliAnaPi0.h:227
TH2F * fhReOpAngleBinMaxClusterColRow[10]
! Column and row location of main cell of highest energy cluster in pair, depending on opening angle ...
Definition: AliAnaPi0.h:506
const char * title
Definition: MakeQAPdf.C:26
TH2F ** fhMiDiffSectorDCALPHOSMod
[6]
Definition: AliAnaPi0.h:278
Int_t fCellNCuts[10]
Array with different cell number cluster cuts.
Definition: AliAnaPi0.h:198
TH2F * fhMCOrgDeltaEta[13]
! Delta Eta vs pt of real pairs, check common origin of pair
Definition: AliAnaPi0.h:452
TH2F * fhReSS[3]
[fNPtCuts*fNAsymCuts*fNCellNCuts]
Definition: AliAnaPi0.h:360
TH2F * fhReOpAngleBinMaxClusterNCellPerSM[10]
! N cells of highest energy cluster in pair, depending on opening angle cut, y axis is SM number ...
Definition: AliAnaPi0.h:512
TH2F * fhPrimEtaCosOpeningAngle
! Cosinus of opening angle of pair version pair energy, eta primaries
Definition: AliAnaPi0.h:433
TH2F * fhPrimPi0PtStatus
! Spectrum of generated pi0 vs pi0 status
Definition: AliAnaPi0.h:444
Bool_t fCheckConversion
Fill histograms with tagged photons as conversion.
Definition: AliAnaPi0.h:209
virtual Int_t GetModuleNumberCellIndexesAbsCaloMap(Int_t absId, Int_t calo, Int_t &icol, Int_t &irow, Int_t &iRCU, Int_t &icolAbs, Int_t &irowAbs) const
TH2F * fhPtBinClusterEtaPhi[10]
! Eta-Phi location of cluster in different energy bins.
Definition: AliAnaPi0.h:536
virtual Float_t GetPairTimeCut() const
Time cut in ns.
virtual void GetVertex(Double_t vertex[3]) const
TH1F * fhPrimPi0E
! Spectrum of Primary
Definition: AliAnaPi0.h:393
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:184
TH2F ** fhMiSameSectorEMCALMod
[fNModules-2]
Definition: AliAnaPi0.h:269
TH2F * fhMixedOpeningAnglePerSM[20]
! Opening angle of pair versus pair energy, per SM
Definition: AliAnaPi0.h:388
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:375
TH2F * fhMiSecondaryCellOutTimeWindow
! Combine clusters when at least one significant cells in cluster has t > 50 ns, different events ...
Definition: AliAnaPi0.h:500
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:222
Float_t DegToRad(Float_t deg) const
TH2F * fhPrimPi0PtCentrality
! primary pi0 reconstructed centrality vs pT
Definition: AliAnaPi0.h:410
TH2F * fhPrimEtaPtCentrality
! primary eta reconstructed centrality vs pT
Definition: AliAnaPi0.h:434
virtual Float_t GetZvertexCut() const
Maximal number of events for mixin.
virtual Double_t GetEventPlaneAngle() const
TH2F * fhReOpAngleBinPairClusterMassPerSM[10]
! cluster pair mass, depending on opening angle cut, y axis is SM number
Definition: AliAnaPi0.h:515
TH2F * fhReMCFromConversion
! Invariant mass of 2 clusters originated in conversions
Definition: AliAnaPi0.h:485
TH2F ** fhReSameSideEMCALMod
[fNModules]
Definition: AliAnaPi0.h:248
TH1F * fhPrimPi0AccPtOpAngCuts[10]
! Spectrum of primary with accepted daughters, different opening angles
Definition: AliAnaPi0.h:398
TH1F * fhPrimPi0AccPtPhotonCuts
! Spectrum of primary with accepted daughters, photon pt or angle cuts
Definition: AliAnaPi0.h:397
TH2F * fhPrimPi0OpeningAngle
! Opening angle of pair versus pair energy, primaries
Definition: AliAnaPi0.h:406
virtual AliNeutralMesonSelection * GetNeutralMesonSelection()
TH2F * fhRealOpeningAngle
! Opening angle of pair versus pair energy
Definition: AliAnaPi0.h:382
TH2F * fhMCOrgMass[13]
! Mass vs pt of real pairs, check common origin of pair
Definition: AliAnaPi0.h:450
TH2F * fhMiOpAngleBinMinClusterEtaPhi[10]
! Eta-Phi location of lowest energy cluster in pair, depending on opening angle cut, mixed event
Definition: AliAnaPi0.h:521
TH2F * fhMiOpAngleBinPairClusterMass[10]
! cluster pair mass vs pT, depending on opening angle cut, mixed event
Definition: AliAnaPi0.h:532
Int_t GetHistoMassBins() const
Int_t GetHistoPhiBins() const
TH2F * fhRePtAsymPi0
! REAL two-photon pt vs asymmetry, close to pi0 mass
Definition: AliAnaPi0.h:365
TLorentzVector fPhotonMom1
! Photon cluster momentum, temporary array
Definition: AliAnaPi0.h:225
TH2F * fhRealCosOpeningAngle
! Cosinus of opening angle of pair version pair energy
Definition: AliAnaPi0.h:383
TH2F ** fhRe1
REAL two-photon invariant mass distribution for different centralities and Asymmetry.
Definition: AliAnaPi0.h:288
Float_t GetHistoMassMin() const
TH2F * fhPrimPi0ProdVertex
! Spectrum of primary pi0 vs production vertex
Definition: AliAnaPi0.h:482
TH2F ** fhReDiffPHOSMod
[fNModules/2]
Definition: AliAnaPi0.h:254
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:283
TH2F * fhReOpAngleBinPairClusterMass[10]
! cluster pair mass vs pT, depending on opening angle cut
Definition: AliAnaPi0.h:514
Float_t fAsymCuts[10]
Array with different assymetry cuts.
Definition: AliAnaPi0.h:196
TH2F * fhMCOrgAsym[13]
! Asymmetry vs pt of real pairs, check common origin of pair
Definition: AliAnaPi0.h:451
TH2F ** fhRePtNCellAsymCuts
[GetNCentrBin()*fNPIDBits*fNAsymCuts]
Definition: AliAnaPi0.h:336
TH2F * fhPrimPi0PtOrigin
! Spectrum of generated pi0 vs mother
Definition: AliAnaPi0.h:441
Float_t GetHistoDiffTimeMin() const
Float_t fPtCuts[11]
Array with different pt cuts, minimum.
Definition: AliAnaPi0.h:193
TH2F * fhRealOpeningAnglePerSM[20]
! Opening angle of pair versus pair energy, per SM
Definition: AliAnaPi0.h:387
TH2F * fhMiOpAngleBinMaxClusterNCellPerSM[10]
! N cells of highest energy cluster in pair, depending on opening angle cut, y axis is SM number...
Definition: AliAnaPi0.h:530
Bool_t IsMCParticleInCalorimeterAcceptance(Int_t calo, TParticle *particle)
TH2F ** fhMCPi0MassPtRec
Real pi0 pairs, reconstructed mass vs reconstructed pt of original pair.
Definition: AliAnaPi0.h:458
Bool_t fFillOpAngleCutHisto
Fill histograms depending on opening angle of pair.
Definition: AliAnaPi0.h:218
TH2F * fhMiPtAsym
! Mix two-photon pt vs asymmetry
Definition: AliAnaPi0.h:367
TH2F * fhPrimPi0AccPtCentrality
! primary pi0 with accepted daughters reconstructed centrality vs pT
Definition: AliAnaPi0.h:412
Float_t GetHistoPhiMin() const
Int_t fNAsymCuts
Number of assymmetry cuts.
Definition: AliAnaPi0.h:195
Float_t GetHistoDiffTimeMax() const
TH2F * fhMixedCosOpeningAngle
! Cosinus of opening angle of pair version pair energy
Definition: AliAnaPi0.h:385
TVector3 fProdVertex
! Production vertex, temporary array
Definition: AliAnaPi0.h:229
TH2F ** fhReInvPt2
[GetNCentrBin()*fNPIDBits*fNAsymCuts]
Definition: AliAnaPi0.h:319
TH2F ** fhMiDiffPHOSMod
[fNModules/2]
Definition: AliAnaPi0.h:272
TH2F ** fhMi2
[GetNCentrBin()*fNPIDBits*fNAsymCuts]
Definition: AliAnaPi0.h:299
Float_t GetHistoOpeningAngleMax() const
TList * GetCreateOutputObjects()
Definition: AliAnaPi0.cxx:306
TH2F ** fhReInvPt3
[GetNCentrBin()*fNPIDBits*fNAsymCuts]
Definition: AliAnaPi0.h:327
Bool_t fFillOriginHisto
Fill histograms depending on their origin.
Definition: AliAnaPi0.h:214
TH2F * fhPrimPi0Phi
! Azimutal distribution of primary particles vs pT
Definition: AliAnaPi0.h:404
TH2F * fhPrimEtaOpeningAnglePhotonCuts
! Opening angle of pair versus pair energy, eta primaries, photon pT cuts
Definition: AliAnaPi0.h:431
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, Double_t angle)
Definition: AliAnaPi0.cxx:2483
Bool_t IsAngleInWindow(Float_t angle, Float_t e) const
TH2F ** fhReSameSectorDCALPHOSMod
[fNModules]
Definition: AliAnaPi0.h:257
TH2F ** fhMiPtNCellAsymCutsOpAngle
[fNPtCuts*fNAsymCuts*fNCellNCuts]
Definition: AliAnaPi0.h:348
const Double_t etamin
Float_t GetHistoAsymmetryMax() const
TH2F * fhMiOpAngleBinMaxClusterEPerSM[10]
! energy of highest energy cluster in pair, depending on opening angle cut, y axis is SM number...
Definition: AliAnaPi0.h:526
Float_t GetHistoMassMax() const
Base class for CaloTrackCorr analysis algorithms.
Bool_t fFillAngleHisto
Fill histograms with pair opening angle.
Definition: AliAnaPi0.h:212
virtual TString GetCalorimeterString() const
Bool_t fFillBadDistHisto
Do plots for different distances to bad channels.
Definition: AliAnaPi0.h:210
Float_t fAngleCutBinsArray[11]
Array with angle cut bins.
Definition: AliAnaPi0.h:203
virtual Bool_t IsRealCaloAcceptanceOn() const
virtual AliFiducialCut * GetFiducialCut()
int Int_t
Definition: External.C:63
virtual TClonesArray * GetInputAODBranch() const
TH2F ** fhMCEtaPtTruePtRec
[fNPtCuts*fNAsymCuts*fNCellNCuts]
Definition: AliAnaPi0.h:473
Bool_t fMultiCutAna
Do analysis with several or fixed cut.
Definition: AliAnaPi0.h:189
Int_t GetHistoNOpeningAngleBins() const
virtual TClonesArray * GetAODMCParticles() const
Definition: External.C:204
virtual AliHistogramRanges * GetHistogramRanges()
TH2F ** fhReMod
REAL two-photon invariant mass distribution for different calorimeter modules.
Definition: AliAnaPi0.h:245
Bool_t ReadStack() const
TH1F * fhCentralityNoPair
! Histogram with centrality bins with no pair
Definition: AliAnaPi0.h:376
Bool_t fFillAsymmetryHisto
Fill histograms with asymmetry vs pt.
Definition: AliAnaPi0.h:213
Int_t GetHistoDiffTimeBins() const
TH2F * fhReOpAngleBinMinClusterTimePerSM[10]
! time of lowest energy cluster in pair, depending on opening angle cut, y axis is SM number ...
Definition: AliAnaPi0.h:509
virtual TList * GetAODBranchList() const
TH2F * fhReOpAngleBinMinClusterNCellPerSM[10]
! N cells of lowest energy cluster in pair, depending on opening angle cut, y axis is SM number ...
Definition: AliAnaPi0.h:511
Int_t GetEventIndex(AliAODPWG4Particle *part, Double_t *vert)
Definition: AliAnaPi0.cxx:3811
float Float_t
Definition: External.C:68
void FillAcceptanceHistograms()
Fill acceptance histograms if MC data is available.
Definition: AliAnaPi0.cxx:1968
TH2F ** fhReSameSectorEMCALMod
[fNModules-2]
Definition: AliAnaPi0.h:251
TH2F ** fhMCEtaMassPtTrue
[fNPtCuts*fNAsymCuts*fNCellNCuts]
Definition: AliAnaPi0.h:470
TH2F ** fhReInvPt1
[GetNCentrBin()*fNPIDBits*fNAsymCuts]
Definition: AliAnaPi0.h:312
Int_t fNCellNCuts
Number of cuts with number of cells in cluster.
Definition: AliAnaPi0.h:197
TH2F ** fhMiInvPt1
[GetNCentrBin()*fNPIDBits*fNAsymCuts]
Definition: AliAnaPi0.h:315
Int_t GetHistoAsymmetryBins() const
const Double_t ptmax
TH2F * fhPrimEtaAccYeta
! PseudoRapidity distribution of primary with accepted daughters vs pT
Definition: AliAnaPi0.h:427
Float_t fAngleMaxCut
Select pairs with opening angle smaller than a threshold.
Definition: AliAnaPi0.h:186
TH2F * fhReOpAngleBinPairClusterMassMCTruePi0[10]
! cluster pair mass vs pT, depending on opening angle cut, true pi0 decay pairs from MC ...
Definition: AliAnaPi0.h:518
TH1I * fhEventMixBin
! Number of mixed pairs in a particular bin (cen,vz,rp)
Definition: AliAnaPi0.h:374
void FillArmenterosThetaStar(Int_t pdg)
Fill armenteros plots.
Definition: AliAnaPi0.cxx:2420
virtual AliEMCALGeometry * GetEMCALGeometry() const
TH2F * fhPrimPi0CosOpeningAngle
! Cosinus of opening angle of pair version pair energy, pi0 primaries
Definition: AliAnaPi0.h:409
TH2F * fhMCPi0PtOrigin
[fNPtCuts*fNAsymCuts*fNCellNCuts]
Definition: AliAnaPi0.h:475
TH2F ** fhMiSameSideEMCALMod
[fNModules]
Definition: AliAnaPi0.h:266
TH1F * fhPrimEtaAccPtOpAngCuts[10]
! Spectrum of primary with accepted daughters, different opening angles
Definition: AliAnaPi0.h:422
Bool_t fFillOnlyMCAcceptanceHisto
Do analysis only of MC kinematics input.
Definition: AliAnaPi0.h:216
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:206
virtual AliAODEvent * GetOutputEvent() const
TH2F * fhPrimEtaPtEventPlane
! primary eta reconstructed event plane vs pT
Definition: AliAnaPi0.h:435
TH1F * fhPrimPi0AccE
! Spectrum of primary with accepted daughters
Definition: AliAnaPi0.h:395
Int_t fNPtCuts
Number of pt cuts.
Definition: AliAnaPi0.h:192
TH1F * fhPrimEtaE
! Spectrum of Primary
Definition: AliAnaPi0.h:417
TH2F * fhPrimEtaAccPtEventPlane
! primary eta with accepted daughters reconstructed event plane vs pT
Definition: AliAnaPi0.h:437
virtual AliCalorimeterUtils * GetCaloUtils() const
TH2F * fhPrimNotResonancePi0PtOrigin
! Spectrum of generated pi0 vs mother
Definition: AliAnaPi0.h:443
TH2F * fhPrimEtaAccY
! Rapidity distribution of primary with accepted daughters vs pT
Definition: AliAnaPi0.h:424
TH1F * fhPrimEtaPt
! Spectrum of Primary
Definition: AliAnaPi0.h:418
Bool_t fMultiCutAnaSim
Do analysis with several or fixed cut, in the simulation related part.
Definition: AliAnaPi0.h:190
TH1F * fhPrimEtaAccPt
! Spectrum of primary with accepted daughters
Definition: AliAnaPi0.h:420
TH2F * fhArmPrimPi0[4]
! Armenteros plots for primary pi0 in 6 energy bins
Definition: AliAnaPi0.h:489
void Print(const Option_t *opt) const
Print some relevant parameters set for the analysis.
Definition: AliAnaPi0.cxx:1922
TH2F * fhPrimPi0AccPtEventPlane
! primary pi0 with accepted daughters reconstructed event plane vs pT
Definition: AliAnaPi0.h:413
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:430
TH2F * fhReOpAngleBinMinClusterEtaPhi[10]
! Eta-Phi location of lowest energy cluster in pair, depending on opening angle cut ...
Definition: AliAnaPi0.h:503
TH2F ** fhMiPtNCellAsymCuts
[fNPtCuts*fNAsymCuts*fNCellNCuts]
Definition: AliAnaPi0.h:339
TH2F * fhPrimPi0OpeningAngleAsym
! Opening angle of pair versus pair E asymmetry, pi0 primaries
Definition: AliAnaPi0.h:408
TH2F * fhPrimPi0AccY
! Rapidity distribution of primary with accepted daughters vs pT
Definition: AliAnaPi0.h:400
TH2F * fhPrimPi0AccPhi
! Azimutal distribution of primary with accepted daughters vs pT
Definition: AliAnaPi0.h:405
Float_t fAngleCut
Select pairs with opening angle larger than a threshold.
Definition: AliAnaPi0.h:185
virtual Double_t GetEventWeight() const
TString fOtherDetectorInputName
String with name of extra detector data.
Definition: AliAnaPi0.h:223
TH2F * fhPrimPi0AccYeta
! PseudoRapidity distribution of primary with accepted daughters vs pT
Definition: AliAnaPi0.h:403
virtual Int_t GetNCentrBin() const
Number of bins in reaction plain.
TH2F ** fhRePtNCellAsymCutsSM[20]
[fNPtCuts*fNAsymCuts*fNCellNCuts]
Definition: AliAnaPi0.h:342
TH2F ** fhMi1
[GetNCentrBin()*fNPIDBits*fNAsymCuts]
Definition: AliAnaPi0.h:291
TH1F * fhPrimEtaAccPtPhotonCuts
! Spectrum of primary with accepted daughters, photon pt or angle cuts
Definition: AliAnaPi0.h:421
TH2F * fhCosThStarPrimEta
! cos(theta*) plots vs E for primary eta, same as asymmetry ...
Definition: AliAnaPi0.h:492
TH2F * fhReOpAngleBinMaxClusterEtaPhi[10]
! Eta-Phi location of highest energy cluster in pair, depending on opening angle cut ...
Definition: AliAnaPi0.h:504
TObjString * GetAnalysisCuts()
Save parameters used for analysis.
Definition: AliAnaPi0.cxx:258
void MakeAnalysisFillHistograms()
Definition: AliAnaPi0.cxx:2785
Float_t GetHistoEtaMin() const
Bool_t fFillSMCombinations
Fill histograms with different cluster pairs in SM combinations.
Definition: AliAnaPi0.h:208
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:486
TH2F * fhReOpAngleBinPairClusterRatioPerSM[10]
! lowest/highest energy cluster in pair, depending on opening angle cut, y axis is SM number ...
Definition: AliAnaPi0.h:513
TH2F * fhPrimPi0YetaYcut
! PseudoRapidity distribution of primary particles vs pT, Y<1
Definition: AliAnaPi0.h:402
TH2F ** fhRePtNCellAsymCutsSMOpAngle[20]
[fNPtCuts*fNAsymCuts*fNCellNCuts]
Definition: AliAnaPi0.h:351
TH2F * fhReOpAngleBinMaxClusterTimePerSM[10]
! time of highest energy cluster in pair, depending on opening angle cut, y axis is SM number ...
Definition: AliAnaPi0.h:510
virtual Int_t GetNMaxEvMix() const
Number of bins in track multiplicity.
TH2F * fhPrimPi0Yeta
! PseudoRapidity distribution of primary particles vs pT
Definition: AliAnaPi0.h:401
TH2F * fhReOpAngleBinMinClusterEPerSM[10]
! energy of lowest energy cluster in pair, depending on opening angle cut, y axis is SM number ...
Definition: AliAnaPi0.h:507
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:477
TH2F * fhRePtAsym
! REAL two-photon pt vs asymmetry
Definition: AliAnaPi0.h:364
TH2F * fhMiOpAngleBinMaxClusterTimePerSM[10]
! time of highest energy cluster in pair, depending on opening angle cut, y axis is SM number...
Definition: AliAnaPi0.h:528
Float_t GetHistoEtaMax() const
TH2F * fhEPairDiffTime
! E pair vs Pair of clusters time difference vs E
Definition: AliAnaPi0.h:494
void InitParameters()
Definition: AliAnaPi0.cxx:210
TH2F * fhMCEtaPtOrigin
! Mass of reconstructed eta pairs in calorimeter vs mother origin ID.
Definition: AliAnaPi0.h:476
TH2F * fhMCPi0PtStatus
! Mass of reconstructed pi0 pairs in calorimeter vs mother status.
Definition: AliAnaPi0.h:478
Int_t GetHistoPtBins() const
TH2F * fhMiConv2
! MIXED two-photon invariant mass distribution both pair photons recombined from 2 clusters with smal...
Definition: AliAnaPi0.h:285
TH2F ** fhRePtNCellAsymCutsOpAngle
[fNPtCuts*fNAsymCuts*fNCellNCuts]
Definition: AliAnaPi0.h:345
virtual ~AliAnaPi0()
Destructor.
Definition: AliAnaPi0.cxx:184
TH2F * fhPrimEtaProdVertex
! Spectrum of primary eta vs production vertex
Definition: AliAnaPi0.h:483
Bool_t fSameSM
Select only pairs in same SM;.
Definition: AliAnaPi0.h:207
TH2F * fhMiOpAngleBinMinClusterTimePerSM[10]
! time of lowest energy cluster in pair, depending on opening angle cut, y axis is SM number...
Definition: AliAnaPi0.h:527
Bool_t fCheckAccInSector
Check that the decay pi0 falls in the same SM or sector.
Definition: AliAnaPi0.h:220
TH2F * fhReSecondaryCellOutTimeWindow
! Combine clusters when at least one significant cells in cluster has t > 50 ns, same event ...
Definition: AliAnaPi0.h:499
TH1I * fhEventBin
! Number of real pairs in a particular bin (cen,vz,rp)
Definition: AliAnaPi0.h:373
TH1F * fhPrimPi0Pt
! Spectrum of Primary
Definition: AliAnaPi0.h:394
TH2F ** fhRe2
[GetNCentrBin()*fNPIDBits*fNAsymCuts]
Definition: AliAnaPi0.h:295
TH1F * fhPrimEtaAccE
! Spectrum of primary with accepted daughters
Definition: AliAnaPi0.h:419
TH2F ** fhMCEtaMassPtRec
[fNPtCuts*fNAsymCuts*fNCellNCuts]
Definition: AliAnaPi0.h:467
Bool_t fUseAngleCut
Select pairs depending on their opening angle.
Definition: AliAnaPi0.h:183
ClassImp(AliAnalysisTaskCRC) AliAnalysisTaskCRC
TH2F * fhMiOpAngleBinPairClusterMassPerSM[10]
! cluster pair mass, depending on opening angle cut, y axis is SM number, mixed event ...
Definition: AliAnaPi0.h:533
TH2F * fhRePtAsymEta
! REAL two-photon pt vs asymmetry, close to eta mass
Definition: AliAnaPi0.h:366
virtual Int_t GetNRPBin() const
Number of bins in vertex.
TH2F * fhEventPlaneResolution
! Histogram with Event plane resolution vs centrality
Definition: AliAnaPi0.h:378
TH2F * fhMixedOpeningAngle
! Opening angle of pair versus pair energy
Definition: AliAnaPi0.h:384
const Double_t etamax
TH2F ** fhMCPi0PtTruePtRec
[fNPtCuts*fNAsymCuts*fNCellNCuts]
Definition: AliAnaPi0.h:464
TH2F * fhPrimEtaYeta
! PseudoRapidity distribution of primary particles vs pT
Definition: AliAnaPi0.h:425
virtual AliMCAnalysisUtils * GetMCAnalysisUtils()
Float_t GetHistoAsymmetryMin() const
TH2F * fhCosThStarPrimPi0
! cos(theta*) plots vs E for primary pi0, same as asymmetry ...
Definition: AliAnaPi0.h:491
virtual void Print(const Option_t *) const
Print some relevant parameters set for the analysis.
Int_t GetHistoTimeBins() const
TH2F ** fhReDiffSectorDCALPHOSMod
[6]
Definition: AliAnaPi0.h:260
TH2F * fhMCPi0ProdVertex
! Spectrum of selected pi0 vs production vertex
Definition: AliAnaPi0.h:480
const char Option_t
Definition: External.C:48
TH2F * fhPrimPi0PtEventPlane
! primary pi0 reconstructed event plane vs pT
Definition: AliAnaPi0.h:411
TH2F * fhMiSecondaryCellInTimeWindow
! Combine clusters when all significant cells in cluster have t < 50 ns, different events ...
Definition: AliAnaPi0.h:498
Float_t GetHistoTimeMax() const
Float_t GetHistoTimeMin() const
TList ** fEventsList
Containers for photons in stored events.
Definition: AliAnaPi0.h:179
virtual AliVCluster * FindCluster(TObjArray *clusters, Int_t clId, Int_t &iclus, Int_t first=0)
Float_t GetHistoOpeningAngleMin() const
TH2F * fhMCEtaProdVertex
! Spectrum of selected eta vs production vertex
Definition: AliAnaPi0.h:481
TH2F * fhPrimEtaOpeningAngleAsym
! Opening angle of pair versus pair E asymmetry, eta primaries
Definition: AliAnaPi0.h:432
TH2F ** fhMiSameSectorDCALPHOSMod
[fNModules-1]
Definition: AliAnaPi0.h:275
Float_t GetHistoPhiMax() const
TH2F * fhReSecondaryCellInTimeWindow
! Combine clusters when all significant cells in cluster have t < 50 ns, same event ...
Definition: AliAnaPi0.h:497
TH2F * fhReOpAngleBinMinClusterColRow[10]
! Column and row location of main cell of lowest energy cluster in pair, depending on opening angle c...
Definition: AliAnaPi0.h:505
bool Bool_t
Definition: External.C:53
virtual AliCaloTrackReader * GetReader() const
Bool_t fFillSecondaryCellTiming
Fill histograms depending on timing of secondary cells in clusters.
Definition: AliAnaPi0.h:217
TH2F ** fhMiInvPt3
[GetNCentrBin()*fNPIDBits*fNAsymCuts]
Definition: AliAnaPi0.h:331
AliAnaPi0()
Default Constructor. Initialized parameters with default values.
Definition: AliAnaPi0.cxx:53
Int_t GetHistoEtaBins() const
TH2F ** fhMi3
[GetNCentrBin()*fNPIDBits*fNAsymCuts]
Definition: AliAnaPi0.h:307
virtual TObjArray * GetEMCALClusters() const
Bool_t fMultiCutAnaAcc
Do analysis with several or fixed cut, acceptance plots (eta-phi, col-row)
Definition: AliAnaPi0.h:191
TH2F * fhPtBinClusterColRow[10]
! Column and row location of cluster in different energy bins.
Definition: AliAnaPi0.h:537
TH2F * fhArmPrimEta[4]
! Armenteros plots for primary eta in 6 energy bins
Definition: AliAnaPi0.h:490
Bool_t fFillSSCombinations
Do invariant mass for different combination of shower shape clusters.
Definition: AliAnaPi0.h:211
TH2F * fhPrimEtaPtOrigin
! Spectrum of generated eta vs mother
Definition: AliAnaPi0.h:442
TH2F * fhMiOpAngleBinMaxClusterEtaPhi[10]
! Eta-Phi location of highest energy cluster in pair, depending on opening angle cut, mixed event
Definition: AliAnaPi0.h:522
TH2F ** fhMiInvPt2
[GetNCentrBin()*fNPIDBits*fNAsymCuts]
Definition: AliAnaPi0.h:323
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
Int_t GetMaxEnergyCell(AliVCaloCells *cells, AliVCluster *clu, Float_t &fraction) const
For a given CaloCluster, it gets the absId of the cell with maximum energy deposit.
TH2F * fhMCOrgDeltaPhi[13]
! Delta Phi vs pt of real pairs, check common origin of pair
Definition: AliAnaPi0.h:453
Float_t fPtCutsMax[11]
Array with different pt cuts, maximum.
Definition: AliAnaPi0.h:194
Int_t nptbins
TH2F ** fhMiMod
[8]
Definition: AliAnaPi0.h:263
TH2F * fhMiOpAngleBinMinClusterNCellPerSM[10]
! N cells of lowest energy cluster in pair, depending on opening angle cut, y axis is SM number...
Definition: AliAnaPi0.h:529
virtual AliMixedEvent * GetMixedEvent() const
TH2F * fhReConv2
! REAL two-photon invariant mass distribution both pair photons recombined from 2 clusters with small...
Definition: AliAnaPi0.h:284
const Double_t phimin
TH2F * fhMiPtAsymEta
! Mix two-photon pt vs asymmetry, close to eta mass
Definition: AliAnaPi0.h:369
Int_t fNAngleCutBins
Number of angle cuts bins.
Definition: AliAnaPi0.h:202
TH2F * fhPrimEtaPhi
! Azimutal distribution of primary particles vs pT
Definition: AliAnaPi0.h:428
TH2F * fhPrimPi0OpeningAnglePhotonCuts
! Opening angle of pair versus pair energy, primaries, photon pt cuts
Definition: AliAnaPi0.h:407
TH2F ** fhMCPi0MassPtTrue
[fNPtCuts*fNAsymCuts*fNCellNCuts]
Definition: AliAnaPi0.h:461