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