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