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