AliPhysics  323987d (323987d)
 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  fhPairGeneratorsBkgMassMCPi0[igen][itag] = new TH2F
2101  (Form("h%sGeneratorPairMass%s%s_MCPi0",
2102  tagBkgNames[itag].Data(),add.Data(),GetCocktailGenNameToCheck(igen).Data()),
2103  Form("Pair Mass with contribution of true #pi^{0} generator%s, %s ",
2104  GetCocktailGenNameToCheck(igen).Data(),tagBkgTitle[itag].Data()),
2105  nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
2106  fhPairGeneratorsBkgMassMCPi0[igen][itag]->SetYTitle("#it{M} (MeV/#it{c}^{2})");
2107  fhPairGeneratorsBkgMassMCPi0[igen][itag]->SetXTitle("#it{E}_{reco} (GeV)");
2108  outputContainer->Add(fhPairGeneratorsBkgMassMCPi0[igen][itag]) ;
2109 
2110  fhPairGeneratorsBkgMassMCEta[igen][itag] = new TH2F
2111  (Form("h%sGeneratorPairMass%s%s_MCEta",
2112  tagBkgNames[itag].Data(),add.Data(),GetCocktailGenNameToCheck(igen).Data()),
2113  Form("Pair Mass with contribution of true #eta generator%s, %s ",
2114  GetCocktailGenNameToCheck(igen).Data(),tagBkgTitle[itag].Data()),
2115  nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
2116  fhPairGeneratorsBkgMassMCEta[igen][itag]->SetYTitle("#it{M} (MeV/#it{c}^{2})");
2117  fhPairGeneratorsBkgMassMCEta[igen][itag]->SetXTitle("#it{E}_{reco} (GeV)");
2118  outputContainer->Add(fhPairGeneratorsBkgMassMCEta[igen][itag]) ;
2119 
2121  (Form("h%sGeneratorPairEPrimRecoRatio%s%s_MCPi0",
2122  tagBkgNames[itag].Data(),add.Data(),GetCocktailGenNameToCheck(igen).Data()),
2123  Form("#it{E}_{reco}/#it{E}_{gen} pair with contribution of true #pi^{0} generator%s, %s ",
2124  GetCocktailGenNameToCheck(igen).Data(),tagBkgTitle[itag].Data()),
2125  nptbins,ptmin,ptmax,nratbins,ratmin,ratmax);
2126  fhPairGeneratorsBkgEPrimRecoRatioMCPi0[igen][itag]->SetYTitle("#it{E}_{reco}/#it{E}_{gen}");
2127  fhPairGeneratorsBkgEPrimRecoRatioMCPi0[igen][itag]->SetXTitle("#it{E}_{reco} (GeV)");
2128  outputContainer->Add(fhPairGeneratorsBkgEPrimRecoRatioMCPi0[igen][itag]) ;
2129 
2131  (Form("h%sGeneratorPairEPrimRecoRatio%s%s_MCEta",
2132  tagBkgNames[itag].Data(),add.Data(),GetCocktailGenNameToCheck(igen).Data()),
2133  Form("#it{E}_{reco}/#it{E}_{gen} pair with contribution of true #eta generator%s, %s ",
2134  GetCocktailGenNameToCheck(igen).Data(),tagBkgTitle[itag].Data()),
2135  nptbins,ptmin,ptmax,nratbins,ratmin,ratmax);
2136  fhPairGeneratorsBkgEPrimRecoRatioMCEta[igen][itag]->SetYTitle("#it{E}_{reco}/#it{E}_{gen}");
2137  fhPairGeneratorsBkgEPrimRecoRatioMCEta[igen][itag]->SetXTitle("#it{E}_{reco} (GeV)");
2138  outputContainer->Add(fhPairGeneratorsBkgEPrimRecoRatioMCEta[igen][itag]) ;
2139 
2141  (Form("h%sGeneratorPairEPrimRecoDiff%s%s_MCPi0",
2142  tagBkgNames[itag].Data(),add.Data(),GetCocktailGenNameToCheck(igen).Data()),
2143  Form("#it{E}_{reco}-#it{E}_{gen} pair with contribution of true #pi^{0} generator%s, %s ",
2144  GetCocktailGenNameToCheck(igen).Data(),tagBkgTitle[itag].Data()),
2145  nptbins,ptmin,ptmax,ndifbins,difmin,difmax);
2146  fhPairGeneratorsBkgEPrimRecoDiffMCPi0[igen][itag]->SetYTitle("#it{E}_{reco}-#it{E}_{gen} (GeV)");
2147  fhPairGeneratorsBkgEPrimRecoDiffMCPi0[igen][itag]->SetXTitle("#it{E}_{reco} (GeV)");
2148  outputContainer->Add(fhPairGeneratorsBkgEPrimRecoDiffMCPi0[igen][itag]) ;
2149 
2151  (Form("h%sGeneratorPairEPrimRecoDiff%s%s_MCEta",
2152  tagBkgNames[itag].Data(),add.Data(),GetCocktailGenNameToCheck(igen).Data()),
2153  Form("#it{E}_{reco}-#it{E}_{gen} pair with contribution of true #eta generator%s, %s ",
2154  GetCocktailGenNameToCheck(igen).Data(),tagBkgTitle[itag].Data()),
2155  nptbins,ptmin,ptmax,ndifbins,difmin,difmax);
2156  fhPairGeneratorsBkgEPrimRecoDiffMCEta[igen][itag]->SetYTitle("#it{E}_{reco}-#it{E}_{gen} (GeV)");
2157  fhPairGeneratorsBkgEPrimRecoDiffMCEta[igen][itag]->SetXTitle("#it{E}_{reco} (GeV)");
2158  outputContainer->Add(fhPairGeneratorsBkgEPrimRecoDiffMCEta[igen][itag]) ;
2159 
2161  (Form("h%sGeneratorPairEPrimRecoRatio%s%s_MCPi0MassCut",
2162  tagBkgNames[itag].Data(),add.Data(),GetCocktailGenNameToCheck(igen).Data()),
2163  Form("#it{E}_{reco}/#it{E}_{gen} pair with contribution of true #pi^{0} generator%s, %s ",
2164  GetCocktailGenNameToCheck(igen).Data(),tagBkgTitle[itag].Data()),
2165  nptbins,ptmin,ptmax,nratbins,ratmin,ratmax);
2166  fhPairGeneratorsBkgEPrimRecoRatioMCPi0MassCut[igen][itag]->SetYTitle("#it{E}_{reco}/#it{E}_{gen}");
2167  fhPairGeneratorsBkgEPrimRecoRatioMCPi0MassCut[igen][itag]->SetXTitle("#it{E}_{reco} (GeV)");
2168  outputContainer->Add(fhPairGeneratorsBkgEPrimRecoRatioMCPi0MassCut[igen][itag]) ;
2169 
2171  (Form("h%sGeneratorPairEPrimRecoRatio%s%s_MCEtaMassCut",
2172  tagBkgNames[itag].Data(),add.Data(),GetCocktailGenNameToCheck(igen).Data()),
2173  Form("#it{E}_{reco}/#it{E}_{gen} pair with contribution of true #eta generator%s, %s ",
2174  GetCocktailGenNameToCheck(igen).Data(),tagBkgTitle[itag].Data()),
2175  nptbins,ptmin,ptmax,nratbins,ratmin,ratmax);
2176  fhPairGeneratorsBkgEPrimRecoRatioMCEtaMassCut[igen][itag]->SetYTitle("#it{E}_{reco}/#it{E}_{gen}");
2177  fhPairGeneratorsBkgEPrimRecoRatioMCEtaMassCut[igen][itag]->SetXTitle("#it{E}_{reco} (GeV)");
2178  outputContainer->Add(fhPairGeneratorsBkgEPrimRecoRatioMCEtaMassCut[igen][itag]) ;
2179 
2181  (Form("h%sGeneratorPairEPrimRecoDiff%s%s_MCPi0MassCut",
2182  tagBkgNames[itag].Data(),add.Data(),GetCocktailGenNameToCheck(igen).Data()),
2183  Form("#it{E}_{reco}-#it{E}_{gen} pair with contribution of true #pi^{0} generator%s, %s ",
2184  GetCocktailGenNameToCheck(igen).Data(),tagBkgTitle[itag].Data()),
2185  nptbins,ptmin,ptmax,ndifbins,difmin,difmax);
2186  fhPairGeneratorsBkgEPrimRecoDiffMCPi0MassCut[igen][itag]->SetYTitle("#it{E}_{reco}-#it{E}_{gen} (GeV)");
2187  fhPairGeneratorsBkgEPrimRecoDiffMCPi0MassCut[igen][itag]->SetXTitle("#it{E}_{reco} (GeV)");
2188  outputContainer->Add(fhPairGeneratorsBkgEPrimRecoDiffMCPi0MassCut[igen][itag]) ;
2189 
2191  (Form("h%sGeneratorPairEPrimRecoDiff%s%s_MCEtaMassCut",
2192  tagBkgNames[itag].Data(),add.Data(),GetCocktailGenNameToCheck(igen).Data()),
2193  Form("#it{E}_{reco}-#it{E}_{gen} pair with contribution of true #eta generator%s, %s ",
2194  GetCocktailGenNameToCheck(igen).Data(),tagBkgTitle[itag].Data()),
2195  nptbins,ptmin,ptmax,ndifbins,difmin,difmax);
2196  fhPairGeneratorsBkgEPrimRecoDiffMCEtaMassCut[igen][itag]->SetYTitle("#it{E}_{reco}-#it{E}_{gen} (GeV)");
2197  fhPairGeneratorsBkgEPrimRecoDiffMCEtaMassCut[igen][itag]->SetXTitle("#it{E}_{reco} (GeV)");
2198  outputContainer->Add(fhPairGeneratorsBkgEPrimRecoDiffMCEtaMassCut[igen][itag]) ;
2199  }
2200  }
2201  }
2202 
2203 // for(Int_t i = 0; i < outputContainer->GetEntries() ; i++){
2204 //
2205 // printf("Histogram %d, name: %s\n ",i, outputContainer->At(i)->GetName());
2206 //
2207 // }
2208 
2209  return outputContainer;
2210 }
2211 
2212 //___________________________________________________
2214 //___________________________________________________
2215 void AliAnaPi0::Print(const Option_t * /*opt*/) const
2216 {
2217  printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
2219 
2220  printf("Number of bins in Centrality: %d \n",GetNCentrBin()) ;
2221  printf("Number of bins in Z vert. pos: %d \n",GetNZvertBin()) ;
2222  printf("Number of bins in Reac. Plain: %d \n",GetNRPBin()) ;
2223  printf("Depth of event buffer: %d \n",GetNMaxEvMix()) ;
2224  printf("Pair in same Module: %d \n",fSameSM) ;
2225  printf("Cuts: \n") ;
2226  // printf("Z vertex position: -%2.3f < z < %2.3f \n",GetZvertexCut(),GetZvertexCut()) ; //It crashes here, why?
2227  printf("Number of modules: %d \n",fNModules) ;
2228  printf("Select pairs with their angle: %d, edep %d, min angle %2.3f, max angle %2.3f \n",fUseAngleCut, fUseAngleEDepCut, fAngleCut, fAngleMaxCut) ;
2229  printf("Asymmetry cuts: n = %d, \n",fNAsymCuts) ;
2230  printf("\tasymmetry < ");
2231  for(Int_t i = 0; i < fNAsymCuts; i++) printf("%2.2f ",fAsymCuts[i]);
2232  printf("\n");
2233 
2234  printf("PID selection bits: n = %d, \n",fNPIDBits) ;
2235  printf("\tPID bit = ");
2236  for(Int_t i = 0; i < fNPIDBits; i++) printf("%d ",fPIDBits[i]);
2237  printf("\n");
2238 
2240  {
2241  printf("pT cuts: n = %d, \n",fNPtCuts) ;
2242  printf("\tpT > ");
2243  for(Int_t i = 0; i < fNPtCuts; i++) printf("%2.2f ",fPtCuts[i]);
2244  printf("GeV/c\n");
2245  printf("\tpT < ");
2246  for(Int_t i = 0; i < fNPtCuts; i++) printf("%2.2f ",fPtCutsMax[i]);
2247  printf("GeV/c\n");
2248 
2249  printf("N cell in cluster cuts: n = %d, \n",fNCellNCuts) ;
2250  printf("\tnCell > ");
2251  for(Int_t i = 0; i < fNCellNCuts; i++) printf("%d ",fCellNCuts[i]);
2252  printf("\n");
2253 
2254  }
2255  printf("------------------------------------------------------\n") ;
2256 }
2257 
2258 //________________________________________
2260 //________________________________________
2262 {
2263  Double_t mesonY = -100 ;
2264  Double_t mesonE = -1 ;
2265  Double_t mesonPt = -1 ;
2266  Double_t mesonPhi = 100 ;
2267  Double_t mesonYeta= -1 ;
2268 
2269  Int_t pdg = 0 ;
2270  Int_t nprim = 0 ;
2271  Int_t nDaught = 0 ;
2272  Int_t iphot1 = 0 ;
2273  Int_t iphot2 = 0 ;
2274 
2275  Float_t cen = GetEventCentrality();
2276  Float_t ep = GetEventPlaneAngle();
2277 
2278  TParticle * primStack = 0;
2279  AliAODMCParticle * primAOD = 0;
2280 
2281  // Get the ESD MC particles container
2282  AliStack * stack = 0;
2283  if( GetReader()->ReadStack() )
2284  {
2285  stack = GetMCStack();
2286  if(!stack ) return;
2287  nprim = stack->GetNtrack();
2288  }
2289 
2290  // Get the AOD MC particles container
2291  TClonesArray * mcparticles = 0;
2292  if( GetReader()->ReadAODMCParticles() )
2293  {
2294  mcparticles = GetReader()->GetAODMCParticles();
2295  if( !mcparticles ) return;
2296  nprim = mcparticles->GetEntriesFast();
2297  }
2298 
2299  for(Int_t i=0 ; i < nprim; i++)
2300  {
2301  if ( !GetReader()->AcceptParticleMCLabel( i ) ) continue ;
2302 
2303  if(GetReader()->ReadStack())
2304  {
2305  primStack = stack->Particle(i) ;
2306  if(!primStack)
2307  {
2308  AliWarning("ESD primaries pointer not available!!");
2309  continue;
2310  }
2311 
2312  // If too small skip
2313  if( primStack->Energy() < 0.4 ) continue;
2314 
2315  pdg = primStack->GetPdgCode();
2316  nDaught = primStack->GetNDaughters();
2317  iphot1 = primStack->GetDaughter(0) ;
2318  iphot2 = primStack->GetDaughter(1) ;
2319 
2320  // Protection against floating point exception
2321  if ( primStack->Energy() == TMath::Abs(primStack->Pz()) ||
2322  (primStack->Energy() - primStack->Pz()) < 1e-3 ||
2323  (primStack->Energy() + primStack->Pz()) < 0 ) continue ;
2324 
2325  //printf("i %d, %s %d %s %d \n",i, stack->Particle(i)->GetName(), stack->Particle(i)->GetPdgCode(),
2326  // prim->GetName(), prim->GetPdgCode());
2327 
2328  //Photon kinematics
2329  primStack->Momentum(fMCPrimMesonMom);
2330 
2331  mesonY = 0.5*TMath::Log((primStack->Energy()+primStack->Pz())/(primStack->Energy()-primStack->Pz())) ;
2332  }
2333  else
2334  {
2335  primAOD = (AliAODMCParticle *) mcparticles->At(i);
2336  if(!primAOD)
2337  {
2338  AliWarning("AOD primaries pointer not available!!");
2339  continue;
2340  }
2341 
2342  // If too small skip
2343  if( primAOD->E() < 0.4 ) continue;
2344 
2345  pdg = primAOD->GetPdgCode();
2346  nDaught = primAOD->GetNDaughters();
2347  iphot1 = primAOD->GetFirstDaughter() ;
2348  iphot2 = primAOD->GetLastDaughter() ;
2349 
2350  // Protection against floating point exception
2351  if ( primAOD->E() == TMath::Abs(primAOD->Pz()) ||
2352  (primAOD->E() - primAOD->Pz()) < 1e-3 ||
2353  (primAOD->E() + primAOD->Pz()) < 0 ) continue ;
2354 
2355  // Photon kinematics
2356  fMCPrimMesonMom.SetPxPyPzE(primAOD->Px(),primAOD->Py(),primAOD->Pz(),primAOD->E());
2357 
2358  mesonY = 0.5*TMath::Log((primAOD->E()+primAOD->Pz())/(primAOD->E()-primAOD->Pz())) ;
2359  }
2360 
2361  // Select only pi0 or eta
2362  if( pdg != 111 && pdg != 221) continue ;
2363 
2364  mesonPt = fMCPrimMesonMom.Pt () ;
2365  mesonE = fMCPrimMesonMom.E () ;
2366  mesonYeta= fMCPrimMesonMom.Eta() ;
2367  mesonPhi = fMCPrimMesonMom.Phi() ;
2368  if( mesonPhi < 0 ) mesonPhi+=TMath::TwoPi();
2369  mesonPhi *= TMath::RadToDeg();
2370 
2371  if(pdg == 111)
2372  {
2373  if(TMath::Abs(mesonY) < 1.0)
2374  {
2375  fhPrimPi0E ->Fill(mesonE , GetEventWeight()) ;
2376  fhPrimPi0Pt ->Fill(mesonPt, GetEventWeight()) ;
2377  fhPrimPi0Phi->Fill(mesonPt, mesonPhi, GetEventWeight()) ;
2378 
2379  fhPrimPi0YetaYcut->Fill(mesonPt, mesonYeta, GetEventWeight()) ;
2380 
2382  {
2383  fhPrimPi0PtCentrality->Fill(mesonPt, cen, GetEventWeight()) ;
2384  fhPrimPi0PtEventPlane->Fill(mesonPt, ep , GetEventWeight()) ;
2385  }
2386  }
2387 
2388  fhPrimPi0Y ->Fill(mesonPt, mesonY , GetEventWeight()) ;
2389  fhPrimPi0Yeta->Fill(mesonPt, mesonYeta, GetEventWeight()) ;
2390  }
2391  else if(pdg == 221)
2392  {
2393  if(TMath::Abs(mesonY) < 1.0)
2394  {
2395  fhPrimEtaE ->Fill(mesonE , GetEventWeight()) ;
2396  fhPrimEtaPt ->Fill(mesonPt, GetEventWeight()) ;
2397  fhPrimEtaPhi->Fill(mesonPt, mesonPhi, GetEventWeight()) ;
2398 
2399  fhPrimEtaYetaYcut->Fill(mesonPt, mesonYeta, GetEventWeight()) ;
2400 
2402  {
2403  fhPrimEtaPtCentrality->Fill(mesonPt, cen, GetEventWeight()) ;
2404  fhPrimEtaPtEventPlane->Fill(mesonPt, ep , GetEventWeight()) ;
2405  }
2406  }
2407 
2408  fhPrimEtaY ->Fill(mesonPt, mesonY , GetEventWeight()) ;
2409  fhPrimEtaYeta->Fill(mesonPt, mesonYeta, GetEventWeight()) ;
2410  }
2411 
2412  // Origin of meson
2413  if(fFillOriginHisto && TMath::Abs(mesonY) < 0.7)
2414  {
2415  Int_t momindex = -1;
2416  Int_t mompdg = -1;
2417  Int_t momstatus = -1;
2418  Int_t status = -1;
2419  Float_t momR = 0;
2420  if(GetReader()->ReadStack()) momindex = primStack->GetFirstMother();
2421  if(GetReader()->ReadAODMCParticles()) momindex = primAOD ->GetMother();
2422 
2423  if(momindex >= 0 && momindex < nprim)
2424  {
2425  if(GetReader()->ReadStack())
2426  {
2427  status = primStack->GetStatusCode();
2428  TParticle* mother = stack->Particle(momindex);
2429  mompdg = TMath::Abs(mother->GetPdgCode());
2430  momstatus = mother->GetStatusCode();
2431  momR = mother->R();
2432  }
2433 
2434  if(GetReader()->ReadAODMCParticles())
2435  {
2436  status = primAOD->GetStatus();
2437  AliAODMCParticle* mother = (AliAODMCParticle*) mcparticles->At(momindex);
2438  mompdg = TMath::Abs(mother->GetPdgCode());
2439  momstatus = mother->GetStatus();
2440  momR = TMath::Sqrt(mother->Xv()*mother->Xv()+mother->Yv()*mother->Yv());
2441  }
2442 
2443  if(pdg == 111)
2444  {
2445  fhPrimPi0ProdVertex->Fill(mesonPt, momR , GetEventWeight());
2446  fhPrimPi0PtStatus ->Fill(mesonPt, status, GetEventWeight());
2447 
2448 
2449  if (momstatus == 21) fhPrimPi0PtOrigin->Fill(mesonPt, 0.5, GetEventWeight());//parton
2450  else if(mompdg < 22 ) fhPrimPi0PtOrigin->Fill(mesonPt, 1.5, GetEventWeight());//quark
2451  else if(mompdg > 2100 && mompdg < 2210)
2452  fhPrimPi0PtOrigin->Fill(mesonPt, 2.5, GetEventWeight());// resonances
2453  else if(mompdg == 221) fhPrimPi0PtOrigin->Fill(mesonPt, 8.5, GetEventWeight());//eta
2454  else if(mompdg == 331) fhPrimPi0PtOrigin->Fill(mesonPt, 9.5, GetEventWeight());//eta prime
2455  else if(mompdg == 213) fhPrimPi0PtOrigin->Fill(mesonPt, 4.5, GetEventWeight());//rho
2456  else if(mompdg == 223) fhPrimPi0PtOrigin->Fill(mesonPt, 5.5, GetEventWeight());//omega
2457  else if(mompdg >= 310 && mompdg <= 323)
2458  fhPrimPi0PtOrigin->Fill(mesonPt, 6.5, GetEventWeight());//k0S, k+-,k*
2459  else if(mompdg == 130) fhPrimPi0PtOrigin->Fill(mesonPt, 6.5, GetEventWeight());//k0L
2460  else if(momstatus == 11 || momstatus == 12 )
2461  fhPrimPi0PtOrigin->Fill(mesonPt, 3.5, GetEventWeight());//resonances
2462  else fhPrimPi0PtOrigin->Fill(mesonPt, 7.5, GetEventWeight());//other?
2463 
2464  //printf("Prim Pi0: index %d, pt %2.2f Prod vertex %3.3f, origin pdg %d, origin status %d, origin UI %d\n",
2465  // momindex, mesonPt,mother->R(),mompdg,momstatus,mother->GetUniqueID());
2466 
2467  if(status!=11)
2468  {
2469  if (momstatus == 21) fhPrimNotResonancePi0PtOrigin->Fill(mesonPt, 0.5, GetEventWeight());//parton
2470  else if(mompdg < 22 ) fhPrimNotResonancePi0PtOrigin->Fill(mesonPt, 1.5, GetEventWeight());//quark
2471  else if(mompdg > 2100 && mompdg < 2210)
2472  fhPrimNotResonancePi0PtOrigin->Fill(mesonPt, 2.5, GetEventWeight());// resonances
2473  else if(mompdg == 221) fhPrimNotResonancePi0PtOrigin->Fill(mesonPt, 8.5, GetEventWeight());//eta
2474  else if(mompdg == 331) fhPrimNotResonancePi0PtOrigin->Fill(mesonPt, 9.5, GetEventWeight());//eta prime
2475  else if(mompdg == 213) fhPrimNotResonancePi0PtOrigin->Fill(mesonPt, 4.5, GetEventWeight());//rho
2476  else if(mompdg == 223) fhPrimNotResonancePi0PtOrigin->Fill(mesonPt, 5.5, GetEventWeight());//omega
2477  else if(mompdg >= 310 && mompdg <= 323)
2478  fhPrimNotResonancePi0PtOrigin->Fill(mesonPt, 6.5, GetEventWeight());//k0S, k+-,k*
2479  else if(mompdg == 130) fhPrimNotResonancePi0PtOrigin->Fill(mesonPt, 6.5, GetEventWeight());//k0L
2480  else if(momstatus == 11 || momstatus == 12 )
2481  fhPrimNotResonancePi0PtOrigin->Fill(mesonPt, 3.5, GetEventWeight());//resonances
2482  else fhPrimNotResonancePi0PtOrigin->Fill(mesonPt, 7.5, GetEventWeight());//other?
2483  }
2484 
2485  }//pi0
2486  else
2487  {
2488  if (momstatus == 21 ) fhPrimEtaPtOrigin->Fill(mesonPt, 0.5, GetEventWeight());//parton
2489  else if(mompdg < 22 ) fhPrimEtaPtOrigin->Fill(mesonPt, 1.5, GetEventWeight());//quark
2490  else if(mompdg > 2100 && mompdg < 2210)
2491  fhPrimEtaPtOrigin->Fill(mesonPt, 2.5, GetEventWeight());//qq resonances
2492  else if(mompdg == 331) fhPrimEtaPtOrigin->Fill(mesonPt, 5.5, GetEventWeight());//eta prime
2493  else if(momstatus == 11 || momstatus == 12 )
2494  fhPrimEtaPtOrigin->Fill(mesonPt, 3.5, GetEventWeight());//resonances
2495  else fhPrimEtaPtOrigin->Fill(mesonPt, 4.5, GetEventWeight());//stable, conversions?
2496  //printf("Other Meson pdg %d, Mother %s, pdg %d, status %d\n",pdg, TDatabasePDG::Instance()->GetParticle(mompdg)->GetName(),mompdg, momstatus );
2497 
2498  fhPrimEtaProdVertex->Fill(mesonPt, momR, GetEventWeight());
2499  }
2500  } // pi0 has mother
2501  }
2502 
2503  // Check if both photons hit calorimeter or a given fiducial region
2504  // only if those settings are specified.
2505  if ( !IsFiducialCutOn() && !IsRealCaloAcceptanceOn() ) continue ;
2506 
2507  if ( nDaught != 2 ) continue; //Only interested in 2 gamma decay
2508 
2509  if ( iphot1 < 0 || iphot1 >= nprim || iphot2 < 0 || iphot2 >= nprim ) continue ;
2510 
2511  Int_t pdg1 = 0;
2512  Int_t pdg2 = 0;
2513  Bool_t inacceptance1 = kTRUE;
2514  Bool_t inacceptance2 = kTRUE;
2515 
2516  if(GetReader()->ReadStack())
2517  {
2518  TParticle * phot1 = stack->Particle(iphot1) ;
2519  TParticle * phot2 = stack->Particle(iphot2) ;
2520 
2521  if(!phot1 || !phot2) continue ;
2522 
2523  pdg1 = phot1->GetPdgCode();
2524  pdg2 = phot2->GetPdgCode();
2525 
2526  phot1->Momentum(fPhotonMom1);
2527  phot2->Momentum(fPhotonMom2);
2528 
2529  // Check if photons hit the Calorimeter acceptance
2531  {
2532  if( !GetCaloUtils()->IsMCParticleInCalorimeterAcceptance( GetCalorimeter(), phot1 )) inacceptance1 = kFALSE ;
2533  if( !GetCaloUtils()->IsMCParticleInCalorimeterAcceptance( GetCalorimeter(), phot2 )) inacceptance2 = kFALSE ;
2534  }
2535  }
2536 
2537  if(GetReader()->ReadAODMCParticles())
2538  {
2539  AliAODMCParticle * phot1 = (AliAODMCParticle *) mcparticles->At(iphot1) ;
2540  AliAODMCParticle * phot2 = (AliAODMCParticle *) mcparticles->At(iphot2) ;
2541 
2542  if(!phot1 || !phot2) continue ;
2543 
2544  pdg1 = phot1->GetPdgCode();
2545  pdg2 = phot2->GetPdgCode();
2546 
2547  fPhotonMom1.SetPxPyPzE(phot1->Px(),phot1->Py(),phot1->Pz(),phot1->E());
2548  fPhotonMom2.SetPxPyPzE(phot2->Px(),phot2->Py(),phot2->Pz(),phot2->E());
2549 
2550  // Check if photons hit the Calorimeter acceptance
2552  {
2553  if( !GetCaloUtils()->IsMCParticleInCalorimeterAcceptance( GetCalorimeter(), phot1 )) inacceptance1 = kFALSE ;
2554  if( !GetCaloUtils()->IsMCParticleInCalorimeterAcceptance( GetCalorimeter(), phot2 )) inacceptance2 = kFALSE ;
2555  }
2556  }
2557 
2558  if( pdg1 != 22 || pdg2 !=22) continue ;
2559 
2560  // Check if photons hit desired acceptance in the fidutial borders fixed in the analysis
2561  if(IsFiducialCutOn())
2562  {
2563  if( inacceptance1 && !GetFiducialCut()->IsInFiducialCut(fPhotonMom1.Eta(), fPhotonMom1.Phi(), GetCalorimeter()) ) inacceptance1 = kFALSE ;
2564  if( inacceptance2 && !GetFiducialCut()->IsInFiducialCut(fPhotonMom2.Eta(), fPhotonMom2.Phi(), GetCalorimeter()) ) inacceptance2 = kFALSE ;
2565  }
2566 
2568 
2569  if(GetCalorimeter()==kEMCAL && inacceptance1 && inacceptance2 && fCheckAccInSector)
2570  {
2571  Int_t absID1=0;
2572  Int_t absID2=0;
2573 
2574  Float_t photonPhi1 = fPhotonMom1.Phi();
2575  Float_t photonPhi2 = fPhotonMom2.Phi();
2576 
2577  if(photonPhi1 < 0) photonPhi1+=TMath::TwoPi();
2578  if(photonPhi2 < 0) photonPhi2+=TMath::TwoPi();
2579 
2580  GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(fPhotonMom1.Eta(),photonPhi1,absID1);
2581  GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(fPhotonMom2.Eta(),photonPhi2,absID2);
2582 
2583  Int_t sm1 = GetEMCALGeometry()->GetSuperModuleNumber(absID1);
2584  Int_t sm2 = GetEMCALGeometry()->GetSuperModuleNumber(absID2);
2585 
2586  Int_t j=0;
2587  Bool_t sameSector = kFALSE;
2588  for(Int_t isector = 0; isector < fNModules/2; isector++)
2589  {
2590  j=2*isector;
2591  if((sm1==j && sm2==j+1) || (sm1==j+1 && sm2==j)) sameSector = kTRUE;
2592  }
2593 
2594  if(sm1!=sm2 && !sameSector)
2595  {
2596  inacceptance1 = kFALSE;
2597  inacceptance2 = kFALSE;
2598  }
2599  //if(sm1!=sm2)printf("sm1 %d, sm2 %d, same sector %d, in acceptance %d\n",sm1,sm2,sameSector,inacceptance);
2600  // if(GetEMCALGeometry()->Impact(phot1) && GetEMCALGeometry()->Impact(phot2))
2601  // inacceptance = kTRUE;
2602  }
2603 
2604  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",
2606  mesonPt,mesonYeta,mesonPhi,
2607  fPhotonMom1.Pt(),fPhotonMom1.Eta(),fPhotonMom1.Phi()*TMath::RadToDeg(),
2608  fPhotonMom2.Pt(),fPhotonMom2.Eta(),fPhotonMom2.Phi()*TMath::RadToDeg(),
2609  inacceptance1, inacceptance2));
2610 
2611  if(inacceptance1 && inacceptance2)
2612  {
2613  Float_t asym = TMath::Abs((fPhotonMom1.E()-fPhotonMom2.E()) / (fPhotonMom1.E()+fPhotonMom2.E()));
2614  Double_t angle = fPhotonMom1.Angle(fPhotonMom2.Vect());
2615 
2616  Bool_t cutAngle = kFALSE;
2617  if(fUseAngleCut && (angle < fAngleCut || angle > fAngleMaxCut)) cutAngle = kTRUE;
2618 
2619  AliDebug(2,Form("\t ACCEPTED pdg %d: pt %2.2f, phi %2.2f, eta %2.2f",pdg,mesonPt,mesonPhi,mesonYeta));
2620 
2621  if(pdg==111)
2622  {
2623  fhPrimPi0AccE ->Fill(mesonE , GetEventWeight()) ;
2624  fhPrimPi0AccPt ->Fill(mesonPt, GetEventWeight()) ;
2625  fhPrimPi0AccPhi ->Fill(mesonPt, mesonPhi , GetEventWeight()) ;
2626  fhPrimPi0AccY ->Fill(mesonPt, mesonY , GetEventWeight()) ;
2627  fhPrimPi0AccYeta->Fill(mesonPt, mesonYeta, GetEventWeight()) ;
2628 
2630  {
2631  fhPrimPi0AccPtCentrality->Fill(mesonPt, cen, GetEventWeight()) ;
2632  fhPrimPi0AccPtEventPlane->Fill(mesonPt, ep , GetEventWeight()) ;
2633  }
2634 
2635  if(fFillAngleHisto)
2636  {
2637  fhPrimPi0OpeningAngle ->Fill(mesonPt, angle, GetEventWeight());
2638  if(mesonPt > 5)fhPrimPi0OpeningAngleAsym->Fill(asym, angle, GetEventWeight());
2639  fhPrimPi0CosOpeningAngle ->Fill(mesonPt, TMath::Cos(angle), GetEventWeight());
2640  }
2641 
2642  if(fPhotonMom1.Pt() > GetMinPt() && fPhotonMom2.Pt() > GetMinPt() && !cutAngle)
2643  {
2644  fhPrimPi0AccPtPhotonCuts->Fill(mesonPt, GetEventWeight()) ;
2645 
2646  if(fFillAngleHisto)
2647  fhPrimPi0OpeningAnglePhotonCuts->Fill(mesonPt, angle, GetEventWeight());
2648 
2649  if(fNAngleCutBins > 0)
2650  {
2651  Int_t angleBin = -1;
2652  for(Int_t ibin = 0; ibin < fNAngleCutBins; ibin++)
2653  {
2654  if(angle >= fAngleCutBinsArray[ibin] &&
2655  angle < fAngleCutBinsArray[ibin+1]) angleBin = ibin;
2656  }
2657 
2658  if( angleBin >= 0 && angleBin < fNAngleCutBins)
2659  fhPrimPi0AccPtOpAngCuts[angleBin]->Fill(mesonPt,GetEventWeight());
2660  }
2661  }
2662  }
2663  else if(pdg==221)
2664  {
2665  fhPrimEtaAccE ->Fill(mesonE , GetEventWeight()) ;
2666  fhPrimEtaAccPt ->Fill(mesonPt, GetEventWeight()) ;
2667  fhPrimEtaAccPhi ->Fill(mesonPt, mesonPhi , GetEventWeight()) ;
2668  fhPrimEtaAccY ->Fill(mesonPt, mesonY , GetEventWeight()) ;
2669  fhPrimEtaAccYeta->Fill(mesonPt, mesonYeta, GetEventWeight()) ;
2670 
2672  {
2673  fhPrimEtaAccPtCentrality->Fill(mesonPt, cen, GetEventWeight()) ;
2674  fhPrimEtaAccPtEventPlane->Fill(mesonPt, ep , GetEventWeight()) ;
2675  }
2676 
2677  if(fFillAngleHisto)
2678  {
2679  fhPrimEtaOpeningAngle ->Fill(mesonPt, angle, GetEventWeight());
2680  if(mesonPt > 5)fhPrimEtaOpeningAngleAsym->Fill(asym, angle, GetEventWeight());
2681  fhPrimEtaCosOpeningAngle ->Fill(mesonPt, TMath::Cos(angle), GetEventWeight());
2682 
2683  if(fPhotonMom1.Pt() > GetMinPt() && fPhotonMom2.Pt() > GetMinPt() && !cutAngle)
2684  {
2685  if(fUseAngleCut && angle > fAngleCut && angle < fAngleMaxCut)
2686  fhPrimEtaAccPtPhotonCuts->Fill(mesonPt, GetEventWeight()) ;
2687 
2688  if(fFillAngleHisto)
2689  fhPrimEtaOpeningAnglePhotonCuts->Fill(mesonPt, angle, GetEventWeight());
2690 
2691  if(fNAngleCutBins > 0)
2692  {
2693  Int_t angleBin = -1;
2694  for(Int_t ibin = 0; ibin < fNAngleCutBins; ibin++)
2695  {
2696  if(angle >= fAngleCutBinsArray[ibin] &&
2697  angle < fAngleCutBinsArray[ibin+1]) angleBin = ibin;
2698  }
2699 
2700  if( angleBin >= 0 && angleBin < fNAngleCutBins)
2701  fhPrimEtaAccPtOpAngCuts[angleBin]->Fill(mesonPt,GetEventWeight());
2702  }
2703  }
2704  }
2705  }
2706  } // Accepted
2707  } // loop on primaries
2708 }
2709 
2710 //________________________________________________
2712 //________________________________________________
2714 {
2715  // Get pTArm and AlphaArm
2716  Float_t momentumSquaredMother = fMCPrimMesonMom.P()*fMCPrimMesonMom.P();
2717  Float_t momentumDaughter1AlongMother = 0.;
2718  Float_t momentumDaughter2AlongMother = 0.;
2719 
2720  if (momentumSquaredMother > 0.)
2721  {
2722  momentumDaughter1AlongMother = (fPhotonMom1.Px()*fMCPrimMesonMom.Px() + fPhotonMom1.Py()*fMCPrimMesonMom.Py()+ fPhotonMom1.Pz()*fMCPrimMesonMom.Pz()) / sqrt(momentumSquaredMother);
2723  momentumDaughter2AlongMother = (fPhotonMom2.Px()*fMCPrimMesonMom.Px() + fPhotonMom2.Py()*fMCPrimMesonMom.Py()+ fPhotonMom2.Pz()*fMCPrimMesonMom.Pz()) / sqrt(momentumSquaredMother);
2724  }
2725 
2726  Float_t momentumSquaredDaughter1 = fPhotonMom1.P()*fPhotonMom1.P();
2727  Float_t ptArmSquared = momentumSquaredDaughter1 - momentumDaughter1AlongMother*momentumDaughter1AlongMother;
2728 
2729  Float_t pTArm = 0.;
2730  if (ptArmSquared > 0.)
2731  pTArm = sqrt(ptArmSquared);
2732 
2733  Float_t alphaArm = 0.;
2734  if(momentumDaughter1AlongMother + momentumDaughter2AlongMother > 0)
2735  alphaArm = (momentumDaughter1AlongMother -momentumDaughter2AlongMother) / (momentumDaughter1AlongMother + momentumDaughter2AlongMother);
2736 
2738  fPhotonMom1Boost.Boost(-fMCPrimMesonMom.BoostVector());
2739  Float_t cosThStar=TMath::Cos(fPhotonMom1Boost.Vect().Angle(fMCPrimMesonMom.Vect()));
2740 
2741  Float_t en = fMCPrimMesonMom.Energy();
2742  Int_t ebin = -1;
2743  if(en > 8 && en <= 12) ebin = 0;
2744  if(en > 12 && en <= 16) ebin = 1;
2745  if(en > 16 && en <= 20) ebin = 2;
2746  if(en > 20) ebin = 3;
2747  if(ebin < 0 || ebin > 3) return ;
2748 
2749  if(pdg==111)
2750  {
2751  fhCosThStarPrimPi0->Fill(en, cosThStar, GetEventWeight());
2752  fhArmPrimPi0[ebin]->Fill(alphaArm, pTArm , GetEventWeight());
2753  }
2754  else
2755  {
2756  fhCosThStarPrimEta->Fill(en, cosThStar, GetEventWeight());
2757  fhArmPrimEta[ebin]->Fill(alphaArm,pTArm , GetEventWeight());
2758  }
2759 
2760  // if(GetDebug() > 2 )
2761  // {
2762  // Float_t asym = 0;
2763  // if(fPhotonMom1.E()+fPhotonMom2.E() > 0) asym = TMath::Abs(fPhotonMom1.E()-fPhotonMom2.E())/(fPhotonMom1.E()+fPhotonMom2.E());
2764  //
2765  // printf("AliAnaPi0::FillArmenterosThetaStar() - E %f, alphaArm %f, pTArm %f, cos(theta*) %f, asymmetry %1.3f\n",
2766  // en,alphaArm,pTArm,cosThStar,asym);
2767  // }
2768 }
2769 
2770 //_______________________________________________________________________________________
2775 //_______________________________________________________________________________________
2777  Int_t iclus1, Int_t iclus2,
2778  Int_t mctag1, Int_t mctag2,
2779  Float_t pt1, Float_t pt2,
2780  Int_t ncell1, Int_t ncell2,
2781  Double_t mass, Double_t pt, Double_t asym,
2782  Double_t deta, Double_t dphi, Double_t angle)
2783 {
2784  Int_t ancPDG = 0;
2785  Int_t ancStatus = 0;
2786  Int_t ancLabel = GetMCAnalysisUtils()->CheckCommonAncestor(index1, index2,
2787  GetReader(), ancPDG, ancStatus,fMCPrimMesonMom, fMCProdVertex);
2788 
2789  Int_t momindex = -1;
2790  Int_t mompdg = -1;
2791  Int_t momstatus = -1;
2792  Int_t status = -1;
2793  Float_t prodR = -1;
2794  Int_t mcIndex = -1;
2795  Float_t ptPrim = fMCPrimMesonMom.Pt();
2796 
2797  if(ancLabel > -1)
2798  {
2799  AliDebug(1,Form("Common ancestor label %d, pdg %d, name %s, status %d",
2800  ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus));
2801 
2802  if(ancPDG==22)
2803  {//gamma
2804  mcIndex = 0;
2805  }
2806  else if(TMath::Abs(ancPDG)==11)
2807  {//e
2808  mcIndex = 1;
2809  }
2810  else if(ancPDG==111)
2811  {//Pi0
2812  mcIndex = 2;
2813 
2814  fhMCPi0PtTruePtRecRat->Fill(pt, ptPrim/pt, GetEventWeight());
2815  fhMCPi0PtTruePtRecDif->Fill(pt, pt-ptPrim, GetEventWeight());
2816  fhMCPi0PtRecOpenAngle->Fill(pt, angle , GetEventWeight());
2817 
2818  if ( mass < fPi0MassWindow[1] && mass > fPi0MassWindow[0] )
2819  {
2820  fhMCPi0PtTruePtRecRatMassCut->Fill(pt, ptPrim/pt, GetEventWeight());
2821  fhMCPi0PtTruePtRecDifMassCut->Fill(pt, pt-ptPrim, GetEventWeight());
2822  fhMCPi0PtRecOpenAngleMassCut->Fill(pt, angle , GetEventWeight());
2823  }
2824 
2825  if(fMultiCutAnaSim)
2826  {
2827  for(Int_t ipt=0; ipt<fNPtCuts; ipt++)
2828  {
2829  for(Int_t icell=0; icell<fNCellNCuts; icell++)
2830  {
2831  for(Int_t iasym=0; iasym<fNAsymCuts; iasym++)
2832  {
2833  Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
2834  if(pt1 > fPtCuts[ipt] && pt2 > fPtCuts[ipt] &&
2835  pt1 < fPtCutsMax[ipt] && pt2 < fPtCutsMax[ipt] &&
2836  asym < fAsymCuts[iasym] &&
2837  ncell1 >= fCellNCuts[icell] && ncell2 >= fCellNCuts[icell])
2838  {
2839  fhMCPi0MassPtRec [index]->Fill(pt ,mass, GetEventWeight());
2840  fhMCPi0MassPtTrue[index]->Fill(ptPrim,mass, GetEventWeight());
2841 
2842  if ( mass < fPi0MassWindow[1] && mass > fPi0MassWindow[0] )
2843  fhMCPi0PtTruePtRecMassCut[index]->Fill(ptPrim, pt, GetEventWeight());
2844  }//pass the different cuts
2845  }// pid bit cut loop
2846  }// icell loop
2847  }// pt cut loop
2848  }// Multi cut ana sim
2849  else
2850  {
2851  fhMCPi0MassPtTrue [0]->Fill(ptPrim, mass, GetEventWeight());
2852  fhMCPi0MassPtRec [0]->Fill(pt , mass, GetEventWeight());
2853  fhMCPi0PtTruePtRec[0]->Fill(ptPrim, pt, GetEventWeight());
2854 
2855  if ( mass < fPi0MassWindow[1] && mass > fPi0MassWindow[0] )
2856  {
2857  fhMCPi0PtTruePtRecMassCut[0]->Fill(ptPrim, pt, GetEventWeight());
2858 
2859  Float_t momOK = kFALSE;
2860  //Int_t uniqueId = -1;
2861  if(GetReader()->ReadStack())
2862  {
2863  TParticle* ancestor = GetMCStack()->Particle(ancLabel);
2864  status = ancestor->GetStatusCode();
2865  momindex = ancestor->GetFirstMother();
2866  if(momindex >= 0)
2867  {
2868  TParticle* mother = GetMCStack()->Particle(momindex);
2869  mompdg = TMath::Abs(mother->GetPdgCode());
2870  momstatus = mother->GetStatusCode();
2871  prodR = mother->R();
2872  //uniqueId = mother->GetUniqueID();
2873  momOK = kTRUE;
2874  }
2875  }
2876  else
2877  {
2878  TClonesArray * mcparticles = GetReader()->GetAODMCParticles();
2879  AliAODMCParticle* ancestor = (AliAODMCParticle *) mcparticles->At(ancLabel);
2880  status = ancestor->GetStatus();
2881  momindex = ancestor->GetMother();
2882  if(momindex >= 0)
2883  {
2884  AliAODMCParticle* mother = (AliAODMCParticle *) mcparticles->At(momindex);
2885  mompdg = TMath::Abs(mother->GetPdgCode());
2886  momstatus = mother->GetStatus();
2887  prodR = TMath::Sqrt(mother->Xv()*mother->Xv()+mother->Yv()*mother->Yv());
2888  //uniqueId = mother->GetUniqueID();
2889  momOK = kTRUE;
2890  }
2891  }
2892 
2893  if(momOK)
2894  {
2895  // printf("Reco Pi0: pt %2.2f Prod vertex %3.3f, origin pdg %d, origin status %d, origin UI %d\n",
2896  // pt,prodR,mompdg,momstatus,uniqueId);
2897 
2898  fhMCPi0ProdVertex->Fill(pt, prodR , GetEventWeight());
2899  fhMCPi0PtStatus ->Fill(pt, status, GetEventWeight());
2900 
2901  if (momstatus == 21) fhMCPi0PtOrigin->Fill(pt, 0.5, GetEventWeight());//parton
2902  else if(mompdg < 22 ) fhMCPi0PtOrigin->Fill(pt, 1.5, GetEventWeight());//quark
2903  else if(mompdg > 2100 && mompdg < 2210)
2904  fhMCPi0PtOrigin->Fill(pt, 2.5, GetEventWeight());// resonances
2905  else if(mompdg == 221) fhMCPi0PtOrigin->Fill(pt, 8.5, GetEventWeight());//eta
2906  else if(mompdg == 331) fhMCPi0PtOrigin->Fill(pt, 9.5, GetEventWeight());//eta prime
2907  else if(mompdg == 213) fhMCPi0PtOrigin->Fill(pt, 4.5, GetEventWeight());//rho
2908  else if(mompdg == 223) fhMCPi0PtOrigin->Fill(pt, 5.5, GetEventWeight());//omega
2909  else if(mompdg >= 310 && mompdg <= 323)
2910  fhMCPi0PtOrigin->Fill(pt, 6.5, GetEventWeight());//k0S, k+-,k*
2911  else if(mompdg == 130) fhMCPi0PtOrigin->Fill(pt, 6.5, GetEventWeight());//k0L
2912  else if(momstatus == 11 || momstatus == 12 )
2913  fhMCPi0PtOrigin->Fill(pt, 3.5, GetEventWeight());//resonances
2914  else fhMCPi0PtOrigin->Fill(pt, 7.5, GetEventWeight());//other?
2915 
2916  if(status!=11)
2917  {
2918  if (momstatus == 21) fhMCNotResonancePi0PtOrigin->Fill(pt, 0.5, GetEventWeight());//parton
2919  else if(mompdg < 22 ) fhMCNotResonancePi0PtOrigin->Fill(pt, 1.5, GetEventWeight());//quark
2920  else if(mompdg > 2100 && mompdg < 2210)
2921  fhMCNotResonancePi0PtOrigin->Fill(pt, 2.5, GetEventWeight());// resonances
2922  else if(mompdg == 221) fhMCNotResonancePi0PtOrigin->Fill(pt, 8.5, GetEventWeight());//eta
2923  else if(mompdg == 331) fhMCNotResonancePi0PtOrigin->Fill(pt, 9.5, GetEventWeight());//eta prime
2924  else if(mompdg == 213) fhMCNotResonancePi0PtOrigin->Fill(pt, 4.5, GetEventWeight());//rho
2925  else if(mompdg == 223) fhMCNotResonancePi0PtOrigin->Fill(pt, 5.5, GetEventWeight());//omega
2926  else if(mompdg >= 310 && mompdg <= 323)
2927  fhMCNotResonancePi0PtOrigin->Fill(pt, 6.5, GetEventWeight());//k0S, k+-,k*
2928  else if(mompdg == 130) fhMCNotResonancePi0PtOrigin->Fill(pt, 6.5, GetEventWeight());//k0L
2929  else if(momstatus == 11 || momstatus == 12 )
2930  fhMCNotResonancePi0PtOrigin->Fill(pt, 3.5, GetEventWeight());//resonances
2931  else fhMCNotResonancePi0PtOrigin->Fill(pt, 7.5, GetEventWeight());//other?
2932  }
2933  }
2934 
2935  }//pi0 mass region
2936  }
2937 
2938  if(fNAngleCutBins > 0)
2939  {
2940  Int_t angleBin = -1;
2941  for(Int_t ibin = 0; ibin < fNAngleCutBins; ibin++)
2942  {
2943  if(angle >= fAngleCutBinsArray[ibin] &&
2944  angle < fAngleCutBinsArray[ibin+1]) angleBin = ibin;
2945  }
2946 
2947  if( angleBin >= 0 && angleBin < fNAngleCutBins)
2948  fhReOpAngleBinPairClusterMassMCTruePi0[angleBin]->Fill(pt, mass, GetEventWeight());
2949  }
2950 
2951  }
2952  else if(ancPDG==221)
2953  {//Eta
2954  mcIndex = 3;
2955 
2956  fhMCEtaPtTruePtRecRat->Fill(pt, ptPrim/pt, GetEventWeight());
2957  fhMCEtaPtTruePtRecDif->Fill(pt, pt-ptPrim, GetEventWeight());
2958  fhMCEtaPtRecOpenAngle->Fill(pt, angle , GetEventWeight());
2959 
2960  if ( mass < fEtaMassWindow[1] && mass > fEtaMassWindow[0] )
2961  {
2962  fhMCEtaPtTruePtRecRatMassCut->Fill(pt, ptPrim/pt, GetEventWeight());
2963  fhMCEtaPtTruePtRecDifMassCut->Fill(pt, pt-ptPrim, GetEventWeight());
2964  fhMCEtaPtRecOpenAngleMassCut->Fill(pt, angle , GetEventWeight());
2965  }
2966 
2967  if(fMultiCutAnaSim)
2968  {
2969  for(Int_t ipt=0; ipt<fNPtCuts; ipt++)
2970  {
2971  for(Int_t icell=0; icell<fNCellNCuts; icell++)
2972  {
2973  for(Int_t iasym=0; iasym<fNAsymCuts; iasym++)
2974  {
2975  Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
2976  if(pt1 > fPtCuts[ipt] && pt2 > fPtCuts[ipt] &&
2977  pt1 < fPtCutsMax[ipt] && pt2 < fPtCutsMax[ipt] &&
2978  asym < fAsymCuts[iasym] &&
2979  ncell1 >= fCellNCuts[icell] && ncell2 >= fCellNCuts[icell])
2980  {
2981  fhMCEtaMassPtRec [index]->Fill(pt , mass, GetEventWeight());
2982  fhMCEtaMassPtTrue [index]->Fill(ptPrim, mass, GetEventWeight());
2983  fhMCEtaPtTruePtRec[index]->Fill(ptPrim, pt, GetEventWeight());
2984 
2985  if ( mass < fEtaMassWindow[1] && mass > fEtaMassWindow[0] )
2986  fhMCEtaPtTruePtRecMassCut[index]->Fill(ptPrim, pt, GetEventWeight());
2987  }//pass the different cuts
2988  }// pid bit cut loop
2989  }// icell loop
2990  }// pt cut loop
2991  } //Multi cut ana sim
2992  else
2993  {
2994  fhMCEtaMassPtTrue [0]->Fill(ptPrim, mass, GetEventWeight());
2995  fhMCEtaMassPtRec [0]->Fill(pt , mass, GetEventWeight());
2996  fhMCEtaPtTruePtRec[0]->Fill(ptPrim, pt, GetEventWeight());
2997 
2998  if ( mass < fEtaMassWindow[1] && mass > fEtaMassWindow[0] )
2999  {
3000  fhMCEtaPtTruePtRecMassCut[0]->Fill(ptPrim, pt, GetEventWeight());
3001 
3002  Float_t momOK = kFALSE;
3003 
3004  if(GetReader()->ReadStack())
3005  {
3006  TParticle* ancestor = GetMCStack()->Particle(ancLabel);
3007  momindex = ancestor->GetFirstMother();
3008  if(momindex >= 0)
3009  {
3010  TParticle* mother = GetMCStack()->Particle(momindex);
3011  mompdg = TMath::Abs(mother->GetPdgCode());
3012  momstatus = mother->GetStatusCode();
3013  prodR = mother->R();
3014  momOK = kTRUE;
3015  }
3016  }
3017  else
3018  {
3019  TClonesArray * mcparticles = GetReader()->GetAODMCParticles();
3020  AliAODMCParticle* ancestor = (AliAODMCParticle *) mcparticles->At(ancLabel);
3021  momindex = ancestor->GetMother();
3022  if(momindex >= 0)
3023  {
3024  AliAODMCParticle* mother = (AliAODMCParticle *) mcparticles->At(momindex);
3025  mompdg = TMath::Abs(mother->GetPdgCode());
3026  momstatus = mother->GetStatus();
3027  prodR = TMath::Sqrt(mother->Xv()*mother->Xv()+mother->Yv()*mother->Yv());
3028  momOK = kTRUE;
3029  }
3030  }
3031 
3032  if(momOK)
3033  {
3034  fhMCEtaProdVertex->Fill(pt, prodR, GetEventWeight());
3035 
3036  if (momstatus == 21 ) fhMCEtaPtOrigin->Fill(pt, 0.5, GetEventWeight());//parton
3037  else if(mompdg < 22 ) fhMCEtaPtOrigin->Fill(pt, 1.5, GetEventWeight());//quark
3038  else if(mompdg > 2100 && mompdg < 2210)
3039  fhMCEtaPtOrigin->Fill(pt, 2.5, GetEventWeight());//qq resonances
3040  else if(mompdg == 331) fhMCEtaPtOrigin->Fill(pt, 5.5, GetEventWeight());//eta prime
3041  else if(momstatus == 11 || momstatus == 12 )
3042  fhMCEtaPtOrigin->Fill(pt, 3.5, GetEventWeight());//resonances
3043  else fhMCEtaPtOrigin->Fill(pt, 4.5, GetEventWeight());//stable, conversions?
3044  //printf("Other Meson pdg %d, Mother %s, pdg %d, status %d\n",pdg, TDatabasePDG::Instance()->GetParticle(mompdg)->GetName(),mompdg, momstatus );
3045  }
3046 
3047  }// eta mass region
3048  }
3049  if(fNAngleCutBins > 0)
3050  {
3051  Int_t angleBin = -1;
3052  for(Int_t ibin = 0; ibin < fNAngleCutBins; ibin++)
3053  {
3054  if(angle >= fAngleCutBinsArray[ibin] &&
3055  angle < fAngleCutBinsArray[ibin+1]) angleBin = ibin;
3056  }
3057 
3058  if( angleBin >= 0 && angleBin < fNAngleCutBins)
3059  fhReOpAngleBinPairClusterMassMCTrueEta[angleBin]->Fill(pt, mass, GetEventWeight());
3060  }
3061  }
3062  else if(ancPDG==-2212)
3063  {//AProton
3064  mcIndex = 4;
3065  }
3066  else if(ancPDG==-2112)
3067  {//ANeutron
3068  mcIndex = 5;
3069  }
3070  else if(TMath::Abs(ancPDG)==13)
3071  {//muons
3072  mcIndex = 6;
3073  }
3074  else if (TMath::Abs(ancPDG) > 100 && ancLabel > 7)
3075  {
3076  if(ancStatus==1)
3077  {//Stable particles, converted? not decayed resonances
3078  mcIndex = 6;
3079  }
3080  else
3081  {//resonances and other decays, more hadron conversions?
3082  mcIndex = 7;
3083  }
3084  }
3085  else
3086  {//Partons, colliding protons, strings, intermediate corrections
3087  if(ancStatus==11 || ancStatus==12)
3088  {//String fragmentation
3089  mcIndex = 8;
3090  }
3091  else if (ancStatus==21)
3092  {
3093  if(ancLabel < 2)
3094  {//Colliding protons
3095  mcIndex = 11;
3096  }//colliding protons
3097  else if(ancLabel < 6)
3098  {//partonic initial states interactions
3099  mcIndex = 9;
3100  }
3101  else if(ancLabel < 8)
3102  {//Final state partons radiations?
3103  mcIndex = 10;
3104  }
3105  // else {
3106  // printf("AliAnaPi0::FillMCVersusRecDataHistograms() - Check ** Common ancestor label %d, pdg %d, name %s, status %d; \n",
3107  // ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus);
3108  // }
3109  }//status 21
3110  //else {
3111  // printf("AliAnaPi0::FillMCVersusRecDataHistograms() - Check *** Common ancestor label %d, pdg %d, name %s, status %d; \n",
3112  // ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus);
3113  // }
3114  }
3115  }//ancLabel > -1
3116  else
3117  { //ancLabel <= -1
3118  //printf("Not related at all label = %d\n",ancLabel);
3119  AliDebug(1,"Common ancestor not found");
3120 
3121  mcIndex = 12;
3122  }
3123 
3124  if(mcIndex >= 0 && mcIndex < 13)
3125  {
3126  fhMCOrgMass [mcIndex]->Fill(pt, mass, GetEventWeight());
3127  fhMCOrgAsym [mcIndex]->Fill(pt, asym, GetEventWeight());
3128  fhMCOrgDeltaEta[mcIndex]->Fill(pt, deta, GetEventWeight());
3129  fhMCOrgDeltaPhi[mcIndex]->Fill(pt, dphi, GetEventWeight());
3130  }
3131 
3133  {
3134  // Get the original clusters
3135  //
3136  Int_t first = 0;
3137  AliVCluster * cluster1 = FindCluster(GetEMCALClusters(),iclus1,first);
3138  AliVCluster * cluster2 = FindCluster(GetEMCALClusters(),iclus2,iclus1);
3139  if(!cluster2 || !cluster1)
3140  {
3141  AliWarning(Form("Cluster1 %p or Cluster 2 %p not found!",cluster1,cluster2));
3142  return;
3143  }
3144 
3145  // Get the generators names of each cluster and background generator tag
3146  //
3147  TString genName1 = "", genName2 = "", genNameBkg1 = "", genNameBkg2 = "";
3148  Int_t genBkgTag1 = GetCocktailGeneratorBackgroundTag(cluster1, mctag1, genName1, genNameBkg1);
3149  if (genBkgTag1 == -1) return;
3150  else if(genBkgTag1 > 3) printf("Bkg1 generator tag larger than 3; Main %s Bkg %s\n",genName1.Data(),genNameBkg1.Data());
3151 
3152  Int_t genBkgTag2 = GetCocktailGeneratorBackgroundTag(cluster2, mctag2, genName2, genNameBkg2);
3153  if (genBkgTag2 == -1) return;
3154  else if(genBkgTag2 > 3) printf("Bkg2 generator tag larger than 3; Main %s Bkg %s\n",genName2.Data(),genNameBkg2.Data());
3155 
3156  // if the clusters do not come from the same generator exclude
3157  if(genName1!=genName2) return;
3158 
3159  Int_t genType = GetNCocktailGenNamesToCheck()-1;
3160  for(Int_t igen = 1; igen < GetNCocktailGenNamesToCheck(); igen++)
3161  {
3162  if ( genName1.Contains(GetCocktailGenNameToCheck(igen)) )
3163  {
3164  genType = igen;
3165  break;
3166  }
3167  }
3168 
3169  // Fill histograms
3170  //
3171  // assign histogram index number depending on pair combination
3172  Int_t tag = -1;
3173  if ( genBkgTag1 == genBkgTag2 )
3174  {
3175  tag = genBkgTag1; // 0,1,2,3
3176  }
3177  else
3178  {
3179  Int_t genBkgMin = -1;
3180  Int_t genBkgMax = -1;
3181 
3182  if(genBkgTag1 > genBkgTag2)
3183  {
3184  genBkgMax = genBkgTag1;
3185  genBkgMin = genBkgTag2;
3186  }
3187  else
3188  {
3189  genBkgMax = genBkgTag2;
3190  genBkgMin = genBkgTag1;
3191  }
3192 
3193  if ( genBkgMin == 0 )
3194  {
3195  if (genBkgMax == 1 ) tag = 4; // clean+hijing
3196  else if(genBkgMax == 2 ) tag = 5; // clean+not hijing
3197  else if(genBkgMax == 3 ) tag = 6; // clean+multiple
3198  }
3199  else if ( genBkgMin == 1 )
3200  {
3201  if ( genBkgMax == 2 ) tag = 7; // hijing + not hijing
3202  else if ( genBkgMax == 3 ) tag = 8; // hijing + multiple
3203  }
3204  else if ( genBkgMin == 2 ) tag = 9; // not hijing + multiple
3205  }
3206 
3207  if(tag == -1)
3208  {
3209  printf("Combination not found, bkg1 tag %d, bkg2 tag %d\n",genBkgTag1,genBkgTag2);
3210  return;
3211  }
3212 
3213  fhPairGeneratorsBkgMass[genType][tag]->Fill(pt, mass, GetEventWeight());
3214  fhPairGeneratorsBkgMass [0][tag]->Fill(pt, mass, GetEventWeight());
3215 
3216  //
3217  if ( ptPrim < 0.1 || pt < 0.5 ) return;
3218 
3219  Float_t ratio = pt / ptPrim;
3220  Float_t diff = pt - ptPrim;
3221 
3222  if ( mcIndex==2 ) // Pi0
3223  {
3224  fhPairGeneratorsBkgMassMCPi0 [0][tag]->Fill(pt, mass, GetEventWeight());
3225  fhPairGeneratorsBkgMassMCPi0[genType][tag]->Fill(pt, mass, GetEventWeight());
3226 
3227  fhPairGeneratorsBkgEPrimRecoRatioMCPi0[0][tag]->Fill(pt, ratio, GetEventWeight());
3228  fhPairGeneratorsBkgEPrimRecoDiffMCPi0 [0][tag]->Fill(pt, diff, GetEventWeight());
3229 
3230  fhPairGeneratorsBkgEPrimRecoRatioMCPi0[genType][tag]->Fill(pt, ratio, GetEventWeight());
3231  fhPairGeneratorsBkgEPrimRecoDiffMCPi0 [genType][tag]->Fill(pt, diff, GetEventWeight());
3232 
3233  if ( mass < fPi0MassWindow[1] && mass > fPi0MassWindow[0] )
3234  {
3237 
3238  fhPairGeneratorsBkgEPrimRecoRatioMCPi0MassCut[genType][tag]->Fill(pt, ratio, GetEventWeight());
3239  fhPairGeneratorsBkgEPrimRecoDiffMCPi0MassCut [genType][tag]->Fill(pt, diff, GetEventWeight());
3240  }
3241  }
3242  else if( mcIndex==3 ) // Eta
3243  {
3244  fhPairGeneratorsBkgMassMCEta [0][tag]->Fill(pt, mass, GetEventWeight());
3245  fhPairGeneratorsBkgMassMCEta[genType][tag]->Fill(pt, mass, GetEventWeight());
3246 
3247  fhPairGeneratorsBkgEPrimRecoRatioMCEta[0][tag]->Fill(pt, ratio, GetEventWeight());
3248  fhPairGeneratorsBkgEPrimRecoDiffMCEta [0][tag]->Fill(pt, diff, GetEventWeight());
3249 
3250  fhPairGeneratorsBkgEPrimRecoRatioMCEta[genType][tag]->Fill(pt, ratio, GetEventWeight());
3251  fhPairGeneratorsBkgEPrimRecoDiffMCEta [genType][tag]->Fill(pt, diff, GetEventWeight());
3252 
3253  if ( mass < fEtaMassWindow[1] && mass > fEtaMassWindow[0] )
3254  {
3256  fhPairGeneratorsBkgEPrimRecoDiffMCEtaMassCut [0][tag]->Fill(pt, diff, GetEventWeight());
3257 
3258  fhPairGeneratorsBkgEPrimRecoRatioMCEtaMassCut[genType][tag]->Fill(pt, ratio, GetEventWeight());
3259  fhPairGeneratorsBkgEPrimRecoDiffMCEtaMassCut [genType][tag]->Fill(pt, diff, GetEventWeight());
3260  }
3261  }
3262  } // do cluster overlaps from cocktails
3263 
3264  // carefull adding something else here, "returns" can affect
3265 }
3266 
3267 //__________________________________________
3271 //__________________________________________
3273 {
3274  // In case of simulated data, fill acceptance histograms
3275  if(IsDataMC())
3276  {
3278 
3279  if(fFillOnlyMCAcceptanceHisto) return;
3280  }
3281 
3282 //if (GetReader()->GetEventNumber()%10000 == 0)
3283 // printf("--- Event %d ---\n",GetReader()->GetEventNumber());
3284 
3285  if(!GetInputAODBranch())
3286  {
3287  AliFatal(Form("No input aod photons in AOD with name branch < %s >, STOP",GetInputAODName().Data()));
3288  return;
3289  }
3290 
3291  //
3292  // Init some variables
3293  //
3294 
3295  // Analysis within the same detector:
3296  TClonesArray* secondLoopInputData = GetInputAODBranch();
3297 
3298  Int_t nPhot = GetInputAODBranch()->GetEntriesFast() ;
3299  Int_t nPhot2 = nPhot; // second loop
3300  Int_t minEntries = 2;
3301  Int_t last = 1; // last entry in first loop to be removed
3302 
3303  // Combine 2 detectors:
3305  {
3306  // Input from CaloTrackCorr tasks
3307  secondLoopInputData = (TClonesArray *) GetReader()->GetAODBranchList()->FindObject(fOtherDetectorInputName);
3308 
3309  // In case input from PCM or other external source
3310  if(!secondLoopInputData) secondLoopInputData = (TClonesArray *) GetReader()->GetOutputEvent()->FindObject(fOtherDetectorInputName);
3311  if(!secondLoopInputData) secondLoopInputData = (TClonesArray *) GetReader()->GetInputEvent ()->FindObject(fOtherDetectorInputName);
3312 
3313  if(!secondLoopInputData)
3314  {
3315  AliFatal(Form("Input for other detector not found in branch %s!",fOtherDetectorInputName.Data()));
3316  return ; // coverity
3317  }
3318 
3319  nPhot2 = secondLoopInputData->GetEntriesFast() ; // add one since we remove one for the normal case
3320  minEntries = 1;
3321  last = 0;
3322  }
3323 
3324  AliDebug(1,Form("Photon entries %d", nPhot));
3325 
3326  // If less than photon 2 (1) entries in the list, skip this event
3327  if ( nPhot < minEntries )
3328  {
3329  AliDebug(1,Form("nPhotons %d, cent bin %d continue to next event",nPhot, GetEventCentralityBin()));
3330 
3331  if ( GetNCentrBin() > 1 && IsHighMultiplicityAnalysisOn() )
3333 
3334  return ;
3335  }
3336  else if(fPairWithOtherDetector && nPhot2 < minEntries)
3337  {
3338  AliDebug(1,Form("nPhotons %d, cent bin %d continue to next event",nPhot, GetEventCentralityBin()));
3339 
3340  if ( GetNCentrBin() > 1 && IsHighMultiplicityAnalysisOn() )
3342 
3343  return ;
3344  }
3345 
3346  Int_t ncentr = GetNCentrBin();
3347  if(GetNCentrBin()==0) ncentr = 1;
3348 
3349  //Init variables
3350  Int_t module1 = -1;
3351  Int_t module2 = -1;
3352  Double_t vert[] = {0.0, 0.0, 0.0} ; //vertex
3353  Int_t evtIndex1 = 0 ;
3354  Int_t currentEvtIndex = -1;
3355  Int_t curCentrBin = GetEventCentralityBin();
3356  //Int_t curVzBin = GetEventVzBin();
3357  //Int_t curRPBin = GetEventRPBin();
3358  Int_t eventbin = GetEventMixBin();
3359 
3360  if(eventbin > GetNCentrBin()*GetNZvertBin()*GetNRPBin())
3361  {
3362  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",
3364  return;
3365  }
3366 
3367  // Get shower shape information of clusters
3368 // TObjArray *clusters = 0;
3369 // if (GetCalorimeter()==kEMCAL) clusters = GetEMCALClusters();
3370 // else if(GetCalorimeter()==kPHOS ) clusters = GetPHOSClusters() ;
3371 
3372  //---------------------------------
3373  // First loop on photons/clusters
3374  //---------------------------------
3375  for(Int_t i1 = 0; i1 < nPhot-last; i1++)
3376  {
3377  AliAODPWG4Particle * p1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i1)) ;
3378 
3379  // Select photons within a pT range
3380  if ( p1->Pt() < GetMinPt() || p1->Pt() > GetMaxPt() ) continue ;
3381 
3382  //printf("AliAnaPi0::MakeAnalysisFillHistograms() : cluster1 id %d/%d\n",i1,nPhot-1);
3383 
3384  // get the event index in the mixed buffer where the photon comes from
3385  // in case of mixing with analysis frame, not own mixing
3386  evtIndex1 = GetEventIndex(p1, vert) ;
3387  if ( evtIndex1 == -1 )
3388  return ;
3389  if ( evtIndex1 == -2 )
3390  continue ;
3391 
3392  // Only effective in case of mixed event frame
3393  if(TMath::Abs(vert[2]) > GetZvertexCut()) continue ; //vertex cut
3394 
3395  if (evtIndex1 != currentEvtIndex)
3396  {
3397  // Fill event bin info
3398  if( DoOwnMix() ) fhEventBin->Fill(eventbin, GetEventWeight()) ;
3399 
3401  {
3402  fhCentrality->Fill(curCentrBin, GetEventWeight());
3403 
3404  if( GetEventPlane() )
3405  fhEventPlaneResolution->Fill(curCentrBin, TMath::Cos(2.*GetEventPlane()->GetQsubRes()), GetEventWeight());
3406  }
3407 
3408  currentEvtIndex = evtIndex1 ;
3409  }
3410 
3411  //printf("AliAnaPi0::MakeAnalysisFillHistograms(): Photon 1 Evt %d Vertex : %f,%f,%f\n",evtIndex1, GetVertex(evtIndex1)[0] ,GetVertex(evtIndex1)[1],GetVertex(evtIndex1)[2]);
3412 
3413  //Get the momentum of this cluster
3414  fPhotonMom1.SetPxPyPzE(p1->Px(),p1->Py(),p1->Pz(),p1->E());
3415 
3416  //Get (Super)Module number of this cluster
3417  module1 = p1->GetSModNumber();// GetModuleNumber(p1);
3418 
3419  //------------------------------------------
3420  // Recover original cluster
3421  // Declare variables for absid col-row identification of pair
3422  // Also fill col-row, eta-phi histograms depending on pt bin
3423 
3424  Int_t iclus1 = -1, iclus2 = -1 ;
3425  Float_t maxCellFraction1 = 0, maxCellFraction2 = 0;
3426  Int_t absIdMax1 = -1, absIdMax2 = -1;
3427  Int_t icol1 = -1, icol2 = -1, icolAbs1 = -1, icolAbs2 = -1;
3428  Int_t irow1 = -1, irow2 = -1, irowAbs1 = -1, irowAbs2 = -1;
3429  Int_t iRCU1 = -1, iRCU2 = -1;
3430 
3432  {
3433  AliVCluster * cluster1 = FindCluster(GetEMCALClusters(),p1->GetCaloLabel(0),iclus1);
3434  if(!cluster1) AliWarning("Cluster1 not found!");
3435 
3436  absIdMax1 = GetCaloUtils()->GetMaxEnergyCell(GetEMCALCells(),cluster1,maxCellFraction1);
3437 
3438  GetModuleNumberCellIndexesAbsCaloMap(absIdMax1,GetCalorimeter(), icol1, irow1, iRCU1, icolAbs1, irowAbs1);
3439 
3440  if(fMultiCutAnaAcc)
3441  {
3442  for(Int_t ipt = 0; ipt < fNPtCuts; ipt++)
3443  {
3444  if( p1->Pt() > fPtCuts[ipt] && p1->Pt() < fPtCuts[ipt+1] )
3445  {
3446  fhPtBinClusterEtaPhi[ipt]->Fill(p1->Eta(),GetPhi(p1->Phi()),GetEventWeight()) ;
3447 
3448  fhPtBinClusterColRow[ipt]->Fill(icolAbs1,irowAbs1,GetEventWeight()) ;
3449  }
3450  }
3451  }
3452  }
3453 
3454  //---------------------------------
3455  // Second loop on photons/clusters
3456  //---------------------------------
3457  Int_t first = i1+1;
3458  if(fPairWithOtherDetector) first = 0;
3459 
3460  for(Int_t i2 = first; i2 < nPhot2; i2++)
3461  {
3462  //AliAODPWG4Particle * p2 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i2)) ;
3463  AliAODPWG4Particle * p2 = (AliAODPWG4Particle*) (secondLoopInputData->At(i2)) ;
3464 
3465  // Select photons within a pT range
3466  if ( p2->Pt() < GetMinPt() || p2->Pt() > GetMaxPt() ) continue ;
3467 
3468  //printf("AliAnaPi0::MakeAnalysisFillHistograms() : cluster2 i %d/%d\n",i2,nPhot);
3469 
3470  //In case of mixing frame, check we are not in the same event as the first cluster
3471  Int_t evtIndex2 = GetEventIndex(p2, vert) ;
3472  if ( evtIndex2 == -1 )
3473  return ;
3474  if ( evtIndex2 == -2 )
3475  continue ;
3476  if (GetMixedEvent() && (evtIndex1 == evtIndex2))
3477  continue ;
3478 
3479 // //------------------------------------------
3480 // // Recover original cluster
3481 // Int_t iclus2 = -1;
3482 // AliVCluster * cluster2 = FindCluster(clusters,p2->GetCaloLabel(0),iclus2,iclus1+1);
3483 // // start new loop from iclus1+1 to gain some time
3484 // if(!cluster2) AliWarning("Cluster2 not found!");
3485 //
3486 // // Get the TOF,l0 and ncells from the clusters
3487 // Float_t tof1 = -1;
3488 // Float_t l01 = -1;
3489 // Int_t ncell1 = 0;
3490 // if(cluster1)
3491 // {
3492 // tof1 = cluster1->GetTOF()*1e9;
3493 // l01 = cluster1->GetM02();
3494 // ncell1 = cluster1->GetNCells();
3495 // //printf("cluster1: E %2.2f (%2.2f), l0 %2.2f, tof %2.2f\n",cluster1->E(),p1->E(),l01,tof1);
3496 // }
3497 // //else printf("cluster1 not available: calo label %d / %d, cluster ID %d\n",
3498 // // p1->GetCaloLabel(0),(GetReader()->GetInputEvent())->GetNumberOfCaloClusters()-1,cluster1->GetID());
3499 //
3500 // Float_t tof2 = -1;
3501 // Float_t l02 = -1;
3502 // Int_t ncell2 = 0;
3503 // if(cluster2)
3504 // {
3505 // tof2 = cluster2->GetTOF()*1e9;
3506 // l02 = cluster2->GetM02();
3507 // ncell2 = cluster2->GetNCells();
3508 // //printf("cluster2: E %2.2f (%2.2f), l0 %2.2f, tof %2.2f\n",cluster2->E(),p2->E(),l02,tof2);
3509 // }
3510 // //else printf("cluster2 not available: calo label %d / %d, cluster ID %d\n",
3511 // // p2->GetCaloLabel(0),(GetReader()->GetInputEvent())->GetNumberOfCaloClusters()-1,cluster2->GetID());
3512 //
3513 // if(cluster1 && cluster2)
3514 // {
3515 // Double_t t12diff = tof1-tof2;
3516 // if(TMath::Abs(t12diff) > GetPairTimeCut()) continue;
3517 // }
3518 
3519  Float_t tof1 = p1->GetTime();
3520  Float_t l01 = p1->GetM02();
3521  Int_t ncell1 = p1->GetNCells();
3522  //printf("cluster1: E %2.2f, l0 %2.2f, tof %2.2f\n",p1->E(),l01,tof1);
3523 
3524  Float_t tof2 = p2->GetTime();
3525  Float_t l02 = p2->GetM02();
3526  Int_t ncell2 = p2->GetNCells();
3527  //printf("cluster2: E %2.2f, l0 %2.2f, tof %2.2f\n",p2->E(),l02,tof2);
3528 
3529  Double_t t12diff = tof1-tof2;
3530  fhEPairDiffTime->Fill((fPhotonMom1 + fPhotonMom2).Pt(), t12diff, GetEventWeight());
3531  if(TMath::Abs(t12diff) > GetPairTimeCut()) continue;
3532 
3533  //------------------------------------------
3534 
3535  //printf("AliAnaPi0::MakeAnalysisFillHistograms(): Photon 2 Evt %d Vertex : %f,%f,%f\n",evtIndex2, GetVertex(evtIndex2)[0] ,GetVertex(evtIndex2)[1],GetVertex(evtIndex2)[2]);
3536 
3537  // Get the momentum of this cluster
3538  fPhotonMom2.SetPxPyPzE(p2->Px(),p2->Py(),p2->Pz(),p2->E());
3539 
3540  // Get module number
3541  module2 = p2->GetSModNumber(); //GetModuleNumber(p2);
3542 
3543  //---------------------------------
3544  // Get pair kinematics
3545  //---------------------------------
3546  Double_t m = (fPhotonMom1 + fPhotonMom2).M() ;
3547  Double_t pt = (fPhotonMom1 + fPhotonMom2).Pt();
3548  Double_t deta = fPhotonMom1.Eta() - fPhotonMom2.Eta();
3549  Double_t dphi = fPhotonMom1.Phi() - fPhotonMom2.Phi();
3550  Double_t a = TMath::Abs(p1->E()-p2->E())/(p1->E()+p2->E()) ;
3551 
3552  AliDebug(2,Form("E: fPhotonMom1 %f, fPhotonMom2 %f; Pair: pT %f, mass %f, a %f", p1->E(), p2->E(), (fPhotonMom1 + fPhotonMom2).E(),m,a));
3553 
3554  //--------------------------------
3555  // Opening angle selection
3556  //--------------------------------
3557  // Check if opening angle is too large or too small compared to what is expected
3558  Double_t angle = fPhotonMom1.Angle(fPhotonMom2.Vect());
3560  {
3561  AliDebug(2,Form("Real pair angle %f (deg) not in E %f window",RadToDeg(angle), (fPhotonMom1+fPhotonMom2).E()));
3562  continue;
3563  }
3564 
3565  if(fUseAngleCut && (angle < fAngleCut || angle > fAngleMaxCut))
3566  {
3567  AliDebug(2,Form("Real pair cut %f < angle %f < cut %f (deg)",RadToDeg(fAngleCut), RadToDeg(angle), RadToDeg(fAngleMaxCut)));
3568  continue;
3569  }
3570 
3571  //-------------------------------------------------------------------------------------------------
3572  // Fill module dependent histograms, put a cut on assymmetry on the first available cut in the array
3573  //-------------------------------------------------------------------------------------------------
3574  if(a < fAsymCuts[0] && fFillSMCombinations)
3575  {
3577  {
3578  if(module1==module2 && module1 >=0 && module1<fNModules)
3579  {
3580  fhReMod[module1]->Fill(pt, m, GetEventWeight()) ;
3581  if(fFillAngleHisto) fhRealOpeningAnglePerSM[module1]->Fill(pt, angle, GetEventWeight());
3582  }
3583 
3584  if (GetCalorimeter() == kEMCAL )
3585  {
3586  // Same sector
3587  Int_t j=0;
3588  for(Int_t i = 0; i < fNModules/2; i++)
3589  {
3590  j=2*i;
3591  if((module1==j && module2==j+1) || (module1==j+1 && module2==j)) fhReSameSectorEMCALMod[i]->Fill(pt, m, GetEventWeight()) ;
3592  }
3593 
3594  // Same side
3595  for(Int_t i = 0; i < fNModules-2; i++){
3596  if((module1==i && module2==i+2) || (module1==i+2 && module2==i)) fhReSameSideEMCALMod[i]->Fill(pt, m, GetEventWeight());
3597  }
3598  } // EMCAL
3599  else
3600  { // PHOS
3601  if((module1==0 && module2==1) || (module1==1 && module2==0)) fhReDiffPHOSMod[0]->Fill(pt, m, GetEventWeight()) ;
3602  if((module1==0 && module2==2) || (module1==2 && module2==0)) fhReDiffPHOSMod[1]->Fill(pt, m, GetEventWeight()) ;
3603  if((module1==1 && module2==2) || (module1==2 && module2==1)) fhReDiffPHOSMod[2]->Fill(pt, m, GetEventWeight()) ;
3604  } // PHOS
3605  }
3606  else
3607  {
3608  Float_t phi1 = GetPhi(fPhotonMom1.Phi());
3609  Float_t phi2 = GetPhi(fPhotonMom2.Phi());
3610  Bool_t etaside = 0;
3611  if( (p1->GetDetectorTag()==kEMCAL && fPhotonMom1.Eta() < 0)
3612  || (p2->GetDetectorTag()==kEMCAL && fPhotonMom2.Eta() < 0)) etaside = 1;
3613 
3614  if ( phi1 > DegToRad(260) && phi2 > DegToRad(260) && phi1 < DegToRad(280) && phi2 < DegToRad(280)) fhReSameSectorDCALPHOSMod[0+etaside]->Fill(pt, m, GetEventWeight());
3615  else if ( phi1 > DegToRad(280) && phi2 > DegToRad(280) && phi1 < DegToRad(300) && phi2 < DegToRad(300)) fhReSameSectorDCALPHOSMod[2+etaside]->Fill(pt, m, GetEventWeight());
3616  else if ( phi1 > DegToRad(300) && phi2 > DegToRad(300) && phi1 < DegToRad(320) && phi2 < DegToRad(320)) fhReSameSectorDCALPHOSMod[4+etaside]->Fill(pt, m, GetEventWeight());
3617  else if ( (phi1 > DegToRad(260) && phi2 > DegToRad(280) && phi1 < DegToRad(280) && phi2 < DegToRad(300))
3618  || (phi1 > DegToRad(280) && phi2 > DegToRad(260) && phi1 < DegToRad(300) && phi2 < DegToRad(280))) fhReDiffSectorDCALPHOSMod[0+etaside]->Fill(pt, m, GetEventWeight());
3619  else if ( (phi1 > DegToRad(280) && phi2 > DegToRad(300) && phi1 < DegToRad(300) && phi2 < DegToRad(320))
3620  || (phi1 > DegToRad(300) && phi2 > DegToRad(280) && phi1 < DegToRad(320) && phi2 < DegToRad(300))) fhReDiffSectorDCALPHOSMod[2+etaside]->Fill(pt, m, GetEventWeight());
3621  else if ( (phi1 > DegToRad(260) && phi2 > DegToRad(300) && phi1 < DegToRad(280) && phi2 < DegToRad(320))
3622  || (phi1 > DegToRad(300) && phi2 > DegToRad(260) && phi1 < DegToRad(320) && phi2 < DegToRad(280))) fhReDiffSectorDCALPHOSMod[4+etaside]->Fill(pt, m, GetEventWeight());
3623  else fhReDiffSectorDCALPHOSMod[6+etaside]->Fill(pt, m, GetEventWeight());
3624  }
3625  } // Fill SM combinations
3626 
3627  // In case we want only pairs in same (super) module, check their origin.
3628  Bool_t ok = kTRUE;
3629 
3630  if(fSameSM)
3631  {
3633  {
3634  if(module1!=module2) ok=kFALSE;
3635  }
3636  else // PHOS and DCal in same sector
3637  {
3638  Float_t phi1 = GetPhi(fPhotonMom1.Phi());
3639  Float_t phi2 = GetPhi(fPhotonMom2.Phi());
3640  ok=kFALSE;
3641  if ( phi1 > DegToRad(260) && phi2 > DegToRad(260) && phi1 < DegToRad(280) && phi2 < DegToRad(280)) ok = kTRUE;
3642  else if ( phi1 > DegToRad(280) && phi2 > DegToRad(280) && phi1 < DegToRad(300) && phi2 < DegToRad(300)) ok = kTRUE;
3643  else if ( phi1 > DegToRad(300) && phi2 > DegToRad(300) && phi1 < DegToRad(320) && phi2 < DegToRad(320)) ok = kTRUE;
3644  }
3645  } // Pair only in same SM
3646 
3647  if(!ok) continue;
3648 
3649  //
3650  // Fill histograms with selected cluster pairs
3651  //
3652 
3653  // Check if one of the clusters comes from a conversion
3654  if(fCheckConversion)
3655  {
3656  if (p1->IsTagged() && p2->IsTagged()) fhReConv2->Fill(pt, m, GetEventWeight());
3657  else if(p1->IsTagged() || p2->IsTagged()) fhReConv ->Fill(pt, m, GetEventWeight());
3658  }
3659 
3660  // Fill shower shape cut histograms
3662  {
3663  if ( l01 > 0.01 && l01 < 0.4 &&
3664  l02 > 0.01 && l02 < 0.4 ) fhReSS[0]->Fill(pt, m, GetEventWeight()); // Tight
3665  else if( l01 > 0.4 && l02 > 0.4 ) fhReSS[1]->Fill(pt, m, GetEventWeight()); // Loose
3666  else if( l01 > 0.01 && l01 < 0.4 && l02 > 0.4 ) fhReSS[2]->Fill(pt, m, GetEventWeight()); // Both
3667  else if( l02 > 0.01 && l02 < 0.4 && l01 > 0.4 ) fhReSS[2]->Fill(pt, m, GetEventWeight()); // Both
3668  }
3669 
3670  //
3671  // Main invariant mass histograms.
3672  // Fill histograms for different bad channel distance, centrality, assymmetry cut and pid bit
3673  //
3674  for(Int_t ipid=0; ipid<fNPIDBits; ipid++)
3675  {
3676  if((p1->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton)) && (p2->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton)))
3677  {
3678  for(Int_t iasym=0; iasym < fNAsymCuts; iasym++)
3679  {
3680  if(a < fAsymCuts[iasym])
3681  {
3682  Int_t index = ((curCentrBin*fNPIDBits)+ipid)*fNAsymCuts + iasym;
3683  //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);
3684 
3685  if(index < 0 || index >= ncentr*fNPIDBits*fNAsymCuts) continue ;
3686 
3687  fhRe1 [index]->Fill(pt, m, GetEventWeight());
3688  if(fMakeInvPtPlots)fhReInvPt1[index]->Fill(pt, m, 1./pt * GetEventWeight()) ;
3689 
3690  if(fFillBadDistHisto)
3691  {
3692  if(p1->DistToBad()>0 && p2->DistToBad()>0)
3693  {
3694  fhRe2 [index]->Fill(pt, m, GetEventWeight()) ;
3695  if(fMakeInvPtPlots)fhReInvPt2[index]->Fill(pt, m, 1./pt * GetEventWeight()) ;
3696 
3697  if(p1->DistToBad()>1 && p2->DistToBad()>1)
3698  {
3699  fhRe3 [index]->Fill(pt, m, GetEventWeight()) ;
3700  if(fMakeInvPtPlots)fhReInvPt3[index]->Fill(pt, m, 1./pt * GetEventWeight()) ;
3701  }// bad 3
3702  }// bad2
3703  }// Fill bad dist histos
3704  }//assymetry cut
3705  }// asymmetry cut loop
3706  }// bad 1
3707  }// pid bit loop
3708 
3709  //
3710  // Fill histograms with opening angle
3711  if(fFillAngleHisto)
3712  {
3713  fhRealOpeningAngle ->Fill(pt, angle, GetEventWeight());
3714  fhRealCosOpeningAngle->Fill(pt, TMath::Cos(angle), GetEventWeight());
3715  }
3716 
3717  //
3718  // Fill histograms for different opening angle bins
3720  {
3721  Int_t angleBin = -1;
3722  for(Int_t ibin = 0; ibin < fNAngleCutBins; ibin++)
3723  {
3724  if(angle >= fAngleCutBinsArray[ibin] &&
3725  angle < fAngleCutBinsArray[ibin+1]) angleBin = ibin;
3726  }
3727 
3728  if( angleBin >= 0 && angleBin < fNAngleCutBins)
3729  {
3730  Float_t e1 = fPhotonMom1.E();
3731  Float_t e2 = fPhotonMom2.E();
3732 
3733  Float_t t1 = tof1;
3734  Float_t t2 = tof2;
3735 
3736  Int_t nc1 = ncell1;
3737  Int_t nc2 = ncell2;
3738 
3739  Float_t eta1 = fPhotonMom1.Eta();
3740  Float_t eta2 = fPhotonMom2.Eta();
3741 
3742  Float_t phi1 = GetPhi(fPhotonMom1.Phi());
3743  Float_t phi2 = GetPhi(fPhotonMom2.Phi());
3744 
3745  Int_t mod1 = module1;
3746  Int_t mod2 = module2;
3747 
3748  // Recover original cluster
3749  AliVCluster * cluster2 = FindCluster(GetEMCALClusters(),p2->GetCaloLabel(0),iclus2);
3750  if(!cluster2) AliWarning("Cluster2 not found!");
3751 
3752  absIdMax2 = GetCaloUtils()->GetMaxEnergyCell(GetEMCALCells(),cluster2,maxCellFraction2);
3753 
3754  if(e2 > e1)
3755  {
3756  e1 = fPhotonMom2.E();
3757  e2 = fPhotonMom1.E();
3758 
3759  t1 = tof2;
3760  t2 = tof1;
3761 
3762  nc1 = ncell2;
3763  nc2 = ncell1;
3764 
3765  eta1 = fPhotonMom2.Eta();
3766  eta2 = fPhotonMom1.Eta();
3767 
3768  phi1 = GetPhi(fPhotonMom2.Phi());
3769  phi2 = GetPhi(fPhotonMom1.Phi());
3770 
3771  mod1 = module2;
3772  mod2 = module1;
3773 
3774  Int_t tmp = absIdMax2;
3775  absIdMax2 = absIdMax1;
3776  absIdMax1 = tmp;
3777  }
3778 
3779  fhReOpAngleBinMinClusterEPerSM[angleBin]->Fill(e2,mod2,GetEventWeight()) ;
3780  fhReOpAngleBinMaxClusterEPerSM[angleBin]->Fill(e1,mod1,GetEventWeight()) ;
3781 
3782  fhReOpAngleBinMinClusterTimePerSM[angleBin]->Fill(t2,mod2,GetEventWeight()) ;
3783  fhReOpAngleBinMaxClusterTimePerSM[angleBin]->Fill(t1,mod1,GetEventWeight()) ;
3784 
3785  fhReOpAngleBinMinClusterNCellPerSM[angleBin]->Fill(nc2,mod2,GetEventWeight()) ;
3786  fhReOpAngleBinMaxClusterNCellPerSM[angleBin]->Fill(nc1,mod1,GetEventWeight()) ;
3787 
3788  fhReOpAngleBinPairClusterMass[angleBin]->Fill(pt,m,GetEventWeight()) ;
3789  if(mod2 == mod1) fhReOpAngleBinPairClusterMassPerSM[angleBin]->Fill(m,mod1,GetEventWeight()) ;
3790 
3791  if(e1 > 0.01) fhReOpAngleBinPairClusterRatioPerSM[angleBin]->Fill(e2/e1,mod1,GetEventWeight()) ;
3792 
3793  fhReOpAngleBinMinClusterEtaPhi[angleBin]->Fill(eta2,phi2,GetEventWeight()) ;
3794  fhReOpAngleBinMaxClusterEtaPhi[angleBin]->Fill(eta1,phi1,GetEventWeight()) ;
3795 
3796  GetModuleNumberCellIndexesAbsCaloMap(absIdMax2,GetCalorimeter(), icol2, irow2, iRCU2, icolAbs2, irowAbs2);
3797 
3798  //fhReOpAngleBinPairClusterAbsIdMaxCell[angleBin]->Fill(absIdMax1,absIdMax2,GetEventWeight());
3799 
3800  fhReOpAngleBinMinClusterColRow[angleBin]->Fill(icolAbs2,irowAbs2,GetEventWeight()) ;
3801  fhReOpAngleBinMaxClusterColRow[angleBin]->Fill(icolAbs1,irowAbs1,GetEventWeight()) ;
3802  }
3803  } // fFillOpAngleHisto
3804 
3805 
3806  // Fill histograms with pair assymmetry
3808  {
3809  fhRePtAsym->Fill(pt, a, GetEventWeight());
3810  if ( m > fPi0MassWindow[0] && m < fPi0MassWindow[1] ) fhRePtAsymPi0->Fill(pt, a, GetEventWeight());
3811  if ( m > fEtaMassWindow[0] && m < fEtaMassWindow[1] ) fhRePtAsymEta->Fill(pt, a, GetEventWeight());
3812  }
3813 
3814  // Check cell time content in cluster
3816  {
3817  if ( p1->GetFiducialArea() == 0 && p2->GetFiducialArea() == 0 )
3819 
3820  else if ( p1->GetFiducialArea() != 0 && p2->GetFiducialArea() != 0 )
3822  }
3823 
3824  //---------
3825  // MC data
3826  //---------
3827  // Do some MC checks on the origin of the pair, is there any common ancestor and if there is one, who?
3828  if(IsDataMC())
3829  {
3830  if(GetMCAnalysisUtils()->CheckTagBit(p1->GetTag(),AliMCAnalysisUtils::kMCConversion) &&
3832  {
3833  fhReMCFromConversion->Fill(pt, m, GetEventWeight());
3834  }
3835  else if(!GetMCAnalysisUtils()->CheckTagBit(p1->GetTag(),AliMCAnalysisUtils::kMCConversion) &&
3837  {
3838  fhReMCFromNotConversion->Fill(pt, m, GetEventWeight());
3839  }
3840  else
3841  {
3842  fhReMCFromMixConversion->Fill(pt, m, GetEventWeight());
3843  }
3844 
3845  if(fFillOriginHisto)
3846  FillMCVersusRecDataHistograms(p1->GetLabel(), p2->GetLabel(),
3847  p1->GetCaloLabel(0), p2->GetCaloLabel(0),
3848  p1->GetTag(),p2->GetTag(),
3849  p1->Pt(), p2->Pt(),
3850  ncell1, ncell2, m, pt, a,deta, dphi, angle);
3851  }
3852 
3853  //-----------------------
3854  // Multi cuts analysis
3855  //-----------------------
3856  if(fMultiCutAna)
3857  {
3858  // Several pt,ncell and asymmetry cuts
3859  for(Int_t ipt = 0; ipt < fNPtCuts; ipt++)
3860  {
3861  for(Int_t icell = 0; icell < fNCellNCuts; icell++)
3862  {
3863  for(Int_t iasym = 0; iasym < fNAsymCuts; iasym++)
3864  {
3865  Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
3866  if(p1->Pt() > fPtCuts[ipt] && p2->Pt() > fPtCuts[ipt] &&
3867  p1->Pt() < fPtCutsMax[ipt] && p2->Pt() < fPtCutsMax[ipt] &&
3868  a < fAsymCuts[iasym] &&
3869  ncell1 >= fCellNCuts[icell] && ncell2 >= fCellNCuts[icell])
3870  {
3871  fhRePtNCellAsymCuts[index]->Fill(pt, m, GetEventWeight()) ;
3872  if(fFillAngleHisto) fhRePtNCellAsymCutsOpAngle[index]->Fill(pt, angle, GetEventWeight()) ;
3873 
3874  if(fFillSMCombinations && module1==module2)
3875  {
3876  fhRePtNCellAsymCutsSM[module1][index]->Fill(pt, m, GetEventWeight()) ;
3877  if(fFillAngleHisto) fhRePtNCellAsymCutsSMOpAngle[module1][index]->Fill(pt, angle, GetEventWeight()) ;
3878  }
3879  }
3880  }// pid bit cut loop
3881  }// icell loop
3882  }// pt cut loop
3883  }// multiple cuts analysis
3884 
3885  }// second same event particle
3886  }// first cluster
3887 
3888  //-------------------------------------------------------------
3889  // Mixing
3890  //-------------------------------------------------------------
3891  if(DoOwnMix())
3892  {
3893  // Recover events in with same characteristics as the current event
3894 
3895  // Check that the bin exists, if not (bad determination of RP, centrality or vz bin) do nothing
3896  if(eventbin < 0) return ;
3897 
3898  TList * evMixList=fEventsList[eventbin] ;
3899 
3900  if(!evMixList)
3901  {
3902  AliWarning(Form("Mix event list not available, bin %d",eventbin));
3903  return;
3904  }
3905 
3906  Int_t nMixed = evMixList->GetSize() ;
3907  for(Int_t ii=0; ii<nMixed; ii++)
3908  {
3909  TClonesArray* ev2= (TClonesArray*) (evMixList->At(ii));
3910  Int_t nPhot2=ev2->GetEntriesFast() ;
3911  Double_t m = -999;
3912  AliDebug(1,Form("Mixed event %d photon entries %d, centrality bin %d",ii, nPhot2, GetEventCentralityBin()));
3913 
3914  fhEventMixBin->Fill(eventbin, GetEventWeight()) ;
3915 
3916  //---------------------------------
3917  // First loop on photons/clusters
3918  //---------------------------------
3919  for(Int_t i1 = 0; i1 < nPhot; i1++)
3920  {
3921  AliAODPWG4Particle * p1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i1)) ;
3922 
3923  // Select photons within a pT range
3924  if ( p1->Pt() < GetMinPt() || p1->Pt() > GetMaxPt() ) continue ;
3925 
3926  // Not sure why this line is here
3927  //if(fSameSM && GetModuleNumber(p1)!=module1) continue;
3928 
3929  //Get kinematics of cluster and (super) module of this cluster
3930  fPhotonMom1.SetPxPyPzE(p1->Px(),p1->Py(),p1->Pz(),p1->E());
3931  module1 = GetModuleNumber(p1);
3932 
3933  //---------------------------------
3934  // Second loop on other mixed event photons/clusters
3935  //---------------------------------
3936  for(Int_t i2 = 0; i2 < nPhot2; i2++)
3937  {
3938  AliAODPWG4Particle * p2 = (AliAODPWG4Particle*) (ev2->At(i2)) ;
3939 
3940  // Select photons within a pT range
3941  if ( p2->Pt() < GetMinPt() || p2->Pt() > GetMaxPt() ) continue ;
3942 
3943  // Get kinematics of second cluster and calculate those of the pair
3944  fPhotonMom2.SetPxPyPzE(p2->Px(),p2->Py(),p2->Pz(),p2->E());
3945  m = (fPhotonMom1+fPhotonMom2).M() ;
3946  Double_t pt = (fPhotonMom1 + fPhotonMom2).Pt();
3947  Double_t a = TMath::Abs(p1->E()-p2->E())/(p1->E()+p2->E()) ;
3948 
3949  // Check if opening angle is too large or too small compared to what is expected
3950  Double_t angle = fPhotonMom1.Angle(fPhotonMom2.Vect());
3952  {
3953  AliDebug(2,Form("Mix pair angle %f (deg) not in E %f window",RadToDeg(angle), (fPhotonMom1+fPhotonMom2).E()));
3954  continue;
3955  }
3956 
3957  if(fUseAngleCut && (angle < fAngleCut || angle > fAngleMaxCut))
3958  {
3959  AliDebug(2,Form("Mix pair cut %f < angle %f < cut %f (deg)",RadToDeg(fAngleCut),RadToDeg(angle),RadToDeg(fAngleMaxCut)));
3960  continue;
3961  }
3962 
3963  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));
3964 
3965  // In case we want only pairs in same (super) module, check their origin.
3966  module2 = GetModuleNumber(p2);
3967 
3968  //-------------------------------------------------------------------------------------------------
3969  // Fill module dependent histograms, put a cut on assymmetry on the first available cut in the array
3970  //-------------------------------------------------------------------------------------------------
3971  if(a < fAsymCuts[0] && fFillSMCombinations)
3972  {
3974  {
3975  if(module1==module2 && module1 >=0 && module1<fNModules)
3976  {
3977  fhMiMod[module1]->Fill(pt, m, GetEventWeight()) ;
3978  if(fFillAngleHisto) fhMixedOpeningAnglePerSM[module1]->Fill(pt, angle, GetEventWeight());
3979  }
3980 
3981  if(GetCalorimeter()==kEMCAL)
3982  {
3983  // Same sector
3984  Int_t j=0;
3985  for(Int_t i = 0; i < fNModules/2; i++)
3986  {
3987  j=2*i;
3988  if((module1==j && module2==j+1) || (module1==j+1 && module2==j)) fhMiSameSectorEMCALMod[i]->Fill(pt, m, GetEventWeight()) ;
3989  }
3990 
3991  // Same side
3992  for(Int_t i = 0; i < fNModules-2; i++)
3993  {
3994  if((module1==i && module2==i+2) || (module1==i+2 && module2==i)) fhMiSameSideEMCALMod[i]->Fill(pt, m, GetEventWeight());
3995  }
3996  } // EMCAL
3997  else
3998  { // PHOS
3999  if((module1==0 && module2==1) || (module1==1 && module2==0)) fhMiDiffPHOSMod[0]->Fill(pt, m, GetEventWeight()) ;
4000  if((module1==0 && module2==2) || (module1==2 && module2==0)) fhMiDiffPHOSMod[1]->Fill(pt, m, GetEventWeight()) ;
4001  if((module1==1 && module2==2) || (module1==2 && module2==1)) fhMiDiffPHOSMod[2]->Fill(pt, m, GetEventWeight()) ;
4002  } // PHOS
4003  }
4004  else
4005  {
4006  Float_t phi1 = GetPhi(fPhotonMom1.Phi());
4007  Float_t phi2 = GetPhi(fPhotonMom2.Phi());
4008  Bool_t etaside = 0;
4009  if( (p1->GetDetectorTag()==kEMCAL && fPhotonMom1.Eta() < 0)
4010  || (p2->GetDetectorTag()==kEMCAL && fPhotonMom2.Eta() < 0)) etaside = 1;
4011 
4012  if ( phi1 > DegToRad(260) && phi2 > DegToRad(260) && phi1 < DegToRad(280) && phi2 < DegToRad(280)) fhMiSameSectorDCALPHOSMod[0+etaside]->Fill(pt, m, GetEventWeight());
4013  else if ( phi1 > DegToRad(280) && phi2 > DegToRad(280) && phi1 < DegToRad(300) && phi2 < DegToRad(300)) fhMiSameSectorDCALPHOSMod[2+etaside]->Fill(pt, m, GetEventWeight());
4014  else if ( phi1 > DegToRad(300) && phi2 > DegToRad(300) && phi1 < DegToRad(320) && phi2 < DegToRad(320)) fhMiSameSectorDCALPHOSMod[4+etaside]->Fill(pt, m, GetEventWeight());
4015  else if ( (phi1 > DegToRad(260) && phi2 > DegToRad(280) && phi1 < DegToRad(280) && phi2 < DegToRad(300))
4016  || (phi1 > DegToRad(280) && phi2 > DegToRad(260) && phi1 < DegToRad(300) && phi2 < DegToRad(280))) fhMiDiffSectorDCALPHOSMod[0+etaside]->Fill(pt, m, GetEventWeight());
4017  else if ( (phi1 > DegToRad(280) && phi2 > DegToRad(300) && phi1 < DegToRad(300) && phi2 < DegToRad(320))
4018  || (phi1 > DegToRad(300) && phi2 > DegToRad(280) && phi1 < DegToRad(320) && phi2 < DegToRad(300))) fhMiDiffSectorDCALPHOSMod[2+etaside]->Fill(pt, m, GetEventWeight());
4019  else if ( (phi1 > DegToRad(260) && phi2 > DegToRad(300) && phi1 < DegToRad(280) && phi2 < DegToRad(320))
4020  || (phi1 > DegToRad(300) && phi2 > DegToRad(260) && phi1 < DegToRad(320) && phi2 < DegToRad(280))) fhMiDiffSectorDCALPHOSMod[4+etaside]->Fill(pt, m, GetEventWeight());
4021  else fhMiDiffSectorDCALPHOSMod[6+etaside]->Fill(pt, m, GetEventWeight());
4022  }
4023  }
4024 
4025  Bool_t ok = kTRUE;
4026  if(fSameSM)
4027  {
4029  {
4030  if(module1!=module2) ok=kFALSE;
4031  }
4032  else // PHOS and DCal in same sector
4033  {
4034  Float_t phi1 = GetPhi(fPhotonMom1.Phi());
4035  Float_t phi2 = GetPhi(fPhotonMom2.Phi());
4036  ok=kFALSE;
4037  if ( phi1 > DegToRad(260) && phi2 > DegToRad(260) && phi1 < DegToRad(280) && phi2 < DegToRad(280)) ok = kTRUE;
4038  else if ( phi1 > DegToRad(280) && phi2 > DegToRad(280) && phi1 < DegToRad(300) && phi2 < DegToRad(300)) ok = kTRUE;
4039  else if ( phi1 > DegToRad(300) && phi2 > DegToRad(300) && phi1 < DegToRad(320) && phi2 < DegToRad(320)) ok = kTRUE;
4040  }
4041  } // Pair only in same SM
4042 
4043  if(!ok) continue ;
4044 
4045  //
4046  // Do the final histograms with the selected clusters
4047  //
4048 
4049  // Check if one of the clusters comes from a conversion
4050  if(fCheckConversion)
4051  {
4052  if (p1->IsTagged() && p2->IsTagged()) fhMiConv2->Fill(pt, m, GetEventWeight());
4053  else if(p1->IsTagged() || p2->IsTagged()) fhMiConv ->Fill(pt, m, GetEventWeight());
4054  }
4055 
4056  //
4057  // Main invariant mass histograms
4058  // Fill histograms for different bad channel distance, centrality, assymmetry cut and pid bit
4059  //
4060  for(Int_t ipid=0; ipid<fNPIDBits; ipid++)
4061  {
4062  if((p1->IsPIDOK(ipid,AliCaloPID::kPhoton)) && (p2->IsPIDOK(ipid,AliCaloPID::kPhoton)))
4063  {
4064  for(Int_t iasym=0; iasym < fNAsymCuts; iasym++)
4065  {
4066  if(a < fAsymCuts[iasym])
4067  {
4068  Int_t index = ((curCentrBin*fNPIDBits)+ipid)*fNAsymCuts + iasym;
4069 
4070  if(index < 0 || index >= ncentr*fNPIDBits*fNAsymCuts) continue ;
4071 
4072  fhMi1[index]->Fill(pt, m, GetEventWeight()) ;
4073  if(fMakeInvPtPlots)fhMiInvPt1[index]->Fill(pt, m, 1./pt * GetEventWeight()) ;
4074 
4075  if(fFillBadDistHisto)
4076  {
4077  if(p1->DistToBad()>0 && p2->DistToBad()>0)
4078  {
4079  fhMi2[index]->Fill(pt, m, GetEventWeight()) ;
4080  if(fMakeInvPtPlots)fhMiInvPt2[index]->Fill(pt, m, 1./pt * GetEventWeight()) ;
4081 
4082  if(p1->DistToBad()>1 && p2->DistToBad()>1)
4083  {
4084  fhMi3[index]->Fill(pt, m, GetEventWeight()) ;
4085  if(fMakeInvPtPlots)fhMiInvPt3[index]->Fill(pt, m, 1./pt * GetEventWeight()) ;
4086  }
4087  }
4088  }// Fill bad dist histo
4089 
4090  }//Asymmetry cut
4091  }// Asymmetry loop
4092  }//PID cut
4093  }// PID loop
4094 
4095  //-----------------------
4096  // Multi cuts analysis
4097  //-----------------------
4098  Int_t ncell1 = p1->GetNCells();
4099  Int_t ncell2 = p1->GetNCells();
4100 
4101  if(fMultiCutAna)
4102  {
4103  // Several pt,ncell and asymmetry cuts
4104  for(Int_t ipt=0; ipt<fNPtCuts; ipt++)
4105  {
4106  for(Int_t icell=0; icell<fNCellNCuts; icell++)
4107  {
4108  for(Int_t iasym=0; iasym<fNAsymCuts; iasym++)
4109  {
4110  Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
4111 
4112  if(p1->Pt() > fPtCuts[ipt] && p2->Pt() > fPtCuts[ipt] &&
4113  p1->Pt() < fPtCutsMax[ipt] && p2->Pt() < fPtCutsMax[ipt] &&
4114  a < fAsymCuts[iasym] &&
4115  ncell1 >= fCellNCuts[icell] && ncell2 >= fCellNCuts[icell]
4116  )
4117  {
4118  //printf("MI ipt %d, iasym%d, icell %d, index %d \n",ipt, iasym, icell, index);
4119  //printf("\t %p, %p\n",fhMiPtNCellAsymCuts[index],fhMiPtNCellAsymCutsOpAngle[index]);
4120 
4121  fhMiPtNCellAsymCuts[index]->Fill(pt, m, GetEventWeight()) ;
4122  if(fFillAngleHisto) fhMiPtNCellAsymCutsOpAngle[index]->Fill(pt, angle, GetEventWeight()) ;
4123 
4124  //printf("ipt %d, icell%d, iasym %d, name %s\n",ipt, icell, iasym, fhRePtNCellAsymCuts[((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym]->GetName());
4125  }
4126  }// pid bit cut loop
4127  }// icell loop
4128  }// pt cut loop
4129  } // Multi cut ana
4130 
4131  //
4132  // Fill histograms with opening angle
4133  if(fFillAngleHisto)
4134  {
4135  fhMixedOpeningAngle ->Fill(pt, angle, GetEventWeight());
4136  fhMixedCosOpeningAngle->Fill(pt, TMath::Cos(angle), GetEventWeight());
4137  }
4138 
4139  //
4140  // Fill histograms for different opening angle bins
4142  {
4143  Int_t angleBin = -1;
4144  for(Int_t ibin = 0; ibin < fNAngleCutBins; ibin++)
4145  {
4146  if(angle >= fAngleCutBinsArray[ibin] &&
4147  angle < fAngleCutBinsArray[ibin+1]) angleBin = ibin;
4148  }
4149 
4150  if( angleBin >= 0 && angleBin < fNAngleCutBins)
4151  {
4152  Float_t e1 = fPhotonMom1.E();
4153  Float_t e2 = fPhotonMom2.E();
4154 
4155  Float_t t1 = p1->GetTime();
4156  Float_t t2 = p2->GetTime();
4157 
4158  Int_t nc1 = ncell1;
4159  Int_t nc2 = ncell2;
4160 
4161  Float_t eta1 = fPhotonMom1.Eta();
4162  Float_t eta2 = fPhotonMom2.Eta();
4163 
4164  Float_t phi1 = GetPhi(fPhotonMom1.Phi());
4165  Float_t phi2 = GetPhi(fPhotonMom2.Phi());
4166 
4167  Int_t mod1 = module1;
4168  Int_t mod2 = module2;
4169 
4170  // // Recover original cluster
4171  // Int_t iclus1 = -1, iclus2 = -1 ;
4172  // AliVCluster * cluster1 = FindCluster(GetEMCALClusters(),p1->GetCaloLabel(0),iclus1);
4173  // AliVCluster * cluster2 = FindCluster(GetEMCALClusters(),p2->GetCaloLabel(0),iclus2);
4174  //
4175  // Float_t maxCellFraction1 = 0, maxCellFraction2 = 0;
4176  // Int_t absIdMax1 = GetCaloUtils()->GetMaxEnergyCell(GetEMCALCells(),cluster1,maxCellFraction1);
4177  // Int_t absIdMax2 = GetCaloUtils()->GetMaxEnergyCell(GetEMCALCells(),cluster2,maxCellFraction2);
4178 
4179  if(e2 > e1)
4180  {
4181  e1 = fPhotonMom2.E();
4182  e2 = fPhotonMom1.E();
4183 
4184  t1 = p2->GetTime();
4185  t2 = p1->GetTime();
4186 
4187  nc1 = ncell2;
4188  nc2 = ncell1;
4189 
4190  eta1 = fPhotonMom2.Eta();
4191  eta2 = fPhotonMom1.Eta();
4192 
4193  phi1 = GetPhi(fPhotonMom2.Phi());
4194  phi2 = GetPhi(fPhotonMom1.Phi());
4195 
4196  mod1 = module2;
4197  mod2 = module1;
4198 
4199  // Int_t tmp = absIdMax2;
4200  // absIdMax2 = absIdMax1;
4201  // absIdMax1 = tmp;
4202  }
4203 
4204  fhMiOpAngleBinMinClusterEPerSM[angleBin]->Fill(e2,mod2,GetEventWeight()) ;
4205  fhMiOpAngleBinMaxClusterEPerSM[angleBin]->Fill(e1,mod1,GetEventWeight()) ;
4206 
4207  fhMiOpAngleBinMinClusterTimePerSM[angleBin]->Fill(t2,mod2,GetEventWeight()) ;
4208  fhMiOpAngleBinMaxClusterTimePerSM[angleBin]->Fill(t1,mod1,GetEventWeight()) ;
4209 
4210  fhMiOpAngleBinMinClusterNCellPerSM[angleBin]->Fill(nc2,mod2,GetEventWeight()) ;
4211  fhMiOpAngleBinMaxClusterNCellPerSM[angleBin]->Fill(nc1,mod1,GetEventWeight()) ;
4212 
4213  fhMiOpAngleBinPairClusterMass[angleBin]->Fill(pt,m,GetEventWeight()) ;
4214  if(mod2 == mod1) fhMiOpAngleBinPairClusterMassPerSM[angleBin]->Fill(m,mod1,GetEventWeight()) ;
4215 
4216  if(e1 > 0.01) fhMiOpAngleBinPairClusterRatioPerSM[angleBin]->Fill(e2/e1,mod1,GetEventWeight()) ;
4217 
4218  fhMiOpAngleBinMinClusterEtaPhi[angleBin]->Fill(eta2,phi2,GetEventWeight()) ;
4219  fhMiOpAngleBinMaxClusterEtaPhi[angleBin]->Fill(eta1,phi1,GetEventWeight()) ;
4220 
4221  // Int_t icol1 = -1, icol2 = -1, icolAbs1 = -1, icolAbs2 = -1;
4222  // Int_t irow1 = -1, irow2 = -1, irowAbs1 = -1, irowAbs2 = -1;
4223  // Int_t iRCU1 = -1, iRCU2 = -1;
4224  // GetModuleNumberCellIndexesAbsCaloMap(absIdMax1,GetCalorimeter(), icol1, irow1, iRCU1, icolAbs1, irowAbs1);
4225  // GetModuleNumberCellIndexesAbsCaloMap(absIdMax2,GetCalorimeter(), icol2, irow2, iRCU1, icolAbs2, irowAbs2);
4226  //
4227  // fhMiOpAngleBinPairClusterAbsIdMaxCell[angleBin]->Fill(absIdMax1,absIdMax2,GetEventWeight());
4228  //
4229  // fhMiColRowClusterMinOpAngleBin[angleBin]->Fill(icolAbs2,irowAbs2,GetEventWeight()) ;
4230  // fhMiOpAngleBinMaxClusterColRow[angleBin]->Fill(icolAbs1,irowAbs1,GetEventWeight()) ;
4231  }
4232  }
4233 
4234  // Fill histograms with pair assymmetry
4236  {
4237  fhMiPtAsym->Fill(pt, a, GetEventWeight());
4238  if ( m > fPi0MassWindow[0] && m < fPi0MassWindow[1] ) fhMiPtAsymPi0->Fill(pt, a, GetEventWeight());
4239  if ( m > fEtaMassWindow[0] && m < fEtaMassWindow[1] ) fhMiPtAsymEta->Fill(pt, a, GetEventWeight());
4240  }
4241 
4242  // Check cell time content in cluster
4244  {
4245  if ( p1->GetFiducialArea() == 0 && p2->GetFiducialArea() == 0 )
4247 
4248  else if ( p1->GetFiducialArea() != 0 && p2->GetFiducialArea() != 0 )
4250  }
4251 
4252  }// second cluster loop
4253  }//first cluster loop
4254  }//loop on mixed events
4255 
4256  //--------------------------------------------------------
4257  // Add the current event to the list of events for mixing
4258  //--------------------------------------------------------
4259 
4260  //TClonesArray *currentEvent = new TClonesArray(*GetInputAODBranch());
4261  TClonesArray *currentEvent = new TClonesArray(*secondLoopInputData);
4262 
4263  // Add current event to buffer and Remove redundant events
4264  if( currentEvent->GetEntriesFast() > 0 )
4265  {
4266  evMixList->AddFirst(currentEvent) ;
4267  currentEvent=0 ; //Now list of particles belongs to buffer and it will be deleted with buffer
4268  if( evMixList->GetSize() >= GetNMaxEvMix() )
4269  {
4270  TClonesArray * tmp = (TClonesArray*) (evMixList->Last()) ;
4271  evMixList->RemoveLast() ;
4272  delete tmp ;
4273  }
4274  }
4275  else
4276  { // empty event
4277  delete currentEvent ;
4278  currentEvent=0 ;
4279  }
4280  }// DoOwnMix
4281 
4282  AliDebug(1,"End fill histograms");
4283 }
4284 
4285 //________________________________________________________________________
4289 //________________________________________________________________________
4290 Int_t AliAnaPi0::GetEventIndex(AliAODPWG4Particle * part, Double_t * vert)
4291 {
4292  Int_t evtIndex = -1 ;
4293 
4294  if(GetReader()->GetDataType()!=AliCaloTrackReader::kMC)
4295  {
4296  if (GetMixedEvent())
4297  {
4298  evtIndex = GetMixedEvent()->EventIndexForCaloCluster(part->GetCaloLabel(0)) ;
4299  GetVertex(vert,evtIndex);
4300 
4301  if(TMath::Abs(vert[2])> GetZvertexCut())
4302  evtIndex = -2 ; //Event can not be used (vertex, centrality,... cuts not fulfilled)
4303  }
4304  else
4305  {
4306  // Single event
4307  GetVertex(vert);
4308 
4309  if(TMath::Abs(vert[2])> GetZvertexCut())
4310  evtIndex = -1 ; //Event can not be used (vertex, centrality,... cuts not fulfilled)
4311  else
4312  evtIndex = 0 ;
4313  }
4314  } // No MC reader
4315  else
4316  {
4317  evtIndex = 0;
4318  vert[0] = 0. ;
4319  vert[1] = 0. ;
4320  vert[2] = 0. ;
4321  }
4322 
4323  return evtIndex ;
4324 }
4325 
TH2F * fhPrimEtaY
! Rapidity distribution of primary particles vs pT
Definition: AliAnaPi0.h:415
Float_t GetHistoPtMax() const
TH2F * fhPairGeneratorsBkgMassMCPi0[10][10]
! Mass for a pair of clusters with depending bkg type, pi0 true pairs
Definition: AliAnaPi0.h:554
TH2F ** fhRe3
[GetNCentrBin()*fNPIDBits*fNAsymCuts]
Definition: AliAnaPi0.h:302
TH2F * fhPrimEtaAccPtCentrality
! primary eta with accepted daughters reconstructed centrality vs pT
Definition: AliAnaPi0.h:428
TH2F * fhMiOpAngleBinMinClusterEPerSM[10]
! energy of lowest energy cluster in pair, depending on opening angle cut, y axis is SM number...
Definition: AliAnaPi0.h:539
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:563
TH2F * fhReOpAngleBinMaxClusterEPerSM[10]
! energy of highest energy cluster in pair, depending on opening angle cut, y axis is SM number ...
Definition: AliAnaPi0.h:522
Int_t pdg
TH2F * fhPrimPi0Y
! Rapidity distribution of primary particles vs pT
Definition: AliAnaPi0.h:391
TLorentzVector fPhotonMom1Boost
! Photon cluster momentum, temporary array
Definition: AliAnaPi0.h:234
TH2F * fhMiOpAngleBinPairClusterRatioPerSM[10]
! lowest/highest energy cluster in pair, depending on opening angle cut, y axis is SM number...
Definition: AliAnaPi0.h:545
Int_t fNModules
[GetNCentrBin()*GetNZvertBin()*GetNRPBin()]
Definition: AliAnaPi0.h:186
TH2F * fhPrimEtaYetaYcut
! PseudoRapidity distribution of primary particles vs pT, Y<1
Definition: AliAnaPi0.h:418
TH1F * fhPrimPi0AccPt
! Spectrum of primary with accepted daughters
Definition: AliAnaPi0.h:388
Float_t GetHistoPtMin() const
TH2F * fhPrimEtaAccPhi
! Azimutal distribution of primary with accepted daughters vs pT
Definition: AliAnaPi0.h:421
double Double_t
Definition: External.C:58
TH2F * fhReMCFromMixConversion
! Invariant mass of 2 clusters one from conversion and the other not
Definition: AliAnaPi0.h:501
Int_t fPIDBits[10]
Array with different PID bits.
Definition: AliAnaPi0.h:208
TH2F * fhPairGeneratorsBkgEPrimRecoDiffMCPi0[10][10]
! pT reco - pT primary for a pair of clusters with depending bkg type, pi0 true pairs ...
Definition: AliAnaPi0.h:556
TH2F * fhMiPtAsymPi0
! Mix two-photon pt vs asymmetry, close to pi0 mass
Definition: AliAnaPi0.h:360
TH2F * fhMCPi0PtRecOpenAngleMassCut
! Real pi0 pairs, reco pT vs reco opening angle, inside a mass window
Definition: AliAnaPi0.h:483
virtual void AddToHistogramsName(TString add)
virtual AliVCaloCells * GetEMCALCells() const
TH2F * fhReConv
[8]
Definition: AliAnaPi0.h:281
TH2F * fhReOpAngleBinPairClusterMassMCTrueEta[10]
! cluster pair mass vs pT, depending on opening angle cut, true eta decay pairs from MC ...
Definition: AliAnaPi0.h:533
Definition: External.C:236
TH2F * fhMCPi0PtTruePtRecDif
! Real pi0 pairs, reco pT vs pT difference generated - reco
Definition: AliAnaPi0.h:474
Bool_t fFillArmenterosThetaStar
Fill armenteros histograms.
Definition: AliAnaPi0.h:223
Int_t fNPIDBits
Number of possible PID bit combinations.
Definition: AliAnaPi0.h:207
TLorentzVector fPhotonMom2
! Photon cluster momentum, temporary array
Definition: AliAnaPi0.h:235
TH2F * fhReOpAngleBinMaxClusterColRow[10]
! Column and row location of main cell of highest energy cluster in pair, depending on opening angle ...
Definition: AliAnaPi0.h:520
const char * title
Definition: MakeQAPdf.C:26
TH2F ** fhMiDiffSectorDCALPHOSMod
[6]
Definition: AliAnaPi0.h:277
Int_t fCellNCuts[10]
Array with different cell number cluster cuts.
Definition: AliAnaPi0.h:206
TH2F * fhMCOrgDeltaEta[13]
! Delta Eta vs pt of real pairs, check common origin of pair
Definition: AliAnaPi0.h:444
TH2F * fhReSS[3]
[fNPtCuts*fNAsymCuts*fNCellNCuts]
Definition: AliAnaPi0.h:352
TH2F * fhReOpAngleBinMaxClusterNCellPerSM[10]
! N cells of highest energy cluster in pair, depending on opening angle cut, y axis is SM number ...
Definition: AliAnaPi0.h:526
TH2F * fhPrimEtaCosOpeningAngle
! Cosinus of opening angle of pair version pair energy, eta primaries
Definition: AliAnaPi0.h:425
Int_t GetHistoEDiffBins() const
TH2F * fhPrimPi0PtStatus
! Spectrum of generated pi0 vs pi0 status
Definition: AliAnaPi0.h:436
Bool_t fCheckConversion
Fill histograms with tagged photons as conversion.
Definition: AliAnaPi0.h:217
virtual Int_t GetModuleNumberCellIndexesAbsCaloMap(Int_t absId, Int_t calo, Int_t &icol, Int_t &irow, Int_t &iRCU, Int_t &icolAbs, Int_t &irowAbs) const
TH2F * fhPtBinClusterEtaPhi[10]
! Eta-Phi location of cluster in different energy bins.
Definition: AliAnaPi0.h:550
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:385
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:189
TH2F * fhMCEtaPtTruePtRecRat
! Real pi0 pairs, reco pT vs pT ratio reco / generated
Definition: AliAnaPi0.h:477
TH2F ** fhMiSameSectorEMCALMod
[fNModules-2]
Definition: AliAnaPi0.h:268
TH2F * fhMixedOpeningAnglePerSM[20]
! Opening angle of pair versus pair energy, per SM
Definition: AliAnaPi0.h:380
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:367
TH2F * fhMiSecondaryCellOutTimeWindow
! Combine clusters when at least one significant cells in cluster has t > 50 ns, different events ...
Definition: AliAnaPi0.h:514
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:230
Float_t DegToRad(Float_t deg) const
TH2F * fhPrimPi0PtCentrality
! primary pi0 reconstructed centrality vs pT
Definition: AliAnaPi0.h:402
TH2F * fhPrimEtaPtCentrality
! primary eta reconstructed centrality vs pT
Definition: AliAnaPi0.h:426
TH2F * fhPairGeneratorsBkgEPrimRecoRatioMCPi0[10][10]
! pT reco / pT primary for a pair of clusters with depending bkg type, pi0 true pairs ...
Definition: AliAnaPi0.h:555
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:473
TH2F * fhMCEtaPtRecOpenAngle
! Real pi0 pairs, reco pT vs reco opening angle
Definition: AliAnaPi0.h:479
TH2F * fhReOpAngleBinPairClusterMassPerSM[10]
! cluster pair mass, depending on opening angle cut, y axis is SM number
Definition: AliAnaPi0.h:529
TH2F * fhReMCFromConversion
! Invariant mass of 2 clusters originated in conversions
Definition: AliAnaPi0.h:499
TH2F ** fhReSameSideEMCALMod
[fNModules]
Definition: AliAnaPi0.h:247
TH1F * fhPrimPi0AccPtOpAngCuts[10]
! Spectrum of primary with accepted daughters, different opening angles
Definition: AliAnaPi0.h:390
TH1F * fhPrimPi0AccPtPhotonCuts
! Spectrum of primary with accepted daughters, photon pt or angle cuts
Definition: AliAnaPi0.h:389
TH2F * fhPrimPi0OpeningAngle
! Opening angle of pair versus pair energy, primaries
Definition: AliAnaPi0.h:398
virtual AliNeutralMesonSelection * GetNeutralMesonSelection()
TH2F * fhRealOpeningAngle
! Opening angle of pair versus pair energy
Definition: AliAnaPi0.h:374
TH2F * fhMCOrgMass[13]
! Mass vs pt of real pairs, check common origin of pair
Definition: AliAnaPi0.h:442
TH2F * fhMiOpAngleBinMinClusterEtaPhi[10]
! Eta-Phi location of lowest energy cluster in pair, depending on opening angle cut, mixed event
Definition: AliAnaPi0.h:535
TH2F * fhMiOpAngleBinPairClusterMass[10]
! cluster pair mass vs pT, depending on opening angle cut, mixed event
Definition: AliAnaPi0.h:546
Int_t GetHistoMassBins() const
Int_t GetHistoPhiBins() const
TH2F * fhRePtAsymPi0
! REAL two-photon pt vs asymmetry, close to pi0 mass
Definition: AliAnaPi0.h:357
TLorentzVector fPhotonMom1
! Photon cluster momentum, temporary array
Definition: AliAnaPi0.h:233
TH2F * fhRealCosOpeningAngle
! Cosinus of opening angle of pair version pair energy
Definition: AliAnaPi0.h:375
TH2F ** fhRe1
REAL two-photon invariant mass distribution for different centralities and Asymmetry.
Definition: AliAnaPi0.h:287
Float_t GetHistoMassMin() const
TH2F * fhPrimPi0ProdVertex
! Spectrum of primary pi0 vs production vertex
Definition: AliAnaPi0.h:496
TH2F ** fhReDiffPHOSMod
[fNModules/2]
Definition: AliAnaPi0.h:253
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:282
TH2F * fhReOpAngleBinPairClusterMass[10]
! cluster pair mass vs pT, depending on opening angle cut
Definition: AliAnaPi0.h:528
Float_t fAsymCuts[10]
Array with different assymetry cuts.
Definition: AliAnaPi0.h:204
TString GetCocktailGenNameToCheck(Int_t i) const
TH2F * fhMCOrgAsym[13]
! Asymmetry vs pt of real pairs, check common origin of pair
Definition: AliAnaPi0.h:443
TH2F ** fhRePtNCellAsymCuts
[GetNCentrBin()*fNPIDBits*fNAsymCuts]
Definition: AliAnaPi0.h:335
TH2F * fhPrimPi0PtOrigin
! Spectrum of generated pi0 vs mother
Definition: AliAnaPi0.h:433
Float_t GetHistoDiffTimeMin() const
Float_t fPtCuts[11]
Array with different pt cuts, minimum.
Definition: AliAnaPi0.h:201
TH2F * fhRealOpeningAnglePerSM[20]
! Opening angle of pair versus pair energy, per SM
Definition: AliAnaPi0.h:379
TH2F * fhMiOpAngleBinMaxClusterNCellPerSM[10]
! N cells of highest energy cluster in pair, depending on opening angle cut, y axis is SM number...
Definition: AliAnaPi0.h:544
Bool_t IsMCParticleInCalorimeterAcceptance(Int_t calo, TParticle *particle)
TH2F ** fhMCPi0MassPtRec
Real pi0 pairs, reconstructed mass vs reconstructed pt of original pair.
Definition: AliAnaPi0.h:450
Bool_t fFillOpAngleCutHisto
Fill histograms depending on opening angle of pair.
Definition: AliAnaPi0.h:226
TH2F * fhMiPtAsym
! Mix two-photon pt vs asymmetry
Definition: AliAnaPi0.h:359
TH2F * fhPrimPi0AccPtCentrality
! primary pi0 with accepted daughters reconstructed centrality vs pT
Definition: AliAnaPi0.h:404
Float_t GetHistoPhiMin() const
Int_t fNAsymCuts
Number of assymmetry cuts.
Definition: AliAnaPi0.h:203
Float_t GetHistoDiffTimeMax() const
TH2F * fhMixedCosOpeningAngle
! Cosinus of opening angle of pair version pair energy
Definition: AliAnaPi0.h:377
TH2F ** fhReInvPt2
[GetNCentrBin()*fNPIDBits*fNAsymCuts]
Definition: AliAnaPi0.h:318
TH2F ** fhMiDiffPHOSMod
[fNModules/2]
Definition: AliAnaPi0.h:271
TH2F ** fhMi2
[GetNCentrBin()*fNPIDBits*fNAsymCuts]
Definition: AliAnaPi0.h:298
Float_t GetHistoOpeningAngleMax() const
TH2F * fhMCPi0PtTruePtRecRatMassCut
! Real pi0 pairs, reco pT vs pT ratio reco / generated, inside a mass window
Definition: AliAnaPi0.h:481
TList * GetCreateOutputObjects()
Definition: AliAnaPi0.cxx:333
TH2F ** fhReInvPt3
[GetNCentrBin()*fNPIDBits*fNAsymCuts]
Definition: AliAnaPi0.h:326
Bool_t fFillOriginHisto
Fill histograms depending on their origin.
Definition: AliAnaPi0.h:222
Float_t GetHistoEDiffMin() const
TH2F * fhPrimPi0Phi
! Azimutal distribution of primary particles vs pT
Definition: AliAnaPi0.h:396
TH2F * fhPrimEtaOpeningAnglePhotonCuts
! Opening angle of pair versus pair energy, eta primaries, photon pT cuts
Definition: AliAnaPi0.h:423
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:562
TLorentzVector fMCPrimMesonMom
! Pi0/Eta MC primary momentum, temporary array
Definition: AliAnaPi0.h:236
TH2F ** fhReSameSectorDCALPHOSMod
[fNModules]
Definition: AliAnaPi0.h:256
TH2F ** fhMiPtNCellAsymCutsOpAngle
[fNPtCuts*fNAsymCuts*fNCellNCuts]
Definition: AliAnaPi0.h:347
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:540
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:220
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:558
Bool_t fFillBadDistHisto
Do plots for different distances to bad channels.
Definition: AliAnaPi0.h:218
Float_t fAngleCutBinsArray[11]
Array with angle cut bins.
Definition: AliAnaPi0.h:211
TH2F * fhPairGeneratorsBkgMass[10][10]
! Mass for a pair of clusters depending bkg type
Definition: AliAnaPi0.h:553
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:468
Bool_t fMultiCutAna
Do analysis with several or fixed cut.
Definition: AliAnaPi0.h:197
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:244
TH1F * fhCentralityNoPair
! Histogram with centrality bins with no pair
Definition: AliAnaPi0.h:368
Bool_t fFillAsymmetryHisto
Fill histograms with asymmetry vs pt.
Definition: AliAnaPi0.h:221
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:523
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:525
Int_t GetEventIndex(AliAODPWG4Particle *part, Double_t *vert)
Definition: AliAnaPi0.cxx:4290
float Float_t
Definition: External.C:68
void FillAcceptanceHistograms()
Fill acceptance histograms if MC data is available.
Definition: AliAnaPi0.cxx:2261
TH2F ** fhReSameSectorEMCALMod
[fNModules-2]
Definition: AliAnaPi0.h:250
TH2F ** fhMCEtaMassPtTrue
[fNPtCuts*fNAsymCuts*fNCellNCuts]
Definition: AliAnaPi0.h:465
TH2F ** fhReInvPt1
[GetNCentrBin()*fNPIDBits*fNAsymCuts]
Definition: AliAnaPi0.h:311
Int_t fNCellNCuts
Number of cuts with number of cells in cluster.
Definition: AliAnaPi0.h:205
TH2F ** fhMiInvPt1
[GetNCentrBin()*fNPIDBits*fNAsymCuts]
Definition: AliAnaPi0.h:314
Int_t GetHistoAsymmetryBins() const
const Double_t ptmax
TH2F * fhPrimEtaAccYeta
! PseudoRapidity distribution of primary with accepted daughters vs pT
Definition: AliAnaPi0.h:419
Float_t fAngleMaxCut
Select pairs with opening angle smaller than a threshold.
Definition: AliAnaPi0.h:191
TH2F * fhReOpAngleBinPairClusterMassMCTruePi0[10]
! cluster pair mass vs pT, depending on opening angle cut, true pi0 decay pairs from MC ...
Definition: AliAnaPi0.h:532
TH1I * fhEventMixBin
! Number of mixed pairs in a particular bin (cen,vz,rp)
Definition: AliAnaPi0.h:366
void FillArmenterosThetaStar(Int_t pdg)
Fill armenteros plots.
Definition: AliAnaPi0.cxx:2713
virtual AliEMCALGeometry * GetEMCALGeometry() const
TH2F * fhPrimPi0CosOpeningAngle
! Cosinus of opening angle of pair version pair energy, pi0 primaries
Definition: AliAnaPi0.h:401
TH2F * fhMCPi0PtOrigin
! Mass of reconstructed pi0 pairs in calorimeter vs mother origin ID.
Definition: AliAnaPi0.h:489
TH2F ** fhMiSameSideEMCALMod
[fNModules]
Definition: AliAnaPi0.h:265
TH1F * fhPrimEtaAccPtOpAngCuts[10]
! Spectrum of primary with accepted daughters, different opening angles
Definition: AliAnaPi0.h:414
Bool_t fFillOnlyMCAcceptanceHisto
Do analysis only of MC kinematics input.
Definition: AliAnaPi0.h:224