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