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