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