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