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