AliPhysics  0e0bd91 (0e0bd91)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros
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),
58 fNPtCuts(0), fNAsymCuts(0), fNCellNCuts(0), fNPIDBits(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 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 fhRePIDBits(0x0), fhRePtMult(0x0), fhReSS(),
82 fhRePtAsym(0x0), fhRePtAsymPi0(0x0), fhRePtAsymEta(0x0),
83 fhEventBin(0), fhEventMixBin(0),
84 fhCentrality(0x0), fhCentralityNoPair(0x0),
85 fhEventPlaneResolution(0x0),
86 fhRealOpeningAngle(0x0), fhRealCosOpeningAngle(0x0), fhMixedOpeningAngle(0x0), fhMixedCosOpeningAngle(0x0),
87 
88 // MC histograms
89 fhPrimPi0E(0x0), fhPrimPi0Pt(0x0),
90 fhPrimPi0AccE(0x0), fhPrimPi0AccPt(0x0), fhPrimPi0AccPtPhotonCuts(0x0),
91 fhPrimPi0Y(0x0), fhPrimPi0AccY(0x0),
92 fhPrimPi0Yeta(0x0), fhPrimPi0YetaYcut(0x0), fhPrimPi0AccYeta(0x0),
93 fhPrimPi0Phi(0x0), fhPrimPi0AccPhi(0x0),
94 fhPrimPi0OpeningAngle(0x0), fhPrimPi0OpeningAnglePhotonCuts(0x0),
95 fhPrimPi0OpeningAngleAsym(0x0),fhPrimPi0CosOpeningAngle(0x0),
96 fhPrimPi0PtCentrality(0), fhPrimPi0PtEventPlane(0),
97 fhPrimPi0AccPtCentrality(0), fhPrimPi0AccPtEventPlane(0),
98 fhPrimEtaE(0x0), fhPrimEtaPt(0x0),
99 fhPrimEtaAccE(0x0), fhPrimEtaAccPt(0x0), fhPrimEtaAccPtPhotonCuts(0x0),
100 fhPrimEtaY(0x0), fhPrimEtaAccY(0x0),
101 fhPrimEtaYeta(0x0), fhPrimEtaYetaYcut(0x0), fhPrimEtaAccYeta(0x0),
102 fhPrimEtaPhi(0x0), fhPrimEtaAccPhi(0x0),
103 fhPrimEtaOpeningAngle(0x0), fhPrimEtaOpeningAnglePhotonCuts(0x0),
104 fhPrimEtaOpeningAngleAsym(0x0),fhPrimEtaCosOpeningAngle(0x0),
105 fhPrimEtaPtCentrality(0), fhPrimEtaPtEventPlane(0),
106 fhPrimEtaAccPtCentrality(0), fhPrimEtaAccPtEventPlane(0),
107 fhPrimPi0PtOrigin(0x0), fhPrimEtaPtOrigin(0x0),
108 fhPrimNotResonancePi0PtOrigin(0x0), fhPrimPi0PtStatus(0x0),
109 fhMCOrgMass(), fhMCOrgAsym(), fhMCOrgDeltaEta(), fhMCOrgDeltaPhi(),
110 fhMCPi0MassPtRec(), fhMCPi0MassPtTrue(), fhMCPi0PtTruePtRec(),
111 fhMCEtaMassPtRec(), fhMCEtaMassPtTrue(), fhMCEtaPtTruePtRec(),
112 fhMCPi0PtOrigin(0x0), fhMCEtaPtOrigin(0x0),
113 fhMCNotResonancePi0PtOrigin(0),fhMCPi0PtStatus(0x0),
114 fhMCPi0ProdVertex(0), fhMCEtaProdVertex(0),
115 fhPrimPi0ProdVertex(0), fhPrimEtaProdVertex(0),
116 fhReMCFromConversion(0), fhReMCFromNotConversion(0), fhReMCFromMixConversion(0),
117 fhCosThStarPrimPi0(0), fhCosThStarPrimEta(0),
118 fhEPairDiffTime(0)
119 {
120  InitParameters();
121 
122  for(Int_t i = 0; i < 4; i++)
123  {
124  fhArmPrimEta[i] = 0;
125  fhArmPrimPi0[i] = 0;
126  }
127 }
128 
129 //_____________________
131 //_____________________
133 {
134  // Remove event containers
135 
136  if(DoOwnMix() && fEventsList)
137  {
138  for(Int_t ic=0; ic<GetNCentrBin(); ic++)
139  {
140  for(Int_t iz=0; iz<GetNZvertBin(); iz++)
141  {
142  for(Int_t irp=0; irp<GetNRPBin(); irp++)
143  {
144  Int_t bin = GetEventMixBin(ic,iz,irp);
145  fEventsList[bin]->Delete() ;
146  delete fEventsList[bin] ;
147  }
148  }
149  }
150  delete[] fEventsList;
151  }
152 }
153 
154 //______________________________
157 //______________________________
159 {
160  SetInputAODName("PWG4Particle");
161 
162  AddToHistogramsName("AnaPi0_");
163 
164  fUseAngleEDepCut = kFALSE;
165 
166  fUseAngleCut = kTRUE;
167  fAngleCut = 0.;
168  fAngleMaxCut = DegToRad(80.); // 80 degrees cut, avoid EMCal/DCal combinations
169 
170  fMultiCutAna = kFALSE;
171 
172  fNPtCuts = 3;
173  fPtCuts[0] = 0.; fPtCuts[1] = 0.3; fPtCuts[2] = 0.5;
174  for(Int_t i = fNPtCuts; i < 10; i++)fPtCuts[i] = 0.;
175 
176  fNAsymCuts = 2;
177  fAsymCuts[0] = 1.; fAsymCuts[1] = 0.7; //fAsymCuts[2] = 0.6; // fAsymCuts[3] = 0.1;
178  for(Int_t i = fNAsymCuts; i < 10; i++)fAsymCuts[i] = 0.;
179 
180  fNCellNCuts = 3;
181  fCellNCuts[0] = 0; fCellNCuts[1] = 1; fCellNCuts[2] = 2;
182  for(Int_t i = fNCellNCuts; i < 10; i++)fCellNCuts[i] = 0;
183 
184  fNPIDBits = 2;
185  fPIDBits[0] = 0; fPIDBits[1] = 2; // fPIDBits[2] = 4; fPIDBits[3] = 6;// check, no cut, dispersion, neutral, dispersion&&neutral
186  for(Int_t i = fNPIDBits; i < 10; i++)fPIDBits[i] = 0;
187 }
188 
189 //_______________________________________
191 //_______________________________________
193 {
194  TString parList ; //this will be list of parameters used for this analysis.
195  const Int_t buffersize = 255;
196  char onePar[buffersize] ;
197  snprintf(onePar,buffersize,"--- AliAnaPi0 ---:") ;
198  parList+=onePar ;
199  snprintf(onePar,buffersize,"Number of bins in Centrality: %d;",GetNCentrBin()) ;
200  parList+=onePar ;
201  snprintf(onePar,buffersize,"Number of bins in Z vert. pos: %d;",GetNZvertBin()) ;
202  parList+=onePar ;
203  snprintf(onePar,buffersize,"Number of bins in Reac. Plain: %d;",GetNRPBin()) ;
204  parList+=onePar ;
205  snprintf(onePar,buffersize,"Depth of event buffer: %d;",GetNMaxEvMix()) ;
206  parList+=onePar ;
207  snprintf(onePar,buffersize,"Select pairs with their angle: %d, edep %d, min angle %2.3f, max angle %2.3f;",fUseAngleCut, fUseAngleEDepCut,fAngleCut,fAngleMaxCut) ;
208  parList+=onePar ;
209  snprintf(onePar,buffersize," Asymmetry cuts: n = %d, asymmetry < ",fNAsymCuts) ;
210  for(Int_t i = 0; i < fNAsymCuts; i++) snprintf(onePar,buffersize,"%s %2.2f;",onePar,fAsymCuts[i]);
211  parList+=onePar ;
212  snprintf(onePar,buffersize," PID selection bits: n = %d, PID bit =",fNPIDBits) ;
213  for(Int_t i = 0; i < fNPIDBits; i++) snprintf(onePar,buffersize,"%s %d;",onePar,fPIDBits[i]);
214  parList+=onePar ;
215  snprintf(onePar,buffersize,"Cuts:") ;
216  parList+=onePar ;
217  snprintf(onePar,buffersize,"Z vertex position: -%f < z < %f;",GetZvertexCut(),GetZvertexCut()) ;
218  parList+=onePar ;
219  snprintf(onePar,buffersize,"Calorimeter: %s;",GetCalorimeterString().Data()) ;
220  parList+=onePar ;
221  snprintf(onePar,buffersize,"Number of modules: %d:",fNModules) ;
222  parList+=onePar ;
223  if(fMultiCutAna)
224  {
225  snprintf(onePar, buffersize," pT cuts: n = %d, pt > ",fNPtCuts) ;
226  for(Int_t i = 0; i < fNPtCuts; i++) snprintf(onePar,buffersize,"%s %2.2f;",onePar,fPtCuts[i]);
227  parList+=onePar ;
228  snprintf(onePar,buffersize, " N cell in cluster cuts: n = %d, nCell > ",fNCellNCuts) ;
229  for(Int_t i = 0; i < fNCellNCuts; i++) snprintf(onePar,buffersize,"%s %d;",onePar,fCellNCuts[i]);
230  parList+=onePar ;
231  }
232 
233  return new TObjString(parList) ;
234 }
235 
236 //_________________________________________
239 //_________________________________________
241 {
242  TList * outputContainer = new TList() ;
243  outputContainer->SetName(GetName());
244 
245  // Set sizes and ranges
246  const Int_t buffersize = 255;
247  char key[buffersize] ;
248  char title[buffersize] ;
249 
251  Int_t nphibins = GetHistogramRanges()->GetHistoPhiBins();
252  Int_t netabins = GetHistogramRanges()->GetHistoEtaBins();
253  Float_t ptmax = GetHistogramRanges()->GetHistoPtMax();
254  Float_t phimax = GetHistogramRanges()->GetHistoPhiMax();
256  Float_t ptmin = GetHistogramRanges()->GetHistoPtMin();
259 
260  Int_t nmassbins = GetHistogramRanges()->GetHistoMassBins();
261  Int_t nasymbins = GetHistogramRanges()->GetHistoAsymmetryBins();
262  Float_t massmax = GetHistogramRanges()->GetHistoMassMax();
263  Float_t asymmax = GetHistogramRanges()->GetHistoAsymmetryMax();
264  Float_t massmin = GetHistogramRanges()->GetHistoMassMin();
265  Float_t asymmin = GetHistogramRanges()->GetHistoAsymmetryMin();
269  Int_t tdbins = GetHistogramRanges()->GetHistoDiffTimeBins() ;
270  Float_t tdmax = GetHistogramRanges()->GetHistoDiffTimeMax();
271  Float_t tdmin = GetHistogramRanges()->GetHistoDiffTimeMin();
272 
273  // Start with pure MC kinematics histograms
274  // In case other tasks just need this info like AliAnaPi0EbE
275  if(IsDataMC())
276  {
277  // Pi0
278 
279  fhPrimPi0E = new TH1F("hPrimPi0E","Primary #pi^{0} E, |#it{Y}|<1",nptbins,ptmin,ptmax) ;
280  fhPrimPi0AccE = new TH1F("hPrimPi0AccE","Primary #pi^{0} #it{E} with both photons in acceptance",nptbins,ptmin,ptmax) ;
281  fhPrimPi0E ->SetXTitle("#it{E} (GeV)");
282  fhPrimPi0AccE->SetXTitle("#it{E} (GeV)");
283  outputContainer->Add(fhPrimPi0E) ;
284  outputContainer->Add(fhPrimPi0AccE) ;
285 
286  fhPrimPi0Pt = new TH1F("hPrimPi0Pt","Primary #pi^{0} #it{p}_{T} , |#it{Y}|<1",nptbins,ptmin,ptmax) ;
287  fhPrimPi0AccPt = new TH1F("hPrimPi0AccPt","Primary #pi^{0} #it{p}_{T} with both photons in acceptance",nptbins,ptmin,ptmax) ;
288  fhPrimPi0AccPtPhotonCuts = new TH1F("hPrimPi0AccPtPhotonCuts","Primary #pi^{0} #it{p}_{T} with both photons in acceptance",nptbins,ptmin,ptmax) ;
289  fhPrimPi0Pt ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
290  fhPrimPi0AccPt->SetXTitle("#it{p}_{T} (GeV/#it{c})");
291  fhPrimPi0AccPtPhotonCuts->SetXTitle("#it{p}_{T} (GeV/#it{c})");
292  outputContainer->Add(fhPrimPi0Pt) ;
293  outputContainer->Add(fhPrimPi0AccPt) ;
294  outputContainer->Add(fhPrimPi0AccPtPhotonCuts) ;
295 
296  Int_t netabinsopen = TMath::Nint(netabins*4/(etamax-etamin));
297  fhPrimPi0Y = new TH2F("hPrimPi0Rapidity","Rapidity of primary #pi^{0}",nptbins,ptmin,ptmax,netabinsopen,-2, 2) ;
298  fhPrimPi0Y ->SetYTitle("#it{Rapidity}");
299  fhPrimPi0Y ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
300  outputContainer->Add(fhPrimPi0Y) ;
301 
302  fhPrimPi0Yeta = new TH2F("hPrimPi0PseudoRapidity","PseudoRapidity of primary #pi^{0}",nptbins,ptmin,ptmax,netabinsopen,-2, 2) ;
303  fhPrimPi0Yeta ->SetYTitle("#eta");
304  fhPrimPi0Yeta ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
305  outputContainer->Add(fhPrimPi0Yeta) ;
306 
307  fhPrimPi0YetaYcut = new TH2F("hPrimPi0PseudoRapidityYcut","PseudoRapidity of primary #pi^{0}, |#it{Y}|<1",nptbins,ptmin,ptmax,netabinsopen,-2, 2) ;
308  fhPrimPi0YetaYcut ->SetYTitle("#eta");
309  fhPrimPi0YetaYcut ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
310  outputContainer->Add(fhPrimPi0YetaYcut) ;
311 
312  fhPrimPi0AccY = new TH2F("hPrimPi0AccRapidity","Rapidity of primary #pi^{0} with accepted daughters",nptbins,ptmin,ptmax,netabins,etamin,etamax) ;
313  fhPrimPi0AccY->SetYTitle("Rapidity");
314  fhPrimPi0AccY->SetXTitle("#it{p}_{T} (GeV/#it{c})");
315  outputContainer->Add(fhPrimPi0AccY) ;
316 
317  fhPrimPi0AccYeta = new TH2F("hPrimPi0AccPseudoRapidity","PseudoRapidity of primary #pi^{0} with accepted daughters",nptbins,ptmin,ptmax,netabins,etamin,etamax) ;
318  fhPrimPi0AccYeta ->SetYTitle("#eta");
319  fhPrimPi0AccYeta ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
320  outputContainer->Add(fhPrimPi0AccYeta) ;
321 
322  Int_t nphibinsopen = TMath::Nint(nphibins*TMath::TwoPi()/(phimax-phimin));
323  fhPrimPi0Phi = new TH2F("hPrimPi0Phi","#phi of primary #pi^{0}, |#it{Y}|<1",nptbins,ptmin,ptmax,nphibinsopen,0,360) ;
324  fhPrimPi0Phi->SetYTitle("#phi (deg)");
325  fhPrimPi0Phi->SetXTitle("#it{p}_{T} (GeV/#it{c})");
326  outputContainer->Add(fhPrimPi0Phi) ;
327 
328  fhPrimPi0AccPhi = new TH2F("hPrimPi0AccPhi","#phi of primary #pi^{0} with accepted daughters",nptbins,ptmin,ptmax,
329  nphibins,phimin*TMath::RadToDeg(),phimax*TMath::RadToDeg()) ;
330  fhPrimPi0AccPhi->SetYTitle("#phi (deg)");
331  fhPrimPi0AccPhi->SetXTitle("#it{p}_{T} (GeV/#it{c})");
332  outputContainer->Add(fhPrimPi0AccPhi) ;
333 
334  // Eta
335 
336  fhPrimEtaE = new TH1F("hPrimEtaE","Primary eta E",nptbins,ptmin,ptmax) ;
337  fhPrimEtaAccE = new TH1F("hPrimEtaAccE","Primary #eta #it{E} with both photons in acceptance",nptbins,ptmin,ptmax) ;
338  fhPrimEtaE ->SetXTitle("#it{E} (GeV)");
339  fhPrimEtaAccE->SetXTitle("#it{E} (GeV)");
340  outputContainer->Add(fhPrimEtaE) ;
341  outputContainer->Add(fhPrimEtaAccE) ;
342 
343  fhPrimEtaPt = new TH1F("hPrimEtaPt","Primary #eta #it{p}_{T}",nptbins,ptmin,ptmax) ;
344  fhPrimEtaAccPt = new TH1F("hPrimEtaAccPt","Primary eta #it{p}_{T} with both photons in acceptance",nptbins,ptmin,ptmax) ;
345  fhPrimEtaAccPtPhotonCuts = new TH1F("hPrimEtaAccPtPhotonCuts","Primary eta #it{p}_{T} with both photons in acceptance",nptbins,ptmin,ptmax) ;
346  fhPrimEtaPt ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
347  fhPrimEtaAccPt->SetXTitle("#it{p}_{T} (GeV/#it{c})");
348  fhPrimEtaAccPtPhotonCuts->SetXTitle("#it{p}_{T} (GeV/#it{c})");
349  outputContainer->Add(fhPrimEtaPt) ;
350  outputContainer->Add(fhPrimEtaAccPt) ;
351  outputContainer->Add(fhPrimEtaAccPtPhotonCuts) ;
352 
353  fhPrimEtaY = new TH2F("hPrimEtaRapidity","Rapidity of primary #eta",nptbins,ptmin,ptmax,netabinsopen,-2, 2) ;
354  fhPrimEtaY->SetYTitle("#it{Rapidity}");
355  fhPrimEtaY->SetXTitle("#it{p}_{T} (GeV/#it{c})");
356  outputContainer->Add(fhPrimEtaY) ;
357 
358  fhPrimEtaYeta = new TH2F("hPrimEtaPseudoRapidityEta","PseudoRapidity of primary #eta",nptbins,ptmin,ptmax,netabinsopen,-2, 2) ;
359  fhPrimEtaYeta->SetYTitle("#it{Rapidity}");
360  fhPrimEtaYeta->SetXTitle("#it{p}_{T} (GeV/#it{c})");
361  outputContainer->Add(fhPrimEtaYeta) ;
362 
363  fhPrimEtaYetaYcut = new TH2F("hPrimEtaPseudoRapidityEtaYcut","PseudoRapidity of primary #eta, |#it{Y}|<1",nptbins,ptmin,ptmax,netabinsopen,-2, 2) ;
364  fhPrimEtaYetaYcut->SetYTitle("#it{Pseudorapidity}");
365  fhPrimEtaYetaYcut->SetXTitle("#it{p}_{T} (GeV/#it{c})");
366  outputContainer->Add(fhPrimEtaYetaYcut) ;
367 
368  fhPrimEtaAccY = new TH2F("hPrimEtaAccRapidity","Rapidity of primary #eta",nptbins,ptmin,ptmax, netabins,etamin,etamax) ;
369  fhPrimEtaAccY->SetYTitle("#it{Rapidity}");
370  fhPrimEtaAccY->SetXTitle("#it{p}_{T} (GeV/#it{c})");
371  outputContainer->Add(fhPrimEtaAccY) ;
372 
373  fhPrimEtaAccYeta = new TH2F("hPrimEtaAccPseudoRapidity","PseudoRapidity of primary #eta",nptbins,ptmin,ptmax, netabins,etamin,etamax) ;
374  fhPrimEtaAccYeta->SetYTitle("#it{Pseudorapidity}");
375  fhPrimEtaAccYeta->SetXTitle("#it{p}_{T} (GeV/#it{c})");
376  outputContainer->Add(fhPrimEtaAccYeta) ;
377 
378  fhPrimEtaPhi = new TH2F("hPrimEtaPhi","Azimuthal of primary #eta",nptbins,ptmin,ptmax, nphibinsopen,0,360) ;
379  fhPrimEtaPhi->SetYTitle("#phi (deg)");
380  fhPrimEtaPhi->SetXTitle("#it{p}_{T} (GeV/#it{c})");
381  outputContainer->Add(fhPrimEtaPhi) ;
382 
383  fhPrimEtaAccPhi = new TH2F("hPrimEtaAccPhi","Azimuthal of primary #eta with accepted daughters",nptbins,ptmin,ptmax, nphibins,phimin*TMath::RadToDeg(),phimax*TMath::RadToDeg()) ;
384  fhPrimEtaAccPhi->SetYTitle("#phi (deg)");
385  fhPrimEtaAccPhi->SetXTitle("#it{p}_{T} (GeV/#it{c})");
386  outputContainer->Add(fhPrimEtaAccPhi) ;
387 
388  // Create histograms only for PbPb or high multiplicity analysis analysis
390  {
391  fhPrimPi0PtCentrality = new TH2F("hPrimPi0PtCentrality","Primary #pi^{0} #it{p}_{T} vs reco centrality, |#it{Y}|<1",
392  nptbins,ptmin,ptmax, 100, 0, 100) ;
393  fhPrimPi0AccPtCentrality = new TH2F("hPrimPi0AccPtCentrality","Primary #pi^{0} with both photons in acceptance #it{p}_{T} vs reco centrality",
394  nptbins,ptmin,ptmax, 100, 0, 100) ;
395  fhPrimPi0PtCentrality ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
396  fhPrimPi0AccPtCentrality->SetXTitle("#it{p}_{T} (GeV/#it{c})");
397  fhPrimPi0PtCentrality ->SetYTitle("Centrality");
398  fhPrimPi0AccPtCentrality->SetYTitle("Centrality");
399  outputContainer->Add(fhPrimPi0PtCentrality) ;
400  outputContainer->Add(fhPrimPi0AccPtCentrality) ;
401 
402  fhPrimEtaPtCentrality = new TH2F("hPrimEtaPtCentrality","Primary #eta #it{p}_{T} vs reco centrality, |#it{Y}|<1",nptbins,ptmin,ptmax, 100, 0, 100) ;
403  fhPrimEtaAccPtCentrality = new TH2F("hPrimEtaAccPtCentrality","Primary #eta with both photons in acceptance #it{p}_{T} vs reco centrality",nptbins,ptmin,ptmax, 100, 0, 100) ;
404  fhPrimEtaPtCentrality ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
405  fhPrimEtaAccPtCentrality->SetXTitle("#it{p}_{T} (GeV/#it{c})");
406  fhPrimEtaPtCentrality ->SetYTitle("Centrality");
407  fhPrimEtaAccPtCentrality->SetYTitle("Centrality");
408  outputContainer->Add(fhPrimEtaPtCentrality) ;
409  outputContainer->Add(fhPrimEtaAccPtCentrality) ;
410 
411  fhPrimPi0PtEventPlane = new TH2F("hPrimPi0PtEventPlane","Primary #pi^{0} #it{p}_{T} vs reco event plane angle, |#it{Y}|<1",
412  nptbins,ptmin,ptmax, 100, 0, TMath::Pi()) ;
413  fhPrimPi0AccPtEventPlane = new TH2F("hPrimPi0AccPtEventPlane","Primary #pi^{0} with both photons in acceptance #it{p}_{T} vs reco event plane angle",
414  nptbins,ptmin,ptmax, 100, 0, TMath::Pi()) ;
415  fhPrimPi0PtEventPlane ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
416  fhPrimPi0AccPtEventPlane->SetXTitle("#it{p}_{T} (GeV/#it{c})");
417  fhPrimPi0PtEventPlane ->SetYTitle("Event Plane Angle (rad)");
418  fhPrimPi0AccPtEventPlane->SetYTitle("Event Plane Angle (rad)");
419  outputContainer->Add(fhPrimPi0PtEventPlane) ;
420  outputContainer->Add(fhPrimPi0AccPtEventPlane) ;
421 
422  fhPrimEtaPtEventPlane = new TH2F("hPrimEtaPtEventPlane","Primary #eta #it{p}_{T} vs reco event plane angle, |#it{Y}|<1",nptbins,ptmin,ptmax, 100, 0, TMath::Pi()) ;
423  fhPrimEtaAccPtEventPlane = new TH2F("hPrimEtaAccPtEventPlane","Primary #eta with both #gamma_{decay} in acceptance #it{p}_{T} vs reco event plane angle",nptbins,ptmin,ptmax, 100, 0, TMath::Pi()) ;
424  fhPrimEtaPtEventPlane ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
425  fhPrimEtaAccPtEventPlane->SetXTitle("#it{p}_{T} (GeV/#it{c})");
426  fhPrimEtaPtEventPlane ->SetYTitle("Event Plane Angle (rad)");
427  fhPrimEtaAccPtEventPlane->SetYTitle("Event Plane Angle (rad)");
428  outputContainer->Add(fhPrimEtaPtEventPlane) ;
429  outputContainer->Add(fhPrimEtaAccPtEventPlane) ;
430  }
431 
432  if(fFillAngleHisto)
433  {
434  fhPrimPi0OpeningAngle = new TH2F
435  ("hPrimPi0OpeningAngle","Angle between all primary #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,200,0,0.7);
436  fhPrimPi0OpeningAngle->SetYTitle("#theta(rad)");
437  fhPrimPi0OpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
438  outputContainer->Add(fhPrimPi0OpeningAngle) ;
439 
441  ("hPrimPi0OpeningAnglePhotonCuts","Angle between all primary #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,200,0,0.7);
442  fhPrimPi0OpeningAnglePhotonCuts->SetYTitle("#theta(rad)");
443  fhPrimPi0OpeningAnglePhotonCuts->SetXTitle("E_{ #pi^{0}} (GeV)");
444  outputContainer->Add(fhPrimPi0OpeningAnglePhotonCuts) ;
445 
446  fhPrimPi0OpeningAngleAsym = new TH2F
447  ("hPrimPi0OpeningAngleAsym","Angle between all primary #gamma pair vs #it{Asymmetry}, #it{p}_{T}>5 GeV/#it{c}",100,0,1,200,0,0.7);
448  fhPrimPi0OpeningAngleAsym->SetXTitle("|A|=| (E_{1}-E_{2}) / (E_{1}+E_{2}) |");
449  fhPrimPi0OpeningAngleAsym->SetYTitle("#theta(rad)");
450  outputContainer->Add(fhPrimPi0OpeningAngleAsym) ;
451 
452  fhPrimPi0CosOpeningAngle = new TH2F
453  ("hPrimPi0CosOpeningAngle","Cosinus of angle between all primary #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,100,-1,1);
454  fhPrimPi0CosOpeningAngle->SetYTitle("cos (#theta) ");
455  fhPrimPi0CosOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
456  outputContainer->Add(fhPrimPi0CosOpeningAngle) ;
457 
458  fhPrimEtaOpeningAngle = new TH2F
459  ("hPrimEtaOpeningAngle","Angle between all primary #gamma pair vs E_{#eta}",nptbins,ptmin,ptmax,200,0,0.7);
460  fhPrimEtaOpeningAngle->SetYTitle("#theta(rad)");
461  fhPrimEtaOpeningAngle->SetXTitle("E_{#eta} (GeV)");
462  outputContainer->Add(fhPrimEtaOpeningAngle) ;
463 
465  ("hPrimEtaOpeningAnglePhotonCuts","Angle between all primary #gamma pair vs E_{#eta}",nptbins,ptmin,ptmax,200,0,0.7);
466  fhPrimEtaOpeningAnglePhotonCuts->SetYTitle("#theta(rad)");
467  fhPrimEtaOpeningAnglePhotonCuts->SetXTitle("E_{#eta} (GeV)");
468  outputContainer->Add(fhPrimEtaOpeningAnglePhotonCuts) ;
469 
470  fhPrimEtaOpeningAngleAsym = new TH2F
471  ("hPrimEtaOpeningAngleAsym","Angle between all primary #gamma pair vs #it{Asymmetry}, #it{p}_{T}>5 GeV/#it{c}",100,0,1,200,0,0.7);
472  fhPrimEtaOpeningAngleAsym->SetXTitle("|#it{A}|=| (#it{E}_{1}-#it{E}_{2}) / (#it{E}_{1}+#it{E}_{2}) |");
473  fhPrimEtaOpeningAngleAsym->SetYTitle("#theta(rad)");
474  outputContainer->Add(fhPrimEtaOpeningAngleAsym) ;
475 
476 
477  fhPrimEtaCosOpeningAngle = new TH2F
478  ("hPrimEtaCosOpeningAngle","Cosinus of angle between all primary #gamma pair vs E_{#eta}",nptbins,ptmin,ptmax,100,-1,1);
479  fhPrimEtaCosOpeningAngle->SetYTitle("cos (#theta) ");
480  fhPrimEtaCosOpeningAngle->SetXTitle("#it{E}_{ #eta} (GeV)");
481  outputContainer->Add(fhPrimEtaCosOpeningAngle) ;
482  }
483 
484  // Primary origin
485  if(fFillOriginHisto)
486  {
487  // Pi0
488  fhPrimPi0PtOrigin = new TH2F("hPrimPi0PtOrigin","Primary #pi^{0} #it{p}_{T} vs origin",nptbins,ptmin,ptmax,11,0,11) ;
489  fhPrimPi0PtOrigin->SetXTitle("#it{p}_{T} (GeV/#it{c})");
490  fhPrimPi0PtOrigin->SetYTitle("Origin");
491  fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(1 ,"Status 21");
492  fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(2 ,"Quark");
493  fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(3 ,"qq Resonances ");
494  fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(4 ,"Resonances");
495  fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(5 ,"#rho");
496  fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(6 ,"#omega");
497  fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(7 ,"K");
498  fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(8 ,"Other");
499  fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(9 ,"#eta");
500  fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(10 ,"#eta prime");
501  outputContainer->Add(fhPrimPi0PtOrigin) ;
502 
503  fhPrimNotResonancePi0PtOrigin = new TH2F("hPrimNotResonancePi0PtOrigin","Primary #pi^{0} #it{p}_{T} vs origin",nptbins,ptmin,ptmax,11,0,11) ;
504  fhPrimNotResonancePi0PtOrigin->SetXTitle("#it{p}_{T} (GeV/#it{c})");
505  fhPrimNotResonancePi0PtOrigin->SetYTitle("Origin");
506  fhPrimNotResonancePi0PtOrigin->GetYaxis()->SetBinLabel(1 ,"Status 21");
507  fhPrimNotResonancePi0PtOrigin->GetYaxis()->SetBinLabel(2 ,"Quark");
508  fhPrimNotResonancePi0PtOrigin->GetYaxis()->SetBinLabel(3 ,"qq Resonances");
509  fhPrimNotResonancePi0PtOrigin->GetYaxis()->SetBinLabel(4 ,"Resonances");
510  fhPrimNotResonancePi0PtOrigin->GetYaxis()->SetBinLabel(5 ,"#rho");
511  fhPrimNotResonancePi0PtOrigin->GetYaxis()->SetBinLabel(6 ,"#omega");
512  fhPrimNotResonancePi0PtOrigin->GetYaxis()->SetBinLabel(7 ,"K");
513  fhPrimNotResonancePi0PtOrigin->GetYaxis()->SetBinLabel(8 ,"Other");
514  fhPrimNotResonancePi0PtOrigin->GetYaxis()->SetBinLabel(9 ,"#eta");
515  fhPrimNotResonancePi0PtOrigin->GetYaxis()->SetBinLabel(10 ,"#eta prime");
516  outputContainer->Add(fhPrimNotResonancePi0PtOrigin) ;
517 
518  fhPrimPi0PtStatus = new TH2F("hPrimPi0PtStatus","Primary #pi^{0} #it{p}_{T} vs status",nptbins,ptmin,ptmax,101,-50,50) ;
519  fhPrimPi0PtStatus->SetXTitle("#it{p}_{T} (GeV/#it{c})");
520  fhPrimPi0PtStatus->SetYTitle("Status");
521  outputContainer->Add(fhPrimPi0PtStatus) ;
522 
523  // Eta
524  fhPrimEtaPtOrigin = new TH2F("hPrimEtaPtOrigin","Primary #pi^{0} #it{p}_{T} vs origin",nptbins,ptmin,ptmax,7,0,7) ;
525  fhPrimEtaPtOrigin->SetXTitle("#it{p}_{T} (GeV/#it{c})");
526  fhPrimEtaPtOrigin->SetYTitle("Origin");
527  fhPrimEtaPtOrigin->GetYaxis()->SetBinLabel(1 ,"Status 21");
528  fhPrimEtaPtOrigin->GetYaxis()->SetBinLabel(2 ,"Quark");
529  fhPrimEtaPtOrigin->GetYaxis()->SetBinLabel(3 ,"qq Resonances");
530  fhPrimEtaPtOrigin->GetYaxis()->SetBinLabel(4 ,"Resonances");
531  fhPrimEtaPtOrigin->GetYaxis()->SetBinLabel(5 ,"Other");
532  fhPrimEtaPtOrigin->GetYaxis()->SetBinLabel(6 ,"#eta prime ");
533  outputContainer->Add(fhPrimEtaPtOrigin) ;
534 
535  // Production vertex
536  fhPrimPi0ProdVertex = new TH2F("hPrimPi0ProdVertex","generated #pi^{0} #it{p}_{T} vs production vertex",
537  200,0.,20.,5000,0,500) ;
538  fhPrimPi0ProdVertex->SetXTitle("#it{p}_{T} (GeV/#it{c})");
539  fhPrimPi0ProdVertex->SetYTitle("#it{R} (cm)");
540  outputContainer->Add(fhPrimPi0ProdVertex) ;
541 
542  fhPrimEtaProdVertex = new TH2F("hPrimEtaProdVertex","generated #eta #it{p}_{T} vs production vertex",
543  200,0.,20.,5000,0,500) ;
544  fhPrimEtaProdVertex->SetXTitle("#it{p}_{T} (GeV/#it{c})");
545  fhPrimEtaProdVertex->SetYTitle("#it{R} (cm)");
546  outputContainer->Add(fhPrimEtaProdVertex) ;
547  }
548 
550  {
551  TString ebin[] = {"8 < E < 12 GeV","12 < E < 16 GeV", "16 < E < 20 GeV", "E > 20 GeV" };
552  Int_t narmbins = 400;
553  Float_t armmin = 0;
554  Float_t armmax = 0.4;
555 
556  for(Int_t i = 0; i < 4; i++)
557  {
558  fhArmPrimPi0[i] = new TH2F(Form("hArmenterosPrimPi0EBin%d",i),
559  Form("Armenteros of primary #pi^{0}, %s",ebin[i].Data()),
560  200, -1, 1, narmbins,armmin,armmax);
561  fhArmPrimPi0[i]->SetYTitle("#it{p}_{T}^{Arm}");
562  fhArmPrimPi0[i]->SetXTitle("#alpha^{Arm}");
563  outputContainer->Add(fhArmPrimPi0[i]) ;
564 
565  fhArmPrimEta[i] = new TH2F(Form("hArmenterosPrimEtaEBin%d",i),
566  Form("Armenteros of primary #eta, %s",ebin[i].Data()),
567  200, -1, 1, narmbins,armmin,armmax);
568  fhArmPrimEta[i]->SetYTitle("#it{p}_{T}^{Arm}");
569  fhArmPrimEta[i]->SetXTitle("#alpha^{Arm}");
570  outputContainer->Add(fhArmPrimEta[i]) ;
571  }
572 
573  // Same as asymmetry ...
574  fhCosThStarPrimPi0 = new TH2F
575  ("hCosThStarPrimPi0","cos(#theta *) for primary #pi^{0}",nptbins,ptmin,ptmax,200,-1,1);
576  fhCosThStarPrimPi0->SetYTitle("cos(#theta *)");
577  fhCosThStarPrimPi0->SetXTitle("E_{ #pi^{0}} (GeV)");
578  outputContainer->Add(fhCosThStarPrimPi0) ;
579 
580  fhCosThStarPrimEta = new TH2F
581  ("hCosThStarPrimEta","cos(#theta *) for primary #eta",nptbins,ptmin,ptmax,200,-1,1);
582  fhCosThStarPrimEta->SetYTitle("cos(#theta *)");
583  fhCosThStarPrimEta->SetXTitle("E_{ #eta} (GeV)");
584  outputContainer->Add(fhCosThStarPrimEta) ;
585  }
586 
587  if(fFillOnlyMCAcceptanceHisto) return outputContainer;
588  }
589 
590  //
591  // Create mixed event containers
592  //
593  fEventsList = new TList*[GetNCentrBin()*GetNZvertBin()*GetNRPBin()] ;
594 
595  for(Int_t ic=0; ic<GetNCentrBin(); ic++)
596  {
597  for(Int_t iz=0; iz<GetNZvertBin(); iz++)
598  {
599  for(Int_t irp=0; irp<GetNRPBin(); irp++)
600  {
601  Int_t bin = GetEventMixBin(ic,iz,irp);
602  fEventsList[bin] = new TList() ;
603  fEventsList[bin]->SetOwner(kFALSE);
604  }
605  }
606  }
607 
608  // Init the number of modules, set in the class AliCalorimeterUtils
610  if(GetCalorimeter()==kPHOS && fNModules > 4) fNModules = 4;
611 
612  fhReMod = new TH2F*[fNModules] ;
613  fhMiMod = new TH2F*[fNModules] ;
614 
616  {
617  if(GetCalorimeter() == kPHOS)
618  {
619  fhReDiffPHOSMod = new TH2F*[fNModules] ;
620  fhMiDiffPHOSMod = new TH2F*[fNModules] ;
621  }
622  else
623  {
624  fhReSameSectorEMCALMod = new TH2F*[fNModules/2] ;
625  fhReSameSideEMCALMod = new TH2F*[fNModules-2] ;
626  fhMiSameSectorEMCALMod = new TH2F*[fNModules/2] ;
627  fhMiSameSideEMCALMod = new TH2F*[fNModules-2] ;
628  }
629  }
630  else
631  {
632  fhReSameSectorDCALPHOSMod = new TH2F*[6] ;
633  fhReDiffSectorDCALPHOSMod = new TH2F*[8] ;
634  fhMiSameSectorDCALPHOSMod = new TH2F*[6] ;
635  fhMiDiffSectorDCALPHOSMod = new TH2F*[8] ;
636  }
637 
638  fhRe1 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
639  fhMi1 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
641  {
642  fhRe2 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
643  fhRe3 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
644  fhMi2 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
645  fhMi3 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
646  }
647  if(fMakeInvPtPlots)
648  {
651  if(fFillBadDistHisto){
656  }
657  }
658 
659  if(fCheckConversion)
660  {
661  fhReConv = new TH2F("hReConv","Real Pair with one recombined conversion ",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
662  fhReConv->SetXTitle("#it{p}_{T} (GeV/#it{c})");
663  fhReConv->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
664  outputContainer->Add(fhReConv) ;
665 
666  fhReConv2 = new TH2F("hReConv2","Real Pair with 2 recombined conversion ",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
667  fhReConv2->SetXTitle("#it{p}_{T} (GeV/#it{c})");
668  fhReConv2->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
669  outputContainer->Add(fhReConv2) ;
670 
671  if(DoOwnMix())
672  {
673  fhMiConv = new TH2F("hMiConv","Mixed Pair with one recombined conversion ",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
674  fhMiConv->SetXTitle("#it{p}_{T} (GeV/#it{c})");
675  fhMiConv->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
676  outputContainer->Add(fhMiConv) ;
677 
678  fhMiConv2 = new TH2F("hMiConv2","Mixed Pair with 2 recombined conversion ",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
679  fhMiConv2->SetXTitle("#it{p}_{T} (GeV/#it{c})");
680  fhMiConv2->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
681  outputContainer->Add(fhMiConv2) ;
682  }
683  }
684 
685  for(Int_t ic=0; ic<GetNCentrBin(); ic++)
686  {
687  for(Int_t ipid=0; ipid<fNPIDBits; ipid++)
688  {
689  for(Int_t iasym=0; iasym<fNAsymCuts; iasym++)
690  {
691  Int_t index = ((ic*fNPIDBits)+ipid)*fNAsymCuts + iasym;
692  //printf("cen %d, pid %d, asy %d, Index %d\n",ic,ipid,iasym,index);
693  //Distance to bad module 1
694  snprintf(key, buffersize,"hRe_cen%d_pidbit%d_asy%d_dist1",ic,ipid,iasym) ;
695  snprintf(title, buffersize,"Real #it{M}_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 1",
696  ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
697  fhRe1[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
698  fhRe1[index]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
699  fhRe1[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
700  //printf("name: %s\n ",fhRe1[index]->GetName());
701  outputContainer->Add(fhRe1[index]) ;
702 
704  {
705  //Distance to bad module 2
706  snprintf(key, buffersize,"hRe_cen%d_pidbit%d_asy%d_dist2",ic,ipid,iasym) ;
707  snprintf(title, buffersize,"Real #it{M}_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 2",
708  ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
709  fhRe2[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
710  fhRe2[index]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
711  fhRe2[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
712  outputContainer->Add(fhRe2[index]) ;
713 
714  //Distance to bad module 3
715  snprintf(key, buffersize,"hRe_cen%d_pidbit%d_asy%d_dist3",ic,ipid,iasym) ;
716  snprintf(title, buffersize,"Real #it{M}_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 3",
717  ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
718  fhRe3[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
719  fhRe3[index]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
720  fhRe3[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
721  outputContainer->Add(fhRe3[index]) ;
722  }
723 
724  //Inverse pT
725  if(fMakeInvPtPlots)
726  {
727  //Distance to bad module 1
728  snprintf(key, buffersize,"hReInvPt_cen%d_pidbit%d_asy%d_dist1",ic,ipid,iasym) ;
729  snprintf(title, buffersize,"Real #it{M}_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 1",
730  ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
731  fhReInvPt1[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
732  fhReInvPt1[index]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
733  fhReInvPt1[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
734  outputContainer->Add(fhReInvPt1[index]) ;
735 
736  if(fFillBadDistHisto){
737  //Distance to bad module 2
738  snprintf(key, buffersize,"hReInvPt_cen%d_pidbit%d_asy%d_dist2",ic,ipid,iasym) ;
739  snprintf(title, buffersize,"Real #it{M}_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 2",
740  ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
741  fhReInvPt2[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
742  fhReInvPt2[index]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
743  fhReInvPt2[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
744  outputContainer->Add(fhReInvPt2[index]) ;
745 
746  //Distance to bad module 3
747  snprintf(key, buffersize,"hReInvPt_cen%d_pidbit%d_asy%d_dist3",ic,ipid,iasym) ;
748  snprintf(title, buffersize,"Real #it{M}_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 3",
749  ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
750  fhReInvPt3[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
751  fhReInvPt3[index]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
752  fhReInvPt3[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
753  outputContainer->Add(fhReInvPt3[index]) ;
754  }
755  }
756 
757  if(DoOwnMix())
758  {
759  //Distance to bad module 1
760  snprintf(key, buffersize,"hMi_cen%d_pidbit%d_asy%d_dist1",ic,ipid,iasym) ;
761  snprintf(title, buffersize,"Mixed #it{M}_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 1",
762  ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
763  fhMi1[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
764  fhMi1[index]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
765  fhMi1[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
766  outputContainer->Add(fhMi1[index]) ;
767  if(fFillBadDistHisto){
768  //Distance to bad module 2
769  snprintf(key, buffersize,"hMi_cen%d_pidbit%d_asy%d_dist2",ic,ipid,iasym) ;
770  snprintf(title, buffersize,"Mixed #it{M}_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 2",
771  ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
772  fhMi2[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
773  fhMi2[index]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
774  fhMi2[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
775  outputContainer->Add(fhMi2[index]) ;
776 
777  //Distance to bad module 3
778  snprintf(key, buffersize,"hMi_cen%d_pidbit%d_asy%d_dist3",ic,ipid,iasym) ;
779  snprintf(title, buffersize,"Mixed #it{M}_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 3",
780  ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
781  fhMi3[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
782  fhMi3[index]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
783  fhMi3[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
784  outputContainer->Add(fhMi3[index]) ;
785  }
786 
787  // Inverse pT
788  if(fMakeInvPtPlots)
789  {
790  // Distance to bad module 1
791  snprintf(key, buffersize,"hMiInvPt_cen%d_pidbit%d_asy%d_dist1",ic,ipid,iasym) ;
792  snprintf(title, buffersize,"Mixed #it{M}_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 1",
793  ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
794  fhMiInvPt1[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
795  fhMiInvPt1[index]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
796  fhMiInvPt1[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
797  outputContainer->Add(fhMiInvPt1[index]) ;
798  if(fFillBadDistHisto){
799  // Distance to bad module 2
800  snprintf(key, buffersize,"hMiInvPt_cen%d_pidbit%d_asy%d_dist2",ic,ipid,iasym) ;
801  snprintf(title, buffersize,"Mixed #it{M}_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 2",
802  ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
803  fhMiInvPt2[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
804  fhMiInvPt2[index]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
805  fhMiInvPt2[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
806  outputContainer->Add(fhMiInvPt2[index]) ;
807 
808  // Distance to bad module 3
809  snprintf(key, buffersize,"hMiInvPt_cen%d_pidbit%d_asy%d_dist3",ic,ipid,iasym) ;
810  snprintf(title, buffersize,"Mixed #it{M}_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f,dist bad 3",
811  ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
812  fhMiInvPt3[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
813  fhMiInvPt3[index]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
814  fhMiInvPt3[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
815  outputContainer->Add(fhMiInvPt3[index]) ;
816  }
817  }
818  }
819  }
820  }
821  }
822 
823  fhEPairDiffTime = new TH2F("hEPairDiffTime","cluster pair time difference vs #it{p}_{T}",nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
824  fhEPairDiffTime->SetXTitle("#it{p}_{T,pair} (GeV/#it{c})");
825  fhEPairDiffTime->SetYTitle("#Delta t (ns)");
826  outputContainer->Add(fhEPairDiffTime);
827 
829  {
830  fhRePtAsym = new TH2F("hRePtAsym","#it{Asymmetry} vs #it{p}_{T} , for pairs",nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax) ;
831  fhRePtAsym->SetXTitle("#it{p}_{T} (GeV/#it{c})");
832  fhRePtAsym->SetYTitle("#it{Asymmetry}");
833  outputContainer->Add(fhRePtAsym);
834 
835  fhRePtAsymPi0 = new TH2F("hRePtAsymPi0","#it{Asymmetry} vs #it{p}_{T} , for pairs close to #pi^{0} mass",nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax) ;
836  fhRePtAsymPi0->SetXTitle("#it{p}_{T} (GeV/#it{c})");
837  fhRePtAsymPi0->SetYTitle("Asymmetry");
838  outputContainer->Add(fhRePtAsymPi0);
839 
840  fhRePtAsymEta = new TH2F("hRePtAsymEta","#it{Asymmetry} vs #it{p}_{T} , for pairs close to #eta mass",nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax) ;
841  fhRePtAsymEta->SetXTitle("#it{p}_{T} (GeV/#it{c})");
842  fhRePtAsymEta->SetYTitle("#it{Asymmetry}");
843  outputContainer->Add(fhRePtAsymEta);
844  }
845 
846  if(fMultiCutAna)
847  {
848  fhRePIDBits = new TH2F*[fNPIDBits];
849  for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
850  snprintf(key, buffersize,"hRe_pidbit%d",ipid) ;
851  snprintf(title, buffersize,"Real #it{M}_{#gamma#gamma} distr. for PIDBit=%d",fPIDBits[ipid]) ;
852  fhRePIDBits[ipid] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
853  fhRePIDBits[ipid]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
854  fhRePIDBits[ipid]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
855  outputContainer->Add(fhRePIDBits[ipid]) ;
856  }// pid bit loop
857 
860 
862  for(Int_t iSM = 0; iSM < fNModules; iSM++) fhRePtNCellAsymCutsSM[iSM] = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
863 
864 
865  for(Int_t ipt=0; ipt<fNPtCuts; ipt++)
866  {
867  for(Int_t icell=0; icell<fNCellNCuts; icell++)
868  {
869  for(Int_t iasym=0; iasym<fNAsymCuts; iasym++)
870  {
871  snprintf(key, buffersize,"hRe_pt%d_cell%d_asym%d",ipt,icell,iasym) ;
872  snprintf(title, buffersize,"Real #it{M}_{#gamma#gamma} distr. for pt >%2.2f, ncell>%d and asym >%1.2f ",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]) ;
873  Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
874  //printf("ipt %d, icell %d, iassym %d, index %d\n",ipt, icell, iasym, index);
875  fhRePtNCellAsymCuts[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
876  fhRePtNCellAsymCuts[index]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
877  fhRePtNCellAsymCuts[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
878  outputContainer->Add(fhRePtNCellAsymCuts[index]) ;
879 
880  snprintf(key, buffersize,"hMi_pt%d_cell%d_asym%d",ipt,icell,iasym) ;
881  snprintf(title, buffersize,"Mixed #it{M}_{#gamma#gamma} distr. for pt >%2.2f, ncell>%d and asym >%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]) ;
882  fhMiPtNCellAsymCuts[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
883  fhMiPtNCellAsymCuts[index]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
884  fhMiPtNCellAsymCuts[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
885  outputContainer->Add(fhMiPtNCellAsymCuts[index]) ;
886 
888  {
889  for(Int_t iSM = 0; iSM < fNModules; iSM++)
890  {
891  snprintf(key, buffersize,"hRe_pt%d_cell%d_asym%d_SM%d",ipt,icell,iasym,iSM) ;
892  snprintf(title, buffersize,"Real #it{M}_{#gamma#gamma} distr. for pt >%2.2f, ncell>%d and asym >%1.2f, SM %d ",
893  fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym],iSM) ;
894  fhRePtNCellAsymCutsSM[iSM][index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
895  fhRePtNCellAsymCutsSM[iSM][index]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
896  fhRePtNCellAsymCutsSM[iSM][index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
897  outputContainer->Add(fhRePtNCellAsymCutsSM[iSM][index]) ;
898  }
899  }
900  }
901  }
902  }
903 
904  if(ntrmbins!=0)
905  {
906  fhRePtMult = new TH3F*[fNAsymCuts] ;
907  for(Int_t iasym = 0; iasym<fNAsymCuts; iasym++)
908  {
909  fhRePtMult[iasym] = new TH3F(Form("hRePtMult_asym%d",iasym),Form("(#it{p}_{T},C,M)_{#gamma#gamma}, A<%1.2f",fAsymCuts[iasym]),
910  nptbins,ptmin,ptmax,ntrmbins,ntrmmin,ntrmmax,nmassbins,massmin,massmax);
911  fhRePtMult[iasym]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
912  fhRePtMult[iasym]->SetYTitle("Track multiplicity");
913  fhRePtMult[iasym]->SetZTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
914  outputContainer->Add(fhRePtMult[iasym]) ;
915  }
916  }
917  } // multi cuts analysis
918 
920  {
921  fhReSS[0] = new TH2F("hRe_SS_Tight"," 0.01 < #lambda_{0}^{2} < 0.4",
922  nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
923  fhReSS[0]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
924  fhReSS[0]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
925  outputContainer->Add(fhReSS[0]) ;
926 
927 
928  fhReSS[1] = new TH2F("hRe_SS_Loose"," #lambda_{0}^{2} > 0.4",
929  nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
930  fhReSS[1]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
931  fhReSS[1]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
932  outputContainer->Add(fhReSS[1]) ;
933 
934 
935  fhReSS[2] = new TH2F("hRe_SS_Both"," cluster_{1} #lambda_{0}^{2} > 0.4; cluster_{2} 0.01 < #lambda_{0}^{2} < 0.4",
936  nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
937  fhReSS[2]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
938  fhReSS[2]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
939  outputContainer->Add(fhReSS[2]) ;
940  }
941 
942  if(DoOwnMix())
943  {
944  fhEventBin=new TH1I("hEventBin","Number of real pairs per bin(cen,vz,rp)",
947  fhEventBin->SetXTitle("bin");
948  outputContainer->Add(fhEventBin) ;
949 
950 
951  fhEventMixBin=new TH1I("hEventMixBin","Number of mixed pairs per bin(cen,vz,rp)",
954  fhEventMixBin->SetXTitle("bin");
955  outputContainer->Add(fhEventMixBin) ;
956  }
957 
959  {
960  fhCentrality=new TH1F("hCentralityBin","Number of events in centrality bin",GetNCentrBin(),0.,1.*GetNCentrBin()) ;
961  fhCentrality->SetXTitle("Centrality bin");
962  outputContainer->Add(fhCentrality) ;
963 
964  fhCentralityNoPair=new TH1F("hCentralityBinNoPair","Number of events in centrality bin, with no cluster pairs",GetNCentrBin(),0.,1.*GetNCentrBin()) ;
965  fhCentralityNoPair->SetXTitle("Centrality bin");
966  outputContainer->Add(fhCentralityNoPair) ;
967 
968  fhEventPlaneResolution=new TH2F("hEventPlaneResolution","Event plane resolution",GetNCentrBin(),0,GetNCentrBin(),100,0.,TMath::TwoPi()) ;
969  fhEventPlaneResolution->SetYTitle("Resolution");
970  fhEventPlaneResolution->SetXTitle("Centrality Bin");
971  outputContainer->Add(fhEventPlaneResolution) ;
972  }
973 
974  if(fFillAngleHisto)
975  {
976  fhRealOpeningAngle = new TH2F
977  ("hRealOpeningAngle","Angle between all #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,300,0,TMath::Pi());
978  fhRealOpeningAngle->SetYTitle("#theta(rad)");
979  fhRealOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
980  outputContainer->Add(fhRealOpeningAngle) ;
981 
982  fhRealCosOpeningAngle = new TH2F
983  ("hRealCosOpeningAngle","Cosinus of angle between all #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,100,0,1);
984  fhRealCosOpeningAngle->SetYTitle("cos (#theta) ");
985  fhRealCosOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
986  outputContainer->Add(fhRealCosOpeningAngle) ;
987 
988  if(DoOwnMix())
989  {
990  fhMixedOpeningAngle = new TH2F
991  ("hMixedOpeningAngle","Angle between all #gamma pair vs E_{#pi^{0}}, Mixed pairs",nptbins,ptmin,ptmax,300,0,TMath::Pi());
992  fhMixedOpeningAngle->SetYTitle("#theta(rad)");
993  fhMixedOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
994  outputContainer->Add(fhMixedOpeningAngle) ;
995 
996  fhMixedCosOpeningAngle = new TH2F
997  ("hMixedCosOpeningAngle","Cosinus of angle between all #gamma pair vs E_{#pi^{0}}, Mixed pairs",nptbins,ptmin,ptmax,100,0,1);
998  fhMixedCosOpeningAngle->SetYTitle("cos (#theta) ");
999  fhMixedCosOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
1000  outputContainer->Add(fhMixedCosOpeningAngle) ;
1001  }
1002  }
1003 
1004  // Histograms filled only if MC data is requested
1005  if(IsDataMC())
1006  {
1007  fhReMCFromConversion = new TH2F("hReMCFromConversion","Invariant mass of 2 clusters originated in conversions",
1008  nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
1009  fhReMCFromConversion->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1010  fhReMCFromConversion->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1011  outputContainer->Add(fhReMCFromConversion) ;
1012 
1013  fhReMCFromNotConversion = new TH2F("hReMCNotFromConversion","Invariant mass of 2 clusters not originated in conversions",
1014  nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
1015  fhReMCFromNotConversion->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1016  fhReMCFromNotConversion->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1017  outputContainer->Add(fhReMCFromNotConversion) ;
1018 
1019  fhReMCFromMixConversion = new TH2F("hReMCFromMixConversion","Invariant mass of 2 clusters one from conversion and the other not",
1020  nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
1021  fhReMCFromMixConversion->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1022  fhReMCFromMixConversion->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1023  outputContainer->Add(fhReMCFromMixConversion) ;
1024 
1025  if(fFillOriginHisto)
1026  {
1027 
1028  fhMCPi0PtOrigin = new TH2F("hMCPi0PtOrigin","Reconstructed pair from generated #pi^{0} #it{p}_{T} vs origin",nptbins,ptmin,ptmax,11,0,11) ;
1029  fhMCPi0PtOrigin->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1030  fhMCPi0PtOrigin->SetYTitle("Origin");
1031  fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(1 ,"Status 21");
1032  fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(2 ,"Quark");
1033  fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(3 ,"qq Resonances");
1034  fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(4 ,"Resonances");
1035  fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(5 ,"#rho");
1036  fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(6 ,"#omega");
1037  fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(7 ,"K");
1038  fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(8 ,"Other");
1039  fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(9 ,"#eta");
1040  fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(10 ,"#eta prime");
1041  outputContainer->Add(fhMCPi0PtOrigin) ;
1042 
1043  fhMCNotResonancePi0PtOrigin = new TH2F("hMCNotResonancePi0PtOrigin","Reconstructed pair from generated #pi^{0} #it{p}_{T} vs origin",nptbins,ptmin,ptmax,11,0,11) ;
1044  fhMCNotResonancePi0PtOrigin->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1045  fhMCNotResonancePi0PtOrigin->SetYTitle("Origin");
1046  fhMCNotResonancePi0PtOrigin->GetYaxis()->SetBinLabel(1 ,"Status 21");
1047  fhMCNotResonancePi0PtOrigin->GetYaxis()->SetBinLabel(2 ,"Quark");
1048  fhMCNotResonancePi0PtOrigin->GetYaxis()->SetBinLabel(3 ,"qq Resonances");
1049  fhMCNotResonancePi0PtOrigin->GetYaxis()->SetBinLabel(4 ,"Resonances");
1050  fhMCNotResonancePi0PtOrigin->GetYaxis()->SetBinLabel(5 ,"#rho");
1051  fhMCNotResonancePi0PtOrigin->GetYaxis()->SetBinLabel(6 ,"#omega");
1052  fhMCNotResonancePi0PtOrigin->GetYaxis()->SetBinLabel(7 ,"K");
1053  fhMCNotResonancePi0PtOrigin->GetYaxis()->SetBinLabel(8 ,"Other");
1054  fhMCNotResonancePi0PtOrigin->GetYaxis()->SetBinLabel(9 ,"#eta");
1055  fhMCNotResonancePi0PtOrigin->GetYaxis()->SetBinLabel(10 ,"#eta prime");
1056  outputContainer->Add(fhMCNotResonancePi0PtOrigin) ;
1057 
1058  fhMCPi0PtStatus = new TH2F("hMCPi0PtStatus","Reconstructed pair from generated #pi^{0} #it{p}_{T} vs status",nptbins,ptmin,ptmax,101,-50,50) ;
1059  fhMCPi0PtStatus->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1060  fhMCPi0PtStatus->SetYTitle("Status");
1061  outputContainer->Add(fhMCPi0PtStatus) ;
1062 
1063  // Eta
1064 
1065  fhMCEtaPtOrigin = new TH2F("hMCEtaPtOrigin","Reconstructed pair from generated #pi^{0} #it{p}_{T} vs origin",nptbins,ptmin,ptmax,7,0,7) ;
1066  fhMCEtaPtOrigin->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1067  fhMCEtaPtOrigin->SetYTitle("Origin");
1068  fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(1 ,"Status 21");
1069  fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(2 ,"Quark");
1070  fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(3 ,"qq Resonances");
1071  fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(4 ,"Resonances");
1072  fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(5 ,"Other");
1073  fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(6 ,"#eta prime");
1074  outputContainer->Add(fhMCEtaPtOrigin) ;
1075 
1076  fhMCPi0ProdVertex = new TH2F("hMCPi0ProdVertex","Selected reco pair from generated #pi^{0} #it{p}_{T} vs production vertex",
1077  200,0.,20.,5000,0,500) ;
1078  fhMCPi0ProdVertex->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1079  fhMCPi0ProdVertex->SetYTitle("#it{R} (cm)");
1080  outputContainer->Add(fhMCPi0ProdVertex) ;
1081 
1082  fhMCEtaProdVertex = new TH2F("hMCEtaProdVertex","Selected reco pair from generated #eta #it{p}_{T} vs production vertex",
1083  200,0.,20.,5000,0,500) ;
1084  fhMCEtaProdVertex->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1085  fhMCEtaProdVertex->SetYTitle("#it{R} (cm)");
1086  outputContainer->Add(fhMCEtaProdVertex) ;
1087 
1088  for(Int_t i = 0; i<13; i++)
1089  {
1090  fhMCOrgMass[i] = new TH2F(Form("hMCOrgMass_%d",i),Form("#it{M} vs #it{p}_{T}, origin %d",i),nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1091  fhMCOrgMass[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1092  fhMCOrgMass[i]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1093  outputContainer->Add(fhMCOrgMass[i]) ;
1094 
1095  fhMCOrgAsym[i]= new TH2F(Form("hMCOrgAsym_%d",i),Form("#it{Asymmetry} vs #it{p}_{T}, origin %d",i),nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax) ;
1096  fhMCOrgAsym[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1097  fhMCOrgAsym[i]->SetYTitle("A");
1098  outputContainer->Add(fhMCOrgAsym[i]) ;
1099 
1100  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) ;
1101  fhMCOrgDeltaEta[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1102  fhMCOrgDeltaEta[i]->SetYTitle("#Delta #eta");
1103  outputContainer->Add(fhMCOrgDeltaEta[i]) ;
1104 
1105  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) ;
1106  fhMCOrgDeltaPhi[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1107  fhMCOrgDeltaPhi[i]->SetYTitle("#Delta #phi (rad)");
1108  outputContainer->Add(fhMCOrgDeltaPhi[i]) ;
1109  }
1110 
1111  if(fMultiCutAnaSim)
1112  {
1119 
1120  for(Int_t ipt=0; ipt<fNPtCuts; ipt++)
1121  {
1122  for(Int_t icell=0; icell<fNCellNCuts; icell++)
1123  {
1124  for(Int_t iasym=0; iasym<fNAsymCuts; iasym++)
1125  {
1126  Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
1127 
1128  fhMCPi0MassPtRec[index] = new TH2F(Form("hMCPi0MassPtRec_pt%d_cell%d_asym%d",ipt,icell,iasym),
1129  Form("Reconstructed #it{M} vs reconstructed #it{p}_{T} of true #pi^{0} cluster pairs for #it{p}_{T} >%2.2f, #it{N}^{cluster}_{cell}>%d and |#it{A}|>%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]),
1130  nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1131  fhMCPi0MassPtRec[index]->SetXTitle("#it{p}_{T, reconstructed} (GeV/#it{c})");
1132  fhMCPi0MassPtRec[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1133  outputContainer->Add(fhMCPi0MassPtRec[index]) ;
1134 
1135  fhMCPi0MassPtTrue[index] = new TH2F(Form("hMCPi0MassPtTrue_pt%d_cell%d_asym%d",ipt,icell,iasym),
1136  Form("Reconstructed #it{M} vs generated #it{p}_{T} of true #pi^{0} cluster pairs for #it{p}_{T} >%2.2f, #it{N}^{cluster}_{cell}>%d and |#it{A}|>%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]),
1137  nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1138  fhMCPi0MassPtTrue[index]->SetXTitle("#it{p}_{T, generated} (GeV/#it{c})");
1139  fhMCPi0MassPtTrue[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1140  outputContainer->Add(fhMCPi0MassPtTrue[index]) ;
1141 
1142  fhMCPi0PtTruePtRec[index] = new TH2F(Form("hMCPi0PtTruePtRec_pt%d_cell%d_asym%d",ipt,icell,iasym),
1143  Form("Generated vs reconstructed #it{p}_{T} of true #pi^{0} cluster pairs, 0.01 < rec. mass < 0.17 MeV/#it{c}^{2} for #it{p}_{T} >%2.2f, #it{N}^{cluster}_{cell}>%d and |#it{A}|>%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]),
1144  nptbins,ptmin,ptmax,nptbins,ptmin,ptmax) ;
1145  fhMCPi0PtTruePtRec[index]->SetXTitle("#it{p}_{T, generated} (GeV/#it{c})");
1146  fhMCPi0PtTruePtRec[index]->SetYTitle("#it{p}_{T, reconstructed} (GeV/#it{c})");
1147  outputContainer->Add(fhMCPi0PtTruePtRec[index]) ;
1148 
1149  fhMCEtaMassPtRec[index] = new TH2F(Form("hMCEtaMassPtRec_pt%d_cell%d_asym%d",ipt,icell,iasym),
1150  Form("Reconstructed #it{M} vs reconstructed #it{p}_{T} of true #eta cluster pairs for #it{p}_{T} >%2.2f, #it{N}^{cluster}_{cell}>%d and |#it{A}|>%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]),
1151  nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1152  fhMCEtaMassPtRec[index]->SetXTitle("#it{p}_{T, generated} (GeV/#it{c})");
1153  fhMCEtaMassPtRec[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1154  outputContainer->Add(fhMCEtaMassPtRec[index]) ;
1155 
1156  fhMCEtaMassPtTrue[index] = new TH2F(Form("hMCEtaMassPtTrue_pt%d_cell%d_asym%d",ipt,icell,iasym),
1157  Form("Reconstructed #it{M} vs generated #it{p}_{T} of true #eta cluster pairs for #it{p}_{T} >%2.2f, #it{N}^{cluster}_{cell}>%d and |#it{A}|>%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]),
1158  nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1159  fhMCEtaMassPtTrue[index]->SetXTitle("#it{p}_{T, generated} (GeV/#it{c})");
1160  fhMCEtaMassPtTrue[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1161  outputContainer->Add(fhMCEtaMassPtTrue[index]) ;
1162 
1163  fhMCEtaPtTruePtRec[index] = new TH2F(Form("hMCEtaPtTruePtRec_pt%d_cell%d_asym%d",ipt,icell,iasym),
1164  Form("Generated vs reconstructed #it{p}_{T} of true #eta cluster pairs, 0.01 < rec. mass < 0.17 MeV/#it{c}^{2} for #it{p}_{T} >%2.2f, ncell>%d and asym >%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]),
1165  nptbins,ptmin,ptmax,nptbins,ptmin,ptmax) ;
1166  fhMCEtaPtTruePtRec[index]->SetXTitle("#it{p}_{T, generated} (GeV/#it{c})");
1167  fhMCEtaPtTruePtRec[index]->SetYTitle("#it{p}_{T, reconstructed} (GeV/#it{c})");
1168  outputContainer->Add(fhMCEtaPtTruePtRec[index]) ;
1169  }
1170  }
1171  }
1172  }//multi cut ana
1173  else
1174  {
1175  fhMCPi0MassPtTrue = new TH2F*[1];
1176  fhMCPi0PtTruePtRec = new TH2F*[1];
1177  fhMCEtaMassPtTrue = new TH2F*[1];
1178  fhMCEtaPtTruePtRec = new TH2F*[1];
1179 
1180  fhMCPi0MassPtTrue[0] = new TH2F("hMCPi0MassPtTrue","Reconstructed Mass vs generated p_T of true #pi^{0} cluster pairs",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1181  fhMCPi0MassPtTrue[0]->SetXTitle("#it{p}_{T, generated} (GeV/#it{c})");
1182  fhMCPi0MassPtTrue[0]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1183  outputContainer->Add(fhMCPi0MassPtTrue[0]) ;
1184 
1185  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) ;
1186  fhMCPi0PtTruePtRec[0]->SetXTitle("#it{p}_{T, generated} (GeV/#it{c})");
1187  fhMCPi0PtTruePtRec[0]->SetYTitle("#it{p}_{T, reconstructed} (GeV/#it{c})");
1188  outputContainer->Add(fhMCPi0PtTruePtRec[0]) ;
1189 
1190  fhMCEtaMassPtTrue[0] = new TH2F("hMCEtaMassPtTrue","Reconstructed Mass vs generated p_T of true #eta cluster pairs",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1191  fhMCEtaMassPtTrue[0]->SetXTitle("#it{p}_{T, generated} (GeV/#it{c})");
1192  fhMCEtaMassPtTrue[0]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1193  outputContainer->Add(fhMCEtaMassPtTrue[0]) ;
1194 
1195  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) ;
1196  fhMCEtaPtTruePtRec[0]->SetXTitle("#it{p}_{T, generated} (GeV/#it{c})");
1197  fhMCEtaPtTruePtRec[0]->SetYTitle("#it{p}_{T, reconstructed} (GeV/#it{c})");
1198  outputContainer->Add(fhMCEtaPtTruePtRec[0]) ;
1199  }
1200  }
1201  }
1202 
1204  {
1206  {
1207  TString pairnamePHOS[] = {"(0-1)","(0-2)","(1-2)","(0-3)","(0-4)","(1-3)","(1-4)","(2-3)","(2-4)","(3-4)"};
1208  for(Int_t imod=0; imod<fNModules; imod++)
1209  {
1210  //Module dependent invariant mass
1211  snprintf(key, buffersize,"hReMod_%d",imod) ;
1212  snprintf(title, buffersize,"Real #it{M}_{#gamma#gamma} distr. for Module %d",imod) ;
1213  fhReMod[imod] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1214  fhReMod[imod]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1215  fhReMod[imod]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1216  outputContainer->Add(fhReMod[imod]) ;
1217  if(GetCalorimeter()==kPHOS)
1218  {
1219  snprintf(key, buffersize,"hReDiffPHOSMod_%d",imod) ;
1220  snprintf(title, buffersize,"Real pairs PHOS, clusters in different Modules: %s",(pairnamePHOS[imod]).Data()) ;
1221  fhReDiffPHOSMod[imod] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1222  fhReDiffPHOSMod[imod]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1223  fhReDiffPHOSMod[imod]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1224  outputContainer->Add(fhReDiffPHOSMod[imod]) ;
1225  }
1226  else
1227  {//EMCAL
1228  if(imod<fNModules/2)
1229  {
1230  snprintf(key, buffersize,"hReSameSectorEMCAL_%d",imod) ;
1231  snprintf(title, buffersize,"Real pairs EMCAL, clusters in same sector, SM(%d,%d)",imod*2,imod*2+1) ;
1232  fhReSameSectorEMCALMod[imod] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1233  fhReSameSectorEMCALMod[imod]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1234  fhReSameSectorEMCALMod[imod]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1235  outputContainer->Add(fhReSameSectorEMCALMod[imod]) ;
1236  }
1237  if(imod<fNModules-2)
1238  {
1239  snprintf(key, buffersize,"hReSameSideEMCAL_%d",imod) ;
1240  snprintf(title, buffersize,"Real pairs EMCAL, clusters in same side SM(%d,%d)",imod, imod+2) ;
1241  fhReSameSideEMCALMod[imod] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1242  fhReSameSideEMCALMod[imod]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1243  fhReSameSideEMCALMod[imod]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1244  outputContainer->Add(fhReSameSideEMCALMod[imod]) ;
1245  }
1246  }//EMCAL
1247 
1248  if(DoOwnMix())
1249  {
1250  snprintf(key, buffersize,"hMiMod_%d",imod) ;
1251  snprintf(title, buffersize,"Mixed #it{M}_{#gamma#gamma} distr. for Module %d",imod) ;
1252  fhMiMod[imod] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1253  fhMiMod[imod]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1254  fhMiMod[imod]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1255  outputContainer->Add(fhMiMod[imod]) ;
1256 
1257  if(GetCalorimeter()==kPHOS)
1258  {
1259  snprintf(key, buffersize,"hMiDiffPHOSMod_%d",imod) ;
1260  snprintf(title, buffersize,"Mixed pairs PHOS, clusters in different Modules: %s",(pairnamePHOS[imod]).Data()) ;
1261  fhMiDiffPHOSMod[imod] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1262  fhMiDiffPHOSMod[imod]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1263  fhMiDiffPHOSMod[imod]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1264  outputContainer->Add(fhMiDiffPHOSMod[imod]) ;
1265  }//PHOS
1266  else
1267  {//EMCAL
1268  if(imod<fNModules/2)
1269  {
1270  snprintf(key, buffersize,"hMiSameSectorEMCALMod_%d",imod) ;
1271  snprintf(title, buffersize,"Mixed pairs EMCAL, clusters in same sector, SM(%d,%d)",imod*2,imod*2+1) ;
1272  fhMiSameSectorEMCALMod[imod] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1273  fhMiSameSectorEMCALMod[imod]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1274  fhMiSameSectorEMCALMod[imod]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1275  outputContainer->Add(fhMiSameSectorEMCALMod[imod]) ;
1276  }
1277  if(imod<fNModules-2){
1278 
1279  snprintf(key, buffersize,"hMiSameSideEMCALMod_%d",imod) ;
1280  snprintf(title, buffersize,"Mixed pairs EMCAL, clusters in same side SM(%d,%d)",imod, imod+2) ;
1281  fhMiSameSideEMCALMod[imod] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1282  fhMiSameSideEMCALMod[imod]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1283  fhMiSameSideEMCALMod[imod]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1284  outputContainer->Add(fhMiSameSideEMCALMod[imod]) ;
1285  }
1286  } // EMCAL
1287  } // own mix
1288  } // loop combinations
1289  } // Not pair of detectors
1290  else
1291  {
1292  Int_t dcSameSM[6] = {12,13,14,15,16,17}; // Check eta order
1293  Int_t phSameSM[6] = {3, 3, 2, 2, 1, 1};
1294 
1295  Int_t dcDiffSM[8] = {12,13,14,15,16,17,0,0};
1296  Int_t phDiffSM[8] = {2, 2, 1, 1, 3, 3,0,0};
1297 
1298  for(Int_t icombi = 0; icombi < 8; icombi++)
1299  {
1300  snprintf(key, buffersize,"hReDiffSectorDCALPHOS_%d",icombi) ;
1301  snprintf(title, buffersize,"Real pairs DCAL-PHOS, clusters in different sector, SM(%d,%d)",dcDiffSM[icombi],phDiffSM[icombi]) ;
1302  fhReDiffSectorDCALPHOSMod[icombi] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1303  fhReDiffSectorDCALPHOSMod[icombi]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1304  fhReDiffSectorDCALPHOSMod[icombi]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1305  outputContainer->Add(fhReDiffSectorDCALPHOSMod[icombi]) ;
1306  if(DoOwnMix())
1307  {
1308  snprintf(key, buffersize,"hMiDiffSectorDCALPHOS_%d",icombi) ;
1309  snprintf(title, buffersize,"Mixed pairs DCAL-PHOS, clusters in different sector, SM(%d,%d)",dcDiffSM[icombi],phDiffSM[icombi]) ;
1310  fhMiDiffSectorDCALPHOSMod[icombi] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1311  fhMiDiffSectorDCALPHOSMod[icombi]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1312  fhMiDiffSectorDCALPHOSMod[icombi]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1313  outputContainer->Add(fhMiDiffSectorDCALPHOSMod[icombi]) ;
1314  }
1315 
1316  if(icombi > 5) continue ;
1317 
1318  snprintf(key, buffersize,"hReSameSectorDCALPHOS_%d",icombi) ;
1319  snprintf(title, buffersize,"Real pairs DCAL-PHOS, clusters in same sector, SM(%d,%d)",dcSameSM[icombi],phSameSM[icombi]) ;
1320  fhReSameSectorDCALPHOSMod[icombi] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1321  fhReSameSectorDCALPHOSMod[icombi]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1322  fhReSameSectorDCALPHOSMod[icombi]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1323  outputContainer->Add(fhReSameSectorDCALPHOSMod[icombi]) ;
1324  if(DoOwnMix())
1325  {
1326  snprintf(key, buffersize,"hMiSameSectorDCALPHOS_%d",icombi) ;
1327  snprintf(title, buffersize,"Mixed pairs DCAL-PHOS, clusters in same sector, SM(%d,%d)",dcSameSM[icombi],phSameSM[icombi]) ;
1328  fhMiSameSectorDCALPHOSMod[icombi] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1329  fhMiSameSectorDCALPHOSMod[icombi]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1330  fhMiSameSectorDCALPHOSMod[icombi]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1331  outputContainer->Add(fhMiSameSectorDCALPHOSMod[icombi]) ;
1332  }
1333  }
1334 
1335  }
1336  } // SM combinations
1337 
1338 // for(Int_t i = 0; i < outputContainer->GetEntries() ; i++){
1339 //
1340 // printf("Histogram %d, name: %s\n ",i, outputContainer->At(i)->GetName());
1341 //
1342 // }
1343 
1344  return outputContainer;
1345 }
1346 
1347 //___________________________________________________
1349 //___________________________________________________
1350 void AliAnaPi0::Print(const Option_t * /*opt*/) const
1351 {
1352  printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
1354 
1355  printf("Number of bins in Centrality: %d \n",GetNCentrBin()) ;
1356  printf("Number of bins in Z vert. pos: %d \n",GetNZvertBin()) ;
1357  printf("Number of bins in Reac. Plain: %d \n",GetNRPBin()) ;
1358  printf("Depth of event buffer: %d \n",GetNMaxEvMix()) ;
1359  printf("Pair in same Module: %d \n",fSameSM) ;
1360  printf("Cuts: \n") ;
1361  // printf("Z vertex position: -%2.3f < z < %2.3f \n",GetZvertexCut(),GetZvertexCut()) ; //It crashes here, why?
1362  printf("Number of modules: %d \n",fNModules) ;
1363  printf("Select pairs with their angle: %d, edep %d, min angle %2.3f, max angle %2.3f \n",fUseAngleCut, fUseAngleEDepCut, fAngleCut, fAngleMaxCut) ;
1364  printf("Asymmetry cuts: n = %d, \n",fNAsymCuts) ;
1365  printf("\tasymmetry < ");
1366  for(Int_t i = 0; i < fNAsymCuts; i++) printf("%2.2f ",fAsymCuts[i]);
1367  printf("\n");
1368 
1369  printf("PID selection bits: n = %d, \n",fNPIDBits) ;
1370  printf("\tPID bit = ");
1371  for(Int_t i = 0; i < fNPIDBits; i++) printf("%d ",fPIDBits[i]);
1372  printf("\n");
1373 
1374  if(fMultiCutAna)
1375  {
1376  printf("pT cuts: n = %d, \n",fNPtCuts) ;
1377  printf("\tpT > ");
1378  for(Int_t i = 0; i < fNPtCuts; i++) printf("%2.2f ",fPtCuts[i]);
1379  printf("GeV/c\n");
1380 
1381  printf("N cell in cluster cuts: n = %d, \n",fNCellNCuts) ;
1382  printf("\tnCell > ");
1383  for(Int_t i = 0; i < fNCellNCuts; i++) printf("%d ",fCellNCuts[i]);
1384  printf("\n");
1385 
1386  }
1387  printf("------------------------------------------------------\n") ;
1388 }
1389 
1390 //________________________________________
1392 //________________________________________
1394 {
1395  Double_t mesonY = -100 ;
1396  Double_t mesonE = -1 ;
1397  Double_t mesonPt = -1 ;
1398  Double_t mesonPhi = 100 ;
1399  Double_t mesonYeta= -1 ;
1400 
1401  Int_t pdg = 0 ;
1402  Int_t nprim = 0 ;
1403  Int_t nDaught = 0 ;
1404  Int_t iphot1 = 0 ;
1405  Int_t iphot2 = 0 ;
1406 
1407  Float_t cen = GetEventCentrality();
1408  Float_t ep = GetEventPlaneAngle();
1409 
1410  TParticle * primStack = 0;
1411  AliAODMCParticle * primAOD = 0;
1412 
1413  // Get the ESD MC particles container
1414  AliStack * stack = 0;
1415  if( GetReader()->ReadStack() )
1416  {
1417  stack = GetMCStack();
1418  if(!stack ) return;
1419  nprim = stack->GetNtrack();
1420  }
1421 
1422  // Get the AOD MC particles container
1423  TClonesArray * mcparticles = 0;
1424  if( GetReader()->ReadAODMCParticles() )
1425  {
1426  mcparticles = GetReader()->GetAODMCParticles();
1427  if( !mcparticles ) return;
1428  nprim = mcparticles->GetEntriesFast();
1429  }
1430 
1431  for(Int_t i=0 ; i < nprim; i++)
1432  {
1433  if(GetReader()->AcceptOnlyHIJINGLabels() && !GetReader()->IsHIJINGLabel(i)) continue ;
1434 
1435  if(GetReader()->ReadStack())
1436  {
1437  primStack = stack->Particle(i) ;
1438  if(!primStack)
1439  {
1440  AliWarning("ESD primaries pointer not available!!");
1441  continue;
1442  }
1443 
1444  // If too small skip
1445  if( primStack->Energy() < 0.4 ) continue;
1446 
1447  pdg = primStack->GetPdgCode();
1448  nDaught = primStack->GetNDaughters();
1449  iphot1 = primStack->GetDaughter(0) ;
1450  iphot2 = primStack->GetDaughter(1) ;
1451  if(primStack->Energy() == TMath::Abs(primStack->Pz())) continue ; //Protection against floating point exception
1452 
1453  //printf("i %d, %s %d %s %d \n",i, stack->Particle(i)->GetName(), stack->Particle(i)->GetPdgCode(),
1454  // prim->GetName(), prim->GetPdgCode());
1455 
1456  //Photon kinematics
1457  primStack->Momentum(fPi0Mom);
1458 
1459  mesonY = 0.5*TMath::Log((primStack->Energy()+primStack->Pz())/(primStack->Energy()-primStack->Pz())) ;
1460  }
1461  else
1462  {
1463  primAOD = (AliAODMCParticle *) mcparticles->At(i);
1464  if(!primAOD)
1465  {
1466  AliWarning("AOD primaries pointer not available!!");
1467  continue;
1468  }
1469 
1470  // If too small skip
1471  if( primAOD->E() < 0.4 ) continue;
1472 
1473  pdg = primAOD->GetPdgCode();
1474  nDaught = primAOD->GetNDaughters();
1475  iphot1 = primAOD->GetFirstDaughter() ;
1476  iphot2 = primAOD->GetLastDaughter() ;
1477 
1478  if(primAOD->E() == TMath::Abs(primAOD->Pz())) continue ; //Protection against floating point exception
1479 
1480  //Photon kinematics
1481  fPi0Mom.SetPxPyPzE(primAOD->Px(),primAOD->Py(),primAOD->Pz(),primAOD->E());
1482 
1483  mesonY = 0.5*TMath::Log((primAOD->E()+primAOD->Pz())/(primAOD->E()-primAOD->Pz())) ;
1484  }
1485 
1486  // Select only pi0 or eta
1487  if( pdg != 111 && pdg != 221) continue ;
1488 
1489  mesonPt = fPi0Mom.Pt () ;
1490  mesonE = fPi0Mom.E () ;
1491  mesonYeta= fPi0Mom.Eta() ;
1492  mesonPhi = fPi0Mom.Phi() ;
1493  if( mesonPhi < 0 ) mesonPhi+=TMath::TwoPi();
1494  mesonPhi *= TMath::RadToDeg();
1495 
1496  if(pdg == 111)
1497  {
1498  if(TMath::Abs(mesonY) < 1.0)
1499  {
1500  fhPrimPi0E ->Fill(mesonE , GetEventWeight()) ;
1501  fhPrimPi0Pt ->Fill(mesonPt, GetEventWeight()) ;
1502  fhPrimPi0Phi->Fill(mesonPt, mesonPhi, GetEventWeight()) ;
1503 
1504  fhPrimPi0YetaYcut->Fill(mesonPt, mesonYeta, GetEventWeight()) ;
1505 
1507  {
1508  fhPrimPi0PtCentrality->Fill(mesonPt, cen, GetEventWeight()) ;
1509  fhPrimPi0PtEventPlane->Fill(mesonPt, ep , GetEventWeight()) ;
1510  }
1511  }
1512 
1513  fhPrimPi0Y ->Fill(mesonPt, mesonY , GetEventWeight()) ;
1514  fhPrimPi0Yeta->Fill(mesonPt, mesonYeta, GetEventWeight()) ;
1515  }
1516  else if(pdg == 221)
1517  {
1518  if(TMath::Abs(mesonY) < 1.0)
1519  {
1520  fhPrimEtaE ->Fill(mesonE , GetEventWeight()) ;
1521  fhPrimEtaPt ->Fill(mesonPt, GetEventWeight()) ;
1522  fhPrimEtaPhi->Fill(mesonPt, mesonPhi, GetEventWeight()) ;
1523 
1524  fhPrimEtaYetaYcut->Fill(mesonPt, mesonYeta, GetEventWeight()) ;
1525 
1527  {
1528  fhPrimEtaPtCentrality->Fill(mesonPt, cen, GetEventWeight()) ;
1529  fhPrimEtaPtEventPlane->Fill(mesonPt, ep , GetEventWeight()) ;
1530  }
1531  }
1532 
1533  fhPrimEtaY ->Fill(mesonPt, mesonY , GetEventWeight()) ;
1534  fhPrimEtaYeta->Fill(mesonPt, mesonYeta, GetEventWeight()) ;
1535  }
1536 
1537  //Origin of meson
1538  if(fFillOriginHisto && TMath::Abs(mesonY) < 0.7)
1539  {
1540  Int_t momindex = -1;
1541  Int_t mompdg = -1;
1542  Int_t momstatus = -1;
1543  Int_t status = -1;
1544  Float_t momR = 0;
1545  if(GetReader()->ReadStack()) momindex = primStack->GetFirstMother();
1546  if(GetReader()->ReadAODMCParticles()) momindex = primAOD ->GetMother();
1547 
1548  if(momindex >= 0 && momindex < nprim)
1549  {
1550  if(GetReader()->ReadStack())
1551  {
1552  status = primStack->GetStatusCode();
1553  TParticle* mother = stack->Particle(momindex);
1554  mompdg = TMath::Abs(mother->GetPdgCode());
1555  momstatus = mother->GetStatusCode();
1556  momR = mother->R();
1557  }
1558 
1559  if(GetReader()->ReadAODMCParticles())
1560  {
1561  status = primAOD->GetStatus();
1562  AliAODMCParticle* mother = (AliAODMCParticle*) mcparticles->At(momindex);
1563  mompdg = TMath::Abs(mother->GetPdgCode());
1564  momstatus = mother->GetStatus();
1565  momR = TMath::Sqrt(mother->Xv()*mother->Xv()+mother->Yv()*mother->Yv());
1566  }
1567 
1568  if(pdg == 111)
1569  {
1570  fhPrimPi0ProdVertex->Fill(mesonPt, momR , GetEventWeight());
1571  fhPrimPi0PtStatus ->Fill(mesonPt, status, GetEventWeight());
1572 
1573 
1574  if (momstatus == 21) fhPrimPi0PtOrigin->Fill(mesonPt, 0.5, GetEventWeight());//parton
1575  else if(mompdg < 22 ) fhPrimPi0PtOrigin->Fill(mesonPt, 1.5, GetEventWeight());//quark
1576  else if(mompdg > 2100 && mompdg < 2210)
1577  fhPrimPi0PtOrigin->Fill(mesonPt, 2.5, GetEventWeight());// resonances
1578  else if(mompdg == 221) fhPrimPi0PtOrigin->Fill(mesonPt, 8.5, GetEventWeight());//eta
1579  else if(mompdg == 331) fhPrimPi0PtOrigin->Fill(mesonPt, 9.5, GetEventWeight());//eta prime
1580  else if(mompdg == 213) fhPrimPi0PtOrigin->Fill(mesonPt, 4.5, GetEventWeight());//rho
1581  else if(mompdg == 223) fhPrimPi0PtOrigin->Fill(mesonPt, 5.5, GetEventWeight());//omega
1582  else if(mompdg >= 310 && mompdg <= 323)
1583  fhPrimPi0PtOrigin->Fill(mesonPt, 6.5, GetEventWeight());//k0S, k+-,k*
1584  else if(mompdg == 130) fhPrimPi0PtOrigin->Fill(mesonPt, 6.5, GetEventWeight());//k0L
1585  else if(momstatus == 11 || momstatus == 12 )
1586  fhPrimPi0PtOrigin->Fill(mesonPt, 3.5, GetEventWeight());//resonances
1587  else fhPrimPi0PtOrigin->Fill(mesonPt, 7.5, GetEventWeight());//other?
1588 
1589  //printf("Prim Pi0: index %d, pt %2.2f Prod vertex %3.3f, origin pdg %d, origin status %d, origin UI %d\n",
1590  // momindex, mesonPt,mother->R(),mompdg,momstatus,mother->GetUniqueID());
1591 
1592  if(status!=11)
1593  {
1594  if (momstatus == 21) fhPrimNotResonancePi0PtOrigin->Fill(mesonPt, 0.5, GetEventWeight());//parton
1595  else if(mompdg < 22 ) fhPrimNotResonancePi0PtOrigin->Fill(mesonPt, 1.5, GetEventWeight());//quark
1596  else if(mompdg > 2100 && mompdg < 2210)
1597  fhPrimNotResonancePi0PtOrigin->Fill(mesonPt, 2.5, GetEventWeight());// resonances
1598  else if(mompdg == 221) fhPrimNotResonancePi0PtOrigin->Fill(mesonPt, 8.5, GetEventWeight());//eta
1599  else if(mompdg == 331) fhPrimNotResonancePi0PtOrigin->Fill(mesonPt, 9.5, GetEventWeight());//eta prime
1600  else if(mompdg == 213) fhPrimNotResonancePi0PtOrigin->Fill(mesonPt, 4.5, GetEventWeight());//rho
1601  else if(mompdg == 223) fhPrimNotResonancePi0PtOrigin->Fill(mesonPt, 5.5, GetEventWeight());//omega
1602  else if(mompdg >= 310 && mompdg <= 323)
1603  fhPrimNotResonancePi0PtOrigin->Fill(mesonPt, 6.5, GetEventWeight());//k0S, k+-,k*
1604  else if(mompdg == 130) fhPrimNotResonancePi0PtOrigin->Fill(mesonPt, 6.5, GetEventWeight());//k0L
1605  else if(momstatus == 11 || momstatus == 12 )
1606  fhPrimNotResonancePi0PtOrigin->Fill(mesonPt, 3.5, GetEventWeight());//resonances
1607  else fhPrimNotResonancePi0PtOrigin->Fill(mesonPt, 7.5, GetEventWeight());//other?
1608  }
1609 
1610  }//pi0
1611  else
1612  {
1613  if (momstatus == 21 ) fhPrimEtaPtOrigin->Fill(mesonPt, 0.5, GetEventWeight());//parton
1614  else if(mompdg < 22 ) fhPrimEtaPtOrigin->Fill(mesonPt, 1.5, GetEventWeight());//quark
1615  else if(mompdg > 2100 && mompdg < 2210)
1616  fhPrimEtaPtOrigin->Fill(mesonPt, 2.5, GetEventWeight());//qq resonances
1617  else if(mompdg == 331) fhPrimEtaPtOrigin->Fill(mesonPt, 5.5, GetEventWeight());//eta prime
1618  else if(momstatus == 11 || momstatus == 12 )
1619  fhPrimEtaPtOrigin->Fill(mesonPt, 3.5, GetEventWeight());//resonances
1620  else fhPrimEtaPtOrigin->Fill(mesonPt, 4.5, GetEventWeight());//stable, conversions?
1621  //printf("Other Meson pdg %d, Mother %s, pdg %d, status %d\n",pdg, TDatabasePDG::Instance()->GetParticle(mompdg)->GetName(),mompdg, momstatus );
1622 
1623  fhPrimEtaProdVertex->Fill(mesonPt, momR, GetEventWeight());
1624  }
1625  } // pi0 has mother
1626  }
1627 
1628  //Check if both photons hit Calorimeter
1629  if(nDaught != 2 ) continue; //Only interested in 2 gamma decay
1630 
1631  if(iphot1 < 0 || iphot1 >= nprim || iphot2 < 0 || iphot2 >= nprim) continue ;
1632 
1633  Int_t pdg1 = 0;
1634  Int_t pdg2 = 0;
1635  Bool_t inacceptance1 = kTRUE;
1636  Bool_t inacceptance2 = kTRUE;
1637 
1638  if(GetReader()->ReadStack())
1639  {
1640  TParticle * phot1 = stack->Particle(iphot1) ;
1641  TParticle * phot2 = stack->Particle(iphot2) ;
1642 
1643  if(!phot1 || !phot2) continue ;
1644 
1645  pdg1 = phot1->GetPdgCode();
1646  pdg2 = phot2->GetPdgCode();
1647 
1648  phot1->Momentum(fPhotonMom1);
1649  phot2->Momentum(fPhotonMom2);
1650 
1651  // Check if photons hit the Calorimeter acceptance
1653  {
1654  if( !GetCaloUtils()->IsMCParticleInCalorimeterAcceptance( GetCalorimeter(), phot1 )) inacceptance1 = kFALSE ;
1655  if( !GetCaloUtils()->IsMCParticleInCalorimeterAcceptance( GetCalorimeter(), phot2 )) inacceptance2 = kFALSE ;
1656  }
1657  }
1658 
1659  if(GetReader()->ReadAODMCParticles())
1660  {
1661  AliAODMCParticle * phot1 = (AliAODMCParticle *) mcparticles->At(iphot1) ;
1662  AliAODMCParticle * phot2 = (AliAODMCParticle *) mcparticles->At(iphot2) ;
1663 
1664  if(!phot1 || !phot2) continue ;
1665 
1666  pdg1 = phot1->GetPdgCode();
1667  pdg2 = phot2->GetPdgCode();
1668 
1669  fPhotonMom1.SetPxPyPzE(phot1->Px(),phot1->Py(),phot1->Pz(),phot1->E());
1670  fPhotonMom2.SetPxPyPzE(phot2->Px(),phot2->Py(),phot2->Pz(),phot2->E());
1671 
1672  // Check if photons hit the Calorimeter acceptance
1674  {
1675  if( !GetCaloUtils()->IsMCParticleInCalorimeterAcceptance( GetCalorimeter(), phot1 )) inacceptance1 = kFALSE ;
1676  if( !GetCaloUtils()->IsMCParticleInCalorimeterAcceptance( GetCalorimeter(), phot2 )) inacceptance2 = kFALSE ;
1677  }
1678  }
1679 
1680  if( pdg1 != 22 || pdg2 !=22) continue ;
1681 
1682  // Check if photons hit desired acceptance in the fidutial borders fixed in the analysis
1683  if(IsFiducialCutOn())
1684  {
1685  if( inacceptance1 && !GetFiducialCut()->IsInFiducialCut(fPhotonMom1.Eta(), fPhotonMom1.Phi(), GetCalorimeter()) ) inacceptance1 = kFALSE ;
1686  if( inacceptance2 && !GetFiducialCut()->IsInFiducialCut(fPhotonMom2.Eta(), fPhotonMom2.Phi(), GetCalorimeter()) ) inacceptance2 = kFALSE ;
1687  }
1688 
1690 
1691  if(GetCalorimeter()==kEMCAL && inacceptance1 && inacceptance2 && fCheckAccInSector)
1692  {
1693  Int_t absID1=0;
1694  Int_t absID2=0;
1695 
1696  Float_t photonPhi1 = fPhotonMom1.Phi();
1697  Float_t photonPhi2 = fPhotonMom2.Phi();
1698 
1699  if(photonPhi1 < 0) photonPhi1+=TMath::TwoPi();
1700  if(photonPhi2 < 0) photonPhi2+=TMath::TwoPi();
1701 
1702  GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(fPhotonMom1.Eta(),photonPhi1,absID1);
1703  GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(fPhotonMom2.Eta(),photonPhi2,absID2);
1704 
1705  Int_t sm1 = GetEMCALGeometry()->GetSuperModuleNumber(absID1);
1706  Int_t sm2 = GetEMCALGeometry()->GetSuperModuleNumber(absID2);
1707 
1708  Int_t j=0;
1709  Bool_t sameSector = kFALSE;
1710  for(Int_t isector = 0; isector < fNModules/2; isector++)
1711  {
1712  j=2*isector;
1713  if((sm1==j && sm2==j+1) || (sm1==j+1 && sm2==j)) sameSector = kTRUE;
1714  }
1715 
1716  if(sm1!=sm2 && !sameSector)
1717  {
1718  inacceptance1 = kFALSE;
1719  inacceptance2 = kFALSE;
1720  }
1721  //if(sm1!=sm2)printf("sm1 %d, sm2 %d, same sector %d, in acceptance %d\n",sm1,sm2,sameSector,inacceptance);
1722  // if(GetEMCALGeometry()->Impact(phot1) && GetEMCALGeometry()->Impact(phot2))
1723  // inacceptance = kTRUE;
1724  }
1725 
1726  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",
1728  mesonPt,mesonYeta,mesonPhi,
1729  fPhotonMom1.Pt(),fPhotonMom1.Eta(),fPhotonMom1.Phi()*TMath::RadToDeg(),
1730  fPhotonMom2.Pt(),fPhotonMom2.Eta(),fPhotonMom2.Phi()*TMath::RadToDeg(),
1731  inacceptance1, inacceptance2));
1732 
1733  if(inacceptance1 && inacceptance2)
1734  {
1735  Float_t asym = TMath::Abs((fPhotonMom1.E()-fPhotonMom2.E()) / (fPhotonMom1.E()+fPhotonMom2.E()));
1736  Double_t angle = fPhotonMom1.Angle(fPhotonMom2.Vect());
1737 
1738  Bool_t cutAngle = kFALSE;
1739  if(fUseAngleCut && (angle < fAngleCut || angle > fAngleMaxCut)) cutAngle = kTRUE;
1740 
1741  AliDebug(2,Form("\t ACCEPTED pdg %d: pt %2.2f, phi %2.2f, eta %2.2f",pdg,mesonPt,mesonPhi,mesonYeta));
1742 
1743  if(pdg==111)
1744  {
1745  fhPrimPi0AccE ->Fill(mesonE , GetEventWeight()) ;
1746  fhPrimPi0AccPt ->Fill(mesonPt, GetEventWeight()) ;
1747  fhPrimPi0AccPhi ->Fill(mesonPt, mesonPhi , GetEventWeight()) ;
1748  fhPrimPi0AccY ->Fill(mesonPt, mesonY , GetEventWeight()) ;
1749  fhPrimPi0AccYeta->Fill(mesonPt, mesonYeta, GetEventWeight()) ;
1750 
1752  {
1753  fhPrimPi0AccPtCentrality->Fill(mesonPt, cen, GetEventWeight()) ;
1754  fhPrimPi0AccPtEventPlane->Fill(mesonPt, ep , GetEventWeight()) ;
1755  }
1756 
1757  if(fFillAngleHisto)
1758  {
1759  fhPrimPi0OpeningAngle ->Fill(mesonPt, angle, GetEventWeight());
1760  if(mesonPt > 5)fhPrimPi0OpeningAngleAsym->Fill(asym, angle, GetEventWeight());
1761  fhPrimPi0CosOpeningAngle ->Fill(mesonPt, TMath::Cos(angle), GetEventWeight());
1762  }
1763 
1764  if(fPhotonMom1.Pt() > GetMinPt() && fPhotonMom2.Pt() > GetMinPt() && !cutAngle)
1765  {
1766  fhPrimPi0AccPtPhotonCuts->Fill(mesonPt, GetEventWeight()) ;
1767 
1768  if(fFillAngleHisto)
1769  fhPrimPi0OpeningAnglePhotonCuts->Fill(mesonPt, angle, GetEventWeight());
1770  }
1771  }
1772  else if(pdg==221)
1773  {
1774  fhPrimEtaAccE ->Fill(mesonE , GetEventWeight()) ;
1775  fhPrimEtaAccPt ->Fill(mesonPt, GetEventWeight()) ;
1776  fhPrimEtaAccPhi ->Fill(mesonPt, mesonPhi , GetEventWeight()) ;
1777  fhPrimEtaAccY ->Fill(mesonPt, mesonY , GetEventWeight()) ;
1778  fhPrimEtaAccYeta->Fill(mesonPt, mesonYeta, GetEventWeight()) ;
1779 
1781  {
1782  fhPrimEtaAccPtCentrality->Fill(mesonPt, cen, GetEventWeight()) ;
1783  fhPrimEtaAccPtEventPlane->Fill(mesonPt, ep , GetEventWeight()) ;
1784  }
1785 
1786  if(fFillAngleHisto)
1787  {
1788  fhPrimEtaOpeningAngle ->Fill(mesonPt, angle, GetEventWeight());
1789  if(mesonPt > 5)fhPrimEtaOpeningAngleAsym->Fill(asym, angle, GetEventWeight());
1790  fhPrimEtaCosOpeningAngle ->Fill(mesonPt, TMath::Cos(angle), GetEventWeight());
1791 
1792  if(fPhotonMom1.Pt() > GetMinPt() && fPhotonMom2.Pt() > GetMinPt() && !cutAngle)
1793  {
1794  if(fUseAngleCut && angle > fAngleCut && angle < fAngleMaxCut)
1795  fhPrimEtaAccPtPhotonCuts->Fill(mesonPt, GetEventWeight()) ;
1796 
1797  if(fFillAngleHisto)
1798  fhPrimEtaOpeningAnglePhotonCuts->Fill(mesonPt, angle, GetEventWeight());
1799  }
1800  }
1801  }
1802  } // Accepted
1803  } // loop on primaries
1804 }
1805 
1806 //________________________________________________
1808 //________________________________________________
1810 {
1811  // Get pTArm and AlphaArm
1812  Float_t momentumSquaredMother = fPi0Mom.P()*fPi0Mom.P();
1813  Float_t momentumDaughter1AlongMother = 0.;
1814  Float_t momentumDaughter2AlongMother = 0.;
1815 
1816  if (momentumSquaredMother > 0.)
1817  {
1818  momentumDaughter1AlongMother = (fPhotonMom1.Px()*fPi0Mom.Px() + fPhotonMom1.Py()*fPi0Mom.Py()+ fPhotonMom1.Pz()*fPi0Mom.Pz()) / sqrt(momentumSquaredMother);
1819  momentumDaughter2AlongMother = (fPhotonMom2.Px()*fPi0Mom.Px() + fPhotonMom2.Py()*fPi0Mom.Py()+ fPhotonMom2.Pz()*fPi0Mom.Pz()) / sqrt(momentumSquaredMother);
1820  }
1821 
1822  Float_t momentumSquaredDaughter1 = fPhotonMom1.P()*fPhotonMom1.P();
1823  Float_t ptArmSquared = momentumSquaredDaughter1 - momentumDaughter1AlongMother*momentumDaughter1AlongMother;
1824 
1825  Float_t pTArm = 0.;
1826  if (ptArmSquared > 0.)
1827  pTArm = sqrt(ptArmSquared);
1828 
1829  Float_t alphaArm = 0.;
1830  if(momentumDaughter1AlongMother + momentumDaughter2AlongMother > 0)
1831  alphaArm = (momentumDaughter1AlongMother -momentumDaughter2AlongMother) / (momentumDaughter1AlongMother + momentumDaughter2AlongMother);
1832 
1834  fPhotonMom1Boost.Boost(-fPi0Mom.BoostVector());
1835  Float_t cosThStar=TMath::Cos(fPhotonMom1Boost.Vect().Angle(fPi0Mom.Vect()));
1836 
1837  Float_t en = fPi0Mom.Energy();
1838  Int_t ebin = -1;
1839  if(en > 8 && en <= 12) ebin = 0;
1840  if(en > 12 && en <= 16) ebin = 1;
1841  if(en > 16 && en <= 20) ebin = 2;
1842  if(en > 20) ebin = 3;
1843  if(ebin < 0 || ebin > 3) return ;
1844 
1845  if(pdg==111)
1846  {
1847  fhCosThStarPrimPi0->Fill(en, cosThStar, GetEventWeight());
1848  fhArmPrimPi0[ebin]->Fill(alphaArm, pTArm , GetEventWeight());
1849  }
1850  else
1851  {
1852  fhCosThStarPrimEta->Fill(en, cosThStar, GetEventWeight());
1853  fhArmPrimEta[ebin]->Fill(alphaArm,pTArm , GetEventWeight());
1854  }
1855 
1856  // if(GetDebug() > 2 )
1857  // {
1858  // Float_t asym = 0;
1859  // if(fPhotonMom1.E()+fPhotonMom2.E() > 0) asym = TMath::Abs(fPhotonMom1.E()-fPhotonMom2.E())/(fPhotonMom1.E()+fPhotonMom2.E());
1860  //
1861  // printf("AliAnaPi0::FillArmenterosThetaStar() - E %f, alphaArm %f, pTArm %f, cos(theta*) %f, asymmetry %1.3f\n",
1862  // en,alphaArm,pTArm,cosThStar,asym);
1863  // }
1864 }
1865 
1866 //_______________________________________________________________________________________
1871 //_______________________________________________________________________________________
1872 void AliAnaPi0::FillMCVersusRecDataHistograms(Int_t index1, Int_t index2,
1873  Float_t pt1, Float_t pt2,
1874  Int_t ncell1, Int_t ncell2,
1875  Double_t mass, Double_t pt, Double_t asym,
1876  Double_t deta, Double_t dphi)
1877 {
1878  Int_t ancPDG = 0;
1879  Int_t ancStatus = 0;
1880  Int_t ancLabel = GetMCAnalysisUtils()->CheckCommonAncestor(index1, index2,
1881  GetReader(), ancPDG, ancStatus,fPi0Mom, fProdVertex);
1882 
1883  Int_t momindex = -1;
1884  Int_t mompdg = -1;
1885  Int_t momstatus = -1;
1886  Int_t status = -1;
1887  Float_t prodR = -1;
1888  Int_t mcIndex = -1;
1889 
1890  if(ancLabel > -1)
1891  {
1892  AliDebug(1,Form("Common ancestor label %d, pdg %d, name %s, status %d",
1893  ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus));
1894 
1895  if(ancPDG==22)
1896  {//gamma
1897  mcIndex = 0;
1898  }
1899  else if(TMath::Abs(ancPDG)==11)
1900  {//e
1901  mcIndex = 1;
1902  }
1903  else if(ancPDG==111)
1904  {//Pi0
1905  mcIndex = 2;
1906  if(fMultiCutAnaSim)
1907  {
1908  for(Int_t ipt=0; ipt<fNPtCuts; ipt++)
1909  {
1910  for(Int_t icell=0; icell<fNCellNCuts; icell++)
1911  {
1912  for(Int_t iasym=0; iasym<fNAsymCuts; iasym++)
1913  {
1914  Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
1915  if(pt1 > fPtCuts[ipt] && pt2 > fPtCuts[ipt] &&
1916  asym < fAsymCuts[iasym] &&
1917  ncell1 >= fCellNCuts[icell] && ncell2 >= fCellNCuts[icell])
1918  {
1919  fhMCPi0MassPtRec [index]->Fill(pt, mass, GetEventWeight());
1920  fhMCPi0MassPtTrue[index]->Fill(fPi0Mom.Pt(),mass, GetEventWeight());
1921  if(mass < 0.17 && mass > 0.1)
1922  fhMCPi0PtTruePtRec[index]->Fill(fPi0Mom.Pt(), pt, GetEventWeight());
1923  }//pass the different cuts
1924  }// pid bit cut loop
1925  }// icell loop
1926  }// pt cut loop
1927  }// Multi cut ana sim
1928  else
1929  {
1930  fhMCPi0MassPtTrue[0]->Fill(fPi0Mom.Pt(), mass, GetEventWeight());
1931 
1932  if(mass < 0.17 && mass > 0.1)
1933  {
1934  fhMCPi0PtTruePtRec[0]->Fill(fPi0Mom.Pt(), pt, GetEventWeight());
1935 
1936  //Int_t uniqueId = -1;
1937  if(GetReader()->ReadStack())
1938  {
1939  TParticle* ancestor = GetMCStack()->Particle(ancLabel);
1940  status = ancestor->GetStatusCode();
1941  momindex = ancestor->GetFirstMother();
1942  if(momindex < 0) return;
1943  TParticle* mother = GetMCStack()->Particle(momindex);
1944  mompdg = TMath::Abs(mother->GetPdgCode());
1945  momstatus = mother->GetStatusCode();
1946  prodR = mother->R();
1947  //uniqueId = mother->GetUniqueID();
1948  }
1949  else
1950  {
1951  TClonesArray * mcparticles = GetReader()->GetAODMCParticles();
1952  AliAODMCParticle* ancestor = (AliAODMCParticle *) mcparticles->At(ancLabel);
1953  status = ancestor->GetStatus();
1954  momindex = ancestor->GetMother();
1955  if(momindex < 0) return;
1956  AliAODMCParticle* mother = (AliAODMCParticle *) mcparticles->At(momindex);
1957  mompdg = TMath::Abs(mother->GetPdgCode());
1958  momstatus = mother->GetStatus();
1959  prodR = TMath::Sqrt(mother->Xv()*mother->Xv()+mother->Yv()*mother->Yv());
1960  //uniqueId = mother->GetUniqueID();
1961  }
1962 
1963  // printf("Reco Pi0: pt %2.2f Prod vertex %3.3f, origin pdg %d, origin status %d, origin UI %d\n",
1964  // pt,prodR,mompdg,momstatus,uniqueId);
1965 
1966  fhMCPi0ProdVertex->Fill(pt, prodR , GetEventWeight());
1967  fhMCPi0PtStatus ->Fill(pt, status, GetEventWeight());
1968 
1969  if (momstatus == 21) fhMCPi0PtOrigin->Fill(pt, 0.5, GetEventWeight());//parton
1970  else if(mompdg < 22 ) fhMCPi0PtOrigin->Fill(pt, 1.5, GetEventWeight());//quark
1971  else if(mompdg > 2100 && mompdg < 2210)
1972  fhMCPi0PtOrigin->Fill(pt, 2.5, GetEventWeight());// resonances
1973  else if(mompdg == 221) fhMCPi0PtOrigin->Fill(pt, 8.5, GetEventWeight());//eta
1974  else if(mompdg == 331) fhMCPi0PtOrigin->Fill(pt, 9.5, GetEventWeight());//eta prime
1975  else if(mompdg == 213) fhMCPi0PtOrigin->Fill(pt, 4.5, GetEventWeight());//rho
1976  else if(mompdg == 223) fhMCPi0PtOrigin->Fill(pt, 5.5, GetEventWeight());//omega
1977  else if(mompdg >= 310 && mompdg <= 323)
1978  fhMCPi0PtOrigin->Fill(pt, 6.5, GetEventWeight());//k0S, k+-,k*
1979  else if(mompdg == 130) fhMCPi0PtOrigin->Fill(pt, 6.5, GetEventWeight());//k0L
1980  else if(momstatus == 11 || momstatus == 12 )
1981  fhMCPi0PtOrigin->Fill(pt, 3.5, GetEventWeight());//resonances
1982  else fhMCPi0PtOrigin->Fill(pt, 7.5, GetEventWeight());//other?
1983 
1984  if(status!=11)
1985  {
1986  if (momstatus == 21) fhMCNotResonancePi0PtOrigin->Fill(pt, 0.5, GetEventWeight());//parton
1987  else if(mompdg < 22 ) fhMCNotResonancePi0PtOrigin->Fill(pt, 1.5, GetEventWeight());//quark
1988  else if(mompdg > 2100 && mompdg < 2210)
1989  fhMCNotResonancePi0PtOrigin->Fill(pt, 2.5, GetEventWeight());// resonances
1990  else if(mompdg == 221) fhMCNotResonancePi0PtOrigin->Fill(pt, 8.5, GetEventWeight());//eta
1991  else if(mompdg == 331) fhMCNotResonancePi0PtOrigin->Fill(pt, 9.5, GetEventWeight());//eta prime
1992  else if(mompdg == 213) fhMCNotResonancePi0PtOrigin->Fill(pt, 4.5, GetEventWeight());//rho
1993  else if(mompdg == 223) fhMCNotResonancePi0PtOrigin->Fill(pt, 5.5, GetEventWeight());//omega
1994  else if(mompdg >= 310 && mompdg <= 323)
1995  fhMCNotResonancePi0PtOrigin->Fill(pt, 6.5, GetEventWeight());//k0S, k+-,k*
1996  else if(mompdg == 130) fhMCNotResonancePi0PtOrigin->Fill(pt, 6.5, GetEventWeight());//k0L
1997  else if(momstatus == 11 || momstatus == 12 )
1998  fhMCNotResonancePi0PtOrigin->Fill(pt, 3.5, GetEventWeight());//resonances
1999  else fhMCNotResonancePi0PtOrigin->Fill(pt, 7.5, GetEventWeight());//other?
2000  }
2001  }//pi0 mass region
2002  }
2003  }
2004  else if(ancPDG==221)
2005  {//Eta
2006  mcIndex = 3;
2007  if(fMultiCutAnaSim)
2008  {
2009  for(Int_t ipt=0; ipt<fNPtCuts; ipt++)
2010  {
2011  for(Int_t icell=0; icell<fNCellNCuts; icell++)
2012  {
2013  for(Int_t iasym=0; iasym<fNAsymCuts; iasym++)
2014  {
2015  Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
2016  if(pt1 > fPtCuts[ipt] && pt2 > fPtCuts[ipt] &&
2017  asym < fAsymCuts[iasym] &&
2018  ncell1 >= fCellNCuts[icell] && ncell2 >= fCellNCuts[icell])
2019  {
2020  fhMCEtaMassPtRec [index]->Fill(pt , mass, GetEventWeight());
2021  fhMCEtaMassPtTrue[index]->Fill(fPi0Mom.Pt(), mass, GetEventWeight());
2022  if(mass < 0.65 && mass > 0.45)
2023  fhMCEtaPtTruePtRec[index]->Fill(fPi0Mom.Pt(), pt, GetEventWeight());
2024  }//pass the different cuts
2025  }// pid bit cut loop
2026  }// icell loop
2027  }// pt cut loop
2028  } //Multi cut ana sim
2029  else
2030  {
2031  fhMCEtaMassPtTrue[0]->Fill(fPi0Mom.Pt(), mass, GetEventWeight());
2032  if(mass < 0.65 && mass > 0.45) fhMCEtaPtTruePtRec[0]->Fill(fPi0Mom.Pt(), pt, GetEventWeight());
2033 
2034  if(GetReader()->ReadStack())
2035  {
2036  TParticle* ancestor = GetMCStack()->Particle(ancLabel);
2037  momindex = ancestor->GetFirstMother();
2038  if(momindex < 0) return;
2039  TParticle* mother = GetMCStack()->Particle(momindex);
2040  mompdg = TMath::Abs(mother->GetPdgCode());
2041  momstatus = mother->GetStatusCode();
2042  prodR = mother->R();
2043  }
2044  else
2045  {
2046  TClonesArray * mcparticles = GetReader()->GetAODMCParticles();
2047  AliAODMCParticle* ancestor = (AliAODMCParticle *) mcparticles->At(ancLabel);
2048  momindex = ancestor->GetMother();
2049  if(momindex < 0) return;
2050  AliAODMCParticle* mother = (AliAODMCParticle *) mcparticles->At(momindex);
2051  mompdg = TMath::Abs(mother->GetPdgCode());
2052  momstatus = mother->GetStatus();
2053  prodR = TMath::Sqrt(mother->Xv()*mother->Xv()+mother->Yv()*mother->Yv());
2054  }
2055 
2056  fhMCEtaProdVertex->Fill(pt, prodR, GetEventWeight());
2057 
2058  if (momstatus == 21 ) fhMCEtaPtOrigin->Fill(pt, 0.5, GetEventWeight());//parton
2059  else if(mompdg < 22 ) fhMCEtaPtOrigin->Fill(pt, 1.5, GetEventWeight());//quark
2060  else if(mompdg > 2100 && mompdg < 2210)
2061  fhMCEtaPtOrigin->Fill(pt, 2.5, GetEventWeight());//qq resonances
2062  else if(mompdg == 331) fhMCEtaPtOrigin->Fill(pt, 5.5, GetEventWeight());//eta prime
2063  else if(momstatus == 11 || momstatus == 12 )
2064  fhMCEtaPtOrigin->Fill(pt, 3.5, GetEventWeight());//resonances
2065  else fhMCEtaPtOrigin->Fill(pt, 4.5, GetEventWeight());//stable, conversions?
2066  //printf("Other Meson pdg %d, Mother %s, pdg %d, status %d\n",pdg, TDatabasePDG::Instance()->GetParticle(mompdg)->GetName(),mompdg, momstatus );
2067 
2068  }// eta mass region
2069  }
2070  else if(ancPDG==-2212)
2071  {//AProton
2072  mcIndex = 4;
2073  }
2074  else if(ancPDG==-2112)
2075  {//ANeutron
2076  mcIndex = 5;
2077  }
2078  else if(TMath::Abs(ancPDG)==13)
2079  {//muons
2080  mcIndex = 6;
2081  }
2082  else if (TMath::Abs(ancPDG) > 100 && ancLabel > 7)
2083  {
2084  if(ancStatus==1)
2085  {//Stable particles, converted? not decayed resonances
2086  mcIndex = 6;
2087  }
2088  else
2089  {//resonances and other decays, more hadron conversions?
2090  mcIndex = 7;
2091  }
2092  }
2093  else
2094  {//Partons, colliding protons, strings, intermediate corrections
2095  if(ancStatus==11 || ancStatus==12)
2096  {//String fragmentation
2097  mcIndex = 8;
2098  }
2099  else if (ancStatus==21){
2100  if(ancLabel < 2)
2101  {//Colliding protons
2102  mcIndex = 11;
2103  }//colliding protons
2104  else if(ancLabel < 6)
2105  {//partonic initial states interactions
2106  mcIndex = 9;
2107  }
2108  else if(ancLabel < 8)
2109  {//Final state partons radiations?
2110  mcIndex = 10;
2111  }
2112  // else {
2113  // printf("AliAnaPi0::FillMCVersusRecDataHistograms() - Check ** Common ancestor label %d, pdg %d, name %s, status %d; \n",
2114  // ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus);
2115  // }
2116  }//status 21
2117  //else {
2118  // printf("AliAnaPi0::FillMCVersusRecDataHistograms() - Check *** Common ancestor label %d, pdg %d, name %s, status %d; \n",
2119  // ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus);
2120  // }
2121  }
2122  }//ancLabel > -1
2123  else
2124  { //ancLabel <= -1
2125  //printf("Not related at all label = %d\n",ancLabel);
2126  AliDebug(1,"Common ancestor not found");
2127 
2128  mcIndex = 12;
2129  }
2130 
2131  if(mcIndex >= 0 && mcIndex < 13)
2132  {
2133  fhMCOrgMass [mcIndex]->Fill(pt, mass, GetEventWeight());
2134  fhMCOrgAsym [mcIndex]->Fill(pt, asym, GetEventWeight());
2135  fhMCOrgDeltaEta[mcIndex]->Fill(pt, deta, GetEventWeight());
2136  fhMCOrgDeltaPhi[mcIndex]->Fill(pt, dphi, GetEventWeight());
2137  }
2138 }
2139 
2140 //__________________________________________
2144 //__________________________________________
2146 {
2147  // In case of simulated data, fill acceptance histograms
2148  if(IsDataMC())
2149  {
2151 
2152  if(fFillOnlyMCAcceptanceHisto) return;
2153  }
2154 
2155 //if (GetReader()->GetEventNumber()%10000 == 0)
2156 // printf("--- Event %d ---\n",GetReader()->GetEventNumber());
2157 
2158  if(!GetInputAODBranch())
2159  {
2160  AliFatal(Form("No input aod photons in AOD with name branch < %s >, STOP",GetInputAODName().Data()));
2161  return;
2162  }
2163 
2164  //
2165  // Init some variables
2166  //
2167 
2168  // Analysis within the same detector:
2169  TClonesArray* secondLoopInputData = GetInputAODBranch();
2170 
2171  Int_t nPhot = GetInputAODBranch()->GetEntriesFast() ;
2172  Int_t nPhot2 = nPhot; // second loop
2173  Int_t minEntries = 2;
2174  Int_t last = 1; // last entry in first loop to be removed
2175 
2176  // Combine 2 detectors:
2178  {
2179  // Input from CaloTrackCorr tasks
2180  secondLoopInputData = (TClonesArray *) GetReader()->GetAODBranchList()->FindObject(fOtherDetectorInputName);
2181 
2182  // In case input from PCM or other external source
2183  if(!secondLoopInputData) secondLoopInputData = (TClonesArray *) GetReader()->GetOutputEvent()->FindObject(fOtherDetectorInputName);
2184  if(!secondLoopInputData) secondLoopInputData = (TClonesArray *) GetReader()->GetInputEvent ()->FindObject(fOtherDetectorInputName);
2185 
2186  if(!secondLoopInputData)
2187  {
2188  AliFatal(Form("Input for other detector not found in branch %s!",fOtherDetectorInputName.Data()));
2189  return ; // coverity
2190  }
2191 
2192  nPhot2 = secondLoopInputData->GetEntriesFast() ; // add one since we remove one for the normal case
2193  minEntries = 1;
2194  last = 0;
2195  }
2196 
2197  AliDebug(1,Form("Photon entries %d", nPhot));
2198 
2199  // If less than photon 2 (1) entries in the list, skip this event
2200  if ( nPhot < minEntries )
2201  {
2202  AliDebug(1,Form("nPhotons %d, cent bin %d continue to next event",nPhot, GetEventCentralityBin()));
2203 
2204  if ( GetNCentrBin() > 1 && IsHighMultiplicityAnalysisOn() )
2206 
2207  return ;
2208  }
2209  else if(fPairWithOtherDetector && nPhot2 < minEntries)
2210  {
2211  AliDebug(1,Form("nPhotons %d, cent bin %d continue to next event",nPhot, GetEventCentralityBin()));
2212 
2213  if ( GetNCentrBin() > 1 && IsHighMultiplicityAnalysisOn() )
2215 
2216  return ;
2217  }
2218 
2219  Int_t ncentr = GetNCentrBin();
2220  if(GetNCentrBin()==0) ncentr = 1;
2221 
2222  //Init variables
2223  Int_t module1 = -1;
2224  Int_t module2 = -1;
2225  Double_t vert[] = {0.0, 0.0, 0.0} ; //vertex
2226  Int_t evtIndex1 = 0 ;
2227  Int_t currentEvtIndex = -1;
2228  Int_t curCentrBin = GetEventCentralityBin();
2229  //Int_t curVzBin = GetEventVzBin();
2230  //Int_t curRPBin = GetEventRPBin();
2231  Int_t eventbin = GetEventMixBin();
2232 
2233  if(eventbin > GetNCentrBin()*GetNZvertBin()*GetNRPBin())
2234  {
2235  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",
2237  return;
2238  }
2239 
2240  // Get shower shape information of clusters
2241 // TObjArray *clusters = 0;
2242 // if (GetCalorimeter()==kEMCAL) clusters = GetEMCALClusters();
2243 // else if(GetCalorimeter()==kPHOS ) clusters = GetPHOSClusters() ;
2244 
2245  //---------------------------------
2246  // First loop on photons/clusters
2247  //---------------------------------
2248  for(Int_t i1 = 0; i1 < nPhot-last; i1++)
2249  {
2250  AliAODPWG4Particle * p1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i1)) ;
2251 
2252  // Select photons within a pT range
2253  if ( p1->Pt() < GetMinPt() || p1->Pt() > GetMaxPt() ) continue ;
2254 
2255  //printf("AliAnaPi0::MakeAnalysisFillHistograms() : cluster1 id %d/%d\n",i1,nPhot-1);
2256 
2257  // get the event index in the mixed buffer where the photon comes from
2258  // in case of mixing with analysis frame, not own mixing
2259  evtIndex1 = GetEventIndex(p1, vert) ;
2260  if ( evtIndex1 == -1 )
2261  return ;
2262  if ( evtIndex1 == -2 )
2263  continue ;
2264 
2265  // Only effective in case of mixed event frame
2266  if(TMath::Abs(vert[2]) > GetZvertexCut()) continue ; //vertex cut
2267 
2268  if (evtIndex1 != currentEvtIndex)
2269  {
2270  // Fill event bin info
2271  if( DoOwnMix() ) fhEventBin->Fill(eventbin, GetEventWeight()) ;
2272 
2274  {
2275  fhCentrality->Fill(curCentrBin, GetEventWeight());
2276 
2277  if( GetEventPlane() )
2278  fhEventPlaneResolution->Fill(curCentrBin, TMath::Cos(2.*GetEventPlane()->GetQsubRes()), GetEventWeight());
2279  }
2280 
2281  currentEvtIndex = evtIndex1 ;
2282  }
2283 
2284  //printf("AliAnaPi0::MakeAnalysisFillHistograms(): Photon 1 Evt %d Vertex : %f,%f,%f\n",evtIndex1, GetVertex(evtIndex1)[0] ,GetVertex(evtIndex1)[1],GetVertex(evtIndex1)[2]);
2285 
2286  //Get the momentum of this cluster
2287  fPhotonMom1.SetPxPyPzE(p1->Px(),p1->Py(),p1->Pz(),p1->E());
2288 
2289  //Get (Super)Module number of this cluster
2290  module1 = p1->GetSModNumber();// GetModuleNumber(p1);
2291 
2292  //------------------------------------------
2293  // Recover original cluster
2294  // Int_t iclus1 = -1 ;
2295  // AliVCluster * cluster1 = FindCluster(clusters,p1->GetCaloLabel(0),iclus1);
2296  // if(!cluster1) AliWarning("Cluster1 not found!");
2297 
2298  //---------------------------------
2299  // Second loop on photons/clusters
2300  //---------------------------------
2301  Int_t first = i1+1;
2302  if(fPairWithOtherDetector) first = 0;
2303 
2304  for(Int_t i2 = first; i2 < nPhot2; i2++)
2305  {
2306  //AliAODPWG4Particle * p2 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i2)) ;
2307  AliAODPWG4Particle * p2 = (AliAODPWG4Particle*) (secondLoopInputData->At(i2)) ;
2308 
2309  // Select photons within a pT range
2310  if ( p2->Pt() < GetMinPt() || p2->Pt() > GetMaxPt() ) continue ;
2311 
2312  //printf("AliAnaPi0::MakeAnalysisFillHistograms() : cluster2 i %d/%d\n",i2,nPhot);
2313 
2314  //In case of mixing frame, check we are not in the same event as the first cluster
2315  Int_t evtIndex2 = GetEventIndex(p2, vert) ;
2316  if ( evtIndex2 == -1 )
2317  return ;
2318  if ( evtIndex2 == -2 )
2319  continue ;
2320  if (GetMixedEvent() && (evtIndex1 == evtIndex2))
2321  continue ;
2322 
2323 // //------------------------------------------
2324 // // Recover original cluster
2325 // Int_t iclus2 = -1;
2326 // AliVCluster * cluster2 = FindCluster(clusters,p2->GetCaloLabel(0),iclus2,iclus1+1);
2327 // // start new loop from iclus1+1 to gain some time
2328 // if(!cluster2) AliWarning("Cluster2 not found!");
2329 //
2330 // // Get the TOF,l0 and ncells from the clusters
2331 // Float_t tof1 = -1;
2332 // Float_t l01 = -1;
2333 // Int_t ncell1 = 0;
2334 // if(cluster1)
2335 // {
2336 // tof1 = cluster1->GetTOF()*1e9;
2337 // l01 = cluster1->GetM02();
2338 // ncell1 = cluster1->GetNCells();
2339 // //printf("cluster1: E %2.2f (%2.2f), l0 %2.2f, tof %2.2f\n",cluster1->E(),p1->E(),l01,tof1);
2340 // }
2341 // //else printf("cluster1 not available: calo label %d / %d, cluster ID %d\n",
2342 // // p1->GetCaloLabel(0),(GetReader()->GetInputEvent())->GetNumberOfCaloClusters()-1,cluster1->GetID());
2343 //
2344 // Float_t tof2 = -1;
2345 // Float_t l02 = -1;
2346 // Int_t ncell2 = 0;
2347 // if(cluster2)
2348 // {
2349 // tof2 = cluster2->GetTOF()*1e9;
2350 // l02 = cluster2->GetM02();
2351 // ncell2 = cluster2->GetNCells();
2352 // //printf("cluster2: E %2.2f (%2.2f), l0 %2.2f, tof %2.2f\n",cluster2->E(),p2->E(),l02,tof2);
2353 // }
2354 // //else printf("cluster2 not available: calo label %d / %d, cluster ID %d\n",
2355 // // p2->GetCaloLabel(0),(GetReader()->GetInputEvent())->GetNumberOfCaloClusters()-1,cluster2->GetID());
2356 //
2357 // if(cluster1 && cluster2)
2358 // {
2359 // Double_t t12diff = tof1-tof2;
2360 // if(TMath::Abs(t12diff) > GetPairTimeCut()) continue;
2361 // }
2362 
2363  Float_t tof1 = p1->GetTime();
2364  Float_t l01 = p1->GetM02();
2365  Int_t ncell1 = p1->GetNCells();
2366  //printf("cluster1: E %2.2f, l0 %2.2f, tof %2.2f\n",p1->E(),l01,tof1);
2367 
2368  Float_t tof2 = p2->GetTime();
2369  Float_t l02 = p2->GetM02();
2370  Int_t ncell2 = p2->GetNCells();
2371  //printf("cluster2: E %2.2f, l0 %2.2f, tof %2.2f\n",p2->E(),l02,tof2);
2372 
2373  Double_t t12diff = tof1-tof2;
2374  fhEPairDiffTime->Fill((fPhotonMom1 + fPhotonMom2).Pt(), t12diff, GetEventWeight());
2375  if(TMath::Abs(t12diff) > GetPairTimeCut()) continue;
2376 
2377  //------------------------------------------
2378 
2379  //printf("AliAnaPi0::MakeAnalysisFillHistograms(): Photon 2 Evt %d Vertex : %f,%f,%f\n",evtIndex2, GetVertex(evtIndex2)[0] ,GetVertex(evtIndex2)[1],GetVertex(evtIndex2)[2]);
2380 
2381  // Get the momentum of this cluster
2382  fPhotonMom2.SetPxPyPzE(p2->Px(),p2->Py(),p2->Pz(),p2->E());
2383 
2384  // Get module number
2385  module2 = p2->GetSModNumber(); //GetModuleNumber(p2);
2386 
2387  //---------------------------------
2388  // Get pair kinematics
2389  //---------------------------------
2390  Double_t m = (fPhotonMom1 + fPhotonMom2).M() ;
2391  Double_t pt = (fPhotonMom1 + fPhotonMom2).Pt();
2392  Double_t deta = fPhotonMom1.Eta() - fPhotonMom2.Eta();
2393  Double_t dphi = fPhotonMom1.Phi() - fPhotonMom2.Phi();
2394  Double_t a = TMath::Abs(p1->E()-p2->E())/(p1->E()+p2->E()) ;
2395 
2396  AliDebug(2,Form("E: fPhotonMom1 %f, fPhotonMom2 %f; Pair: pT %f, mass %f, a %f", p1->E(), p2->E(), (fPhotonMom1 + fPhotonMom2).E(),m,a));
2397 
2398  //--------------------------------
2399  // Opening angle selection
2400  //--------------------------------
2401  // Check if opening angle is too large or too small compared to what is expected
2402  Double_t angle = fPhotonMom1.Angle(fPhotonMom2.Vect());
2404  {
2405  AliDebug(2,Form("Real pair angle %f (deg) not in E %f window",RadToDeg(angle), (fPhotonMom1+fPhotonMom2).E()));
2406  continue;
2407  }
2408 
2409  if(fUseAngleCut && (angle < fAngleCut || angle > fAngleMaxCut))
2410  {
2411  AliDebug(2,Form("Real pair cut %f < angle %f < cut %f (deg)",RadToDeg(fAngleCut), RadToDeg(angle), RadToDeg(fAngleMaxCut)));
2412  continue;
2413  }
2414 
2415  //-------------------------------------------------------------------------------------------------
2416  // Fill module dependent histograms, put a cut on assymmetry on the first available cut in the array
2417  //-------------------------------------------------------------------------------------------------
2418  if(a < fAsymCuts[0] && fFillSMCombinations)
2419  {
2421  {
2422  if(module1==module2 && module1 >=0 && module1<fNModules)
2423  fhReMod[module1]->Fill(pt, m, GetEventWeight()) ;
2424 
2425  if (GetCalorimeter() == kEMCAL )
2426  {
2427  // Same sector
2428  Int_t j=0;
2429  for(Int_t i = 0; i < fNModules/2; i++)
2430  {
2431  j=2*i;
2432  if((module1==j && module2==j+1) || (module1==j+1 && module2==j)) fhReSameSectorEMCALMod[i]->Fill(pt, m, GetEventWeight()) ;
2433  }
2434 
2435  // Same side
2436  for(Int_t i = 0; i < fNModules-2; i++){
2437  if((module1==i && module2==i+2) || (module1==i+2 && module2==i)) fhReSameSideEMCALMod[i]->Fill(pt, m, GetEventWeight());
2438  }
2439  } // EMCAL
2440  else
2441  { // PHOS
2442  if((module1==0 && module2==1) || (module1==1 && module2==0)) fhReDiffPHOSMod[0]->Fill(pt, m, GetEventWeight()) ;
2443  if((module1==0 && module2==2) || (module1==2 && module2==0)) fhReDiffPHOSMod[1]->Fill(pt, m, GetEventWeight()) ;
2444  if((module1==1 && module2==2) || (module1==2 && module2==1)) fhReDiffPHOSMod[2]->Fill(pt, m, GetEventWeight()) ;
2445  } // PHOS
2446  }
2447  else
2448  {
2449  Float_t phi1 = GetPhi(fPhotonMom1.Phi());
2450  Float_t phi2 = GetPhi(fPhotonMom2.Phi());
2451  Bool_t etaside = 0;
2452  if( (p1->GetDetectorTag()==kEMCAL && fPhotonMom1.Eta() < 0)
2453  || (p2->GetDetectorTag()==kEMCAL && fPhotonMom2.Eta() < 0)) etaside = 1;
2454 
2455  if ( phi1 > DegToRad(260) && phi2 > DegToRad(260) && phi1 < DegToRad(280) && phi2 < DegToRad(280)) fhReSameSectorDCALPHOSMod[0+etaside]->Fill(pt, m, GetEventWeight());
2456  else if ( phi1 > DegToRad(280) && phi2 > DegToRad(280) && phi1 < DegToRad(300) && phi2 < DegToRad(300)) fhReSameSectorDCALPHOSMod[2+etaside]->Fill(pt, m, GetEventWeight());
2457  else if ( phi1 > DegToRad(300) && phi2 > DegToRad(300) && phi1 < DegToRad(320) && phi2 < DegToRad(320)) fhReSameSectorDCALPHOSMod[4+etaside]->Fill(pt, m, GetEventWeight());
2458  else if ( (phi1 > DegToRad(260) && phi2 > DegToRad(280) && phi1 < DegToRad(280) && phi2 < DegToRad(300))
2459  || (phi1 > DegToRad(280) && phi2 > DegToRad(260) && phi1 < DegToRad(300) && phi2 < DegToRad(280))) fhReDiffSectorDCALPHOSMod[0+etaside]->Fill(pt, m, GetEventWeight());
2460  else if ( (phi1 > DegToRad(280) && phi2 > DegToRad(300) && phi1 < DegToRad(300) && phi2 < DegToRad(320))
2461  || (phi1 > DegToRad(300) && phi2 > DegToRad(280) && phi1 < DegToRad(320) && phi2 < DegToRad(300))) fhReDiffSectorDCALPHOSMod[2+etaside]->Fill(pt, m, GetEventWeight());
2462  else if ( (phi1 > DegToRad(260) && phi2 > DegToRad(300) && phi1 < DegToRad(280) && phi2 < DegToRad(320))
2463  || (phi1 > DegToRad(300) && phi2 > DegToRad(260) && phi1 < DegToRad(320) && phi2 < DegToRad(280))) fhReDiffSectorDCALPHOSMod[4+etaside]->Fill(pt, m, GetEventWeight());
2464  else fhReDiffSectorDCALPHOSMod[6+etaside]->Fill(pt, m, GetEventWeight());
2465  }
2466  } // Fill SM combinations
2467 
2468  // In case we want only pairs in same (super) module, check their origin.
2469  Bool_t ok = kTRUE;
2470 
2471  if(fSameSM)
2472  {
2474  {
2475  if(module1!=module2) ok=kFALSE;
2476  }
2477  else // PHOS and DCal in same sector
2478  {
2479  Float_t phi1 = GetPhi(fPhotonMom1.Phi());
2480  Float_t phi2 = GetPhi(fPhotonMom2.Phi());
2481  ok=kFALSE;
2482  if ( phi1 > DegToRad(260) && phi2 > DegToRad(260) && phi1 < DegToRad(280) && phi2 < DegToRad(280)) ok = kTRUE;
2483  else if ( phi1 > DegToRad(280) && phi2 > DegToRad(280) && phi1 < DegToRad(300) && phi2 < DegToRad(300)) ok = kTRUE;
2484  else if ( phi1 > DegToRad(300) && phi2 > DegToRad(300) && phi1 < DegToRad(320) && phi2 < DegToRad(320)) ok = kTRUE;
2485  }
2486  } // Pair only in same SM
2487 
2488  if(!ok) continue;
2489 
2490  //
2491  // Fill histograms with selected cluster pairs
2492  //
2493 
2494  // Check if one of the clusters comes from a conversion
2495  if(fCheckConversion)
2496  {
2497  if (p1->IsTagged() && p2->IsTagged()) fhReConv2->Fill(pt, m, GetEventWeight());
2498  else if(p1->IsTagged() || p2->IsTagged()) fhReConv ->Fill(pt, m, GetEventWeight());
2499  }
2500 
2501  // Fill shower shape cut histograms
2503  {
2504  if ( l01 > 0.01 && l01 < 0.4 &&
2505  l02 > 0.01 && l02 < 0.4 ) fhReSS[0]->Fill(pt, m, GetEventWeight()); // Tight
2506  else if( l01 > 0.4 && l02 > 0.4 ) fhReSS[1]->Fill(pt, m, GetEventWeight()); // Loose
2507  else if( l01 > 0.01 && l01 < 0.4 && l02 > 0.4 ) fhReSS[2]->Fill(pt, m, GetEventWeight()); // Both
2508  else if( l02 > 0.01 && l02 < 0.4 && l01 > 0.4 ) fhReSS[2]->Fill(pt, m, GetEventWeight()); // Both
2509  }
2510 
2511  // Fill histograms for different bad channel distance, centrality, assymmetry cut and pid bit
2512  for(Int_t ipid=0; ipid<fNPIDBits; ipid++)
2513  {
2514  if((p1->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton)) && (p2->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton)))
2515  {
2516  for(Int_t iasym=0; iasym < fNAsymCuts; iasym++)
2517  {
2518  if(a < fAsymCuts[iasym])
2519  {
2520  Int_t index = ((curCentrBin*fNPIDBits)+ipid)*fNAsymCuts + iasym;
2521  //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);
2522 
2523  if(index < 0 || index >= ncentr*fNPIDBits*fNAsymCuts) continue ;
2524 
2525  fhRe1 [index]->Fill(pt, m, GetEventWeight());
2526  if(fMakeInvPtPlots)fhReInvPt1[index]->Fill(pt, m, 1./pt * GetEventWeight()) ;
2527 
2528  if(fFillBadDistHisto)
2529  {
2530  if(p1->DistToBad()>0 && p2->DistToBad()>0)
2531  {
2532  fhRe2 [index]->Fill(pt, m, GetEventWeight()) ;
2533  if(fMakeInvPtPlots)fhReInvPt2[index]->Fill(pt, m, 1./pt * GetEventWeight()) ;
2534 
2535  if(p1->DistToBad()>1 && p2->DistToBad()>1)
2536  {
2537  fhRe3 [index]->Fill(pt, m, GetEventWeight()) ;
2538  if(fMakeInvPtPlots)fhReInvPt3[index]->Fill(pt, m, 1./pt * GetEventWeight()) ;
2539  }// bad 3
2540  }// bad2
2541  }// Fill bad dist histos
2542  }//assymetry cut
2543  }// asymmetry cut loop
2544  }// bad 1
2545  }// pid bit loop
2546 
2547  // Fill histograms with opening angle
2548  if(fFillAngleHisto)
2549  {
2550  fhRealOpeningAngle ->Fill(pt, angle, GetEventWeight());
2551  fhRealCosOpeningAngle->Fill(pt, TMath::Cos(angle), GetEventWeight());
2552  }
2553 
2554  // Fill histograms with pair assymmetry
2556  {
2557  fhRePtAsym->Fill(pt, a, GetEventWeight());
2558  if(m > 0.10 && m < 0.17) fhRePtAsymPi0->Fill(pt, a, GetEventWeight());
2559  if(m > 0.45 && m < 0.65) fhRePtAsymEta->Fill(pt, a, GetEventWeight());
2560  }
2561 
2562  //---------
2563  // MC data
2564  //---------
2565  // Do some MC checks on the origin of the pair, is there any common ancestor and if there is one, who?
2566  if(IsDataMC())
2567  {
2568  if(GetMCAnalysisUtils()->CheckTagBit(p1->GetTag(),AliMCAnalysisUtils::kMCConversion) &&
2570  {
2571  fhReMCFromConversion->Fill(pt, m, GetEventWeight());
2572  }
2573  else if(!GetMCAnalysisUtils()->CheckTagBit(p1->GetTag(),AliMCAnalysisUtils::kMCConversion) &&
2575  {
2576  fhReMCFromNotConversion->Fill(pt, m, GetEventWeight());
2577  }
2578  else
2579  {
2580  fhReMCFromMixConversion->Fill(pt, m, GetEventWeight());
2581  }
2582 
2583  if(fFillOriginHisto)
2584  FillMCVersusRecDataHistograms(p1->GetLabel(), p2->GetLabel(),p1->Pt(), p2->Pt(),ncell1, ncell2, m, pt, a,deta, dphi);
2585  }
2586 
2587  //-----------------------
2588  // Multi cuts analysis
2589  //-----------------------
2590  if(fMultiCutAna)
2591  {
2592  // Histograms for different PID bits selection
2593  for(Int_t ipid=0; ipid<fNPIDBits; ipid++)
2594  {
2595  if(p1->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton) &&
2596  p2->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton)) fhRePIDBits[ipid]->Fill(pt, m, GetEventWeight()) ;
2597 
2598  //printf("ipt %d, ipid%d, name %s\n",ipt, ipid, fhRePtPIDCuts[ipt*fNPIDBitsBits+ipid]->GetName());
2599  } // pid bit cut loop
2600 
2601  // Several pt,ncell and asymmetry cuts
2602  for(Int_t ipt = 0; ipt < fNPtCuts; ipt++)
2603  {
2604  for(Int_t icell = 0; icell < fNCellNCuts; icell++)
2605  {
2606  for(Int_t iasym = 0; iasym < fNAsymCuts; iasym++)
2607  {
2608  Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
2609  if(p1->E() > fPtCuts[ipt] && p2->E() > fPtCuts[ipt] &&
2610  a < fAsymCuts[iasym] &&
2611  ncell1 >= fCellNCuts[icell] && ncell2 >= fCellNCuts[icell])
2612  {
2613  fhRePtNCellAsymCuts[index]->Fill(pt, m, GetEventWeight()) ;
2614  //printf("ipt %d, icell%d, iasym %d, name %s\n",ipt, icell, iasym, fhRePtNCellAsymCuts[((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym]->GetName());
2615  if(fFillSMCombinations && module1==module2)
2616  {
2617  fhRePtNCellAsymCutsSM[module1][index]->Fill(pt, m, GetEventWeight()) ;
2618  }
2619  }
2620  }// pid bit cut loop
2621  }// icell loop
2622  }// pt cut loop
2623 
2624  if(GetHistogramRanges()->GetHistoTrackMultiplicityBins())
2625  {
2626  for(Int_t iasym = 0; iasym < fNAsymCuts; iasym++)
2627  {
2628  if(a < fAsymCuts[iasym]) fhRePtMult[iasym]->Fill(pt, GetTrackMultiplicity(), m, GetEventWeight()) ;
2629  }
2630  }
2631  }// multiple cuts analysis
2632  }// second same event particle
2633  }// first cluster
2634 
2635  //-------------------------------------------------------------
2636  // Mixing
2637  //-------------------------------------------------------------
2638  if(DoOwnMix())
2639  {
2640  // Recover events in with same characteristics as the current event
2641 
2642  // Check that the bin exists, if not (bad determination of RP, centrality or vz bin) do nothing
2643  if(eventbin < 0) return ;
2644 
2645  TList * evMixList=fEventsList[eventbin] ;
2646 
2647  if(!evMixList)
2648  {
2649  AliWarning(Form("Mix event list not available, bin %d",eventbin));
2650  return;
2651  }
2652 
2653  Int_t nMixed = evMixList->GetSize() ;
2654  for(Int_t ii=0; ii<nMixed; ii++)
2655  {
2656  TClonesArray* ev2= (TClonesArray*) (evMixList->At(ii));
2657  Int_t nPhot2=ev2->GetEntriesFast() ;
2658  Double_t m = -999;
2659  AliDebug(1,Form("Mixed event %d photon entries %d, centrality bin %d",ii, nPhot2, GetEventCentralityBin()));
2660 
2661  fhEventMixBin->Fill(eventbin, GetEventWeight()) ;
2662 
2663  //---------------------------------
2664  // First loop on photons/clusters
2665  //---------------------------------
2666  for(Int_t i1 = 0; i1 < nPhot; i1++)
2667  {
2668  AliAODPWG4Particle * p1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i1)) ;
2669 
2670  // Select photons within a pT range
2671  if ( p1->Pt() < GetMinPt() || p1->Pt() > GetMaxPt() ) continue ;
2672 
2673  // Not sure why this line is here
2674  //if(fSameSM && GetModuleNumber(p1)!=module1) continue;
2675 
2676  //Get kinematics of cluster and (super) module of this cluster
2677  fPhotonMom1.SetPxPyPzE(p1->Px(),p1->Py(),p1->Pz(),p1->E());
2678  module1 = GetModuleNumber(p1);
2679 
2680  //---------------------------------
2681  // Second loop on other mixed event photons/clusters
2682  //---------------------------------
2683  for(Int_t i2 = 0; i2 < nPhot2; i2++)
2684  {
2685  AliAODPWG4Particle * p2 = (AliAODPWG4Particle*) (ev2->At(i2)) ;
2686 
2687  // Select photons within a pT range
2688  if ( p2->Pt() < GetMinPt() || p2->Pt() > GetMaxPt() ) continue ;
2689 
2690  // Get kinematics of second cluster and calculate those of the pair
2691  fPhotonMom2.SetPxPyPzE(p2->Px(),p2->Py(),p2->Pz(),p2->E());
2692  m = (fPhotonMom1+fPhotonMom2).M() ;
2693  Double_t pt = (fPhotonMom1 + fPhotonMom2).Pt();
2694  Double_t a = TMath::Abs(p1->E()-p2->E())/(p1->E()+p2->E()) ;
2695 
2696  // Check if opening angle is too large or too small compared to what is expected
2697  Double_t angle = fPhotonMom1.Angle(fPhotonMom2.Vect());
2699  {
2700  AliDebug(2,Form("Mix pair angle %f (deg) not in E %f window",RadToDeg(angle), (fPhotonMom1+fPhotonMom2).E()));
2701  continue;
2702  }
2703 
2704  if(fUseAngleCut && (angle < fAngleCut || angle > fAngleMaxCut))
2705  {
2706  AliDebug(2,Form("Mix pair cut %f < angle %f < cut %f (deg)",RadToDeg(fAngleCut),RadToDeg(angle),RadToDeg(fAngleMaxCut)));
2707  continue;
2708  }
2709 
2710  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));
2711 
2712  // In case we want only pairs in same (super) module, check their origin.
2713  module2 = GetModuleNumber(p2);
2714 
2715  //-------------------------------------------------------------------------------------------------
2716  // Fill module dependent histograms, put a cut on assymmetry on the first available cut in the array
2717  //-------------------------------------------------------------------------------------------------
2718  if(a < fAsymCuts[0] && fFillSMCombinations)
2719  {
2721  {
2722  if(module1==module2 && module1 >=0 && module1<fNModules)
2723  fhMiMod[module1]->Fill(pt, m, GetEventWeight()) ;
2724 
2725  if(GetCalorimeter()==kEMCAL)
2726  {
2727  // Same sector
2728  Int_t j=0;
2729  for(Int_t i = 0; i < fNModules/2; i++)
2730  {
2731  j=2*i;
2732  if((module1==j && module2==j+1) || (module1==j+1 && module2==j)) fhMiSameSectorEMCALMod[i]->Fill(pt, m, GetEventWeight()) ;
2733  }
2734 
2735  // Same side
2736  for(Int_t i = 0; i < fNModules-2; i++)
2737  {
2738  if((module1==i && module2==i+2) || (module1==i+2 && module2==i)) fhMiSameSideEMCALMod[i]->Fill(pt, m, GetEventWeight());
2739  }
2740  } // EMCAL
2741  else
2742  { // PHOS
2743  if((module1==0 && module2==1) || (module1==1 && module2==0)) fhMiDiffPHOSMod[0]->Fill(pt, m, GetEventWeight()) ;
2744  if((module1==0 && module2==2) || (module1==2 && module2==0)) fhMiDiffPHOSMod[1]->Fill(pt, m, GetEventWeight()) ;
2745  if((module1==1 && module2==2) || (module1==2 && module2==1)) fhMiDiffPHOSMod[2]->Fill(pt, m, GetEventWeight()) ;
2746  } // PHOS
2747  }
2748  else
2749  {
2750  Float_t phi1 = GetPhi(fPhotonMom1.Phi());
2751  Float_t phi2 = GetPhi(fPhotonMom2.Phi());
2752  Bool_t etaside = 0;
2753  if( (p1->GetDetectorTag()==kEMCAL && fPhotonMom1.Eta() < 0)
2754  || (p2->GetDetectorTag()==kEMCAL && fPhotonMom2.Eta() < 0)) etaside = 1;
2755 
2756  if ( phi1 > DegToRad(260) && phi2 > DegToRad(260) && phi1 < DegToRad(280) && phi2 < DegToRad(280)) fhMiSameSectorDCALPHOSMod[0+etaside]->Fill(pt, m, GetEventWeight());
2757  else if ( phi1 > DegToRad(280) && phi2 > DegToRad(280) && phi1 < DegToRad(300) && phi2 < DegToRad(300)) fhMiSameSectorDCALPHOSMod[2+etaside]->Fill(pt, m, GetEventWeight());
2758  else if ( phi1 > DegToRad(300) && phi2 > DegToRad(300) && phi1 < DegToRad(320) && phi2 < DegToRad(320)) fhMiSameSectorDCALPHOSMod[4+etaside]->Fill(pt, m, GetEventWeight());
2759  else if ( (phi1 > DegToRad(260) && phi2 > DegToRad(280) && phi1 < DegToRad(280) && phi2 < DegToRad(300))
2760  || (phi1 > DegToRad(280) && phi2 > DegToRad(260) && phi1 < DegToRad(300) && phi2 < DegToRad(280))) fhMiDiffSectorDCALPHOSMod[0+etaside]->Fill(pt, m, GetEventWeight());
2761  else if ( (phi1 > DegToRad(280) && phi2 > DegToRad(300) && phi1 < DegToRad(300) && phi2 < DegToRad(320))
2762  || (phi1 > DegToRad(300) && phi2 > DegToRad(280) && phi1 < DegToRad(320) && phi2 < DegToRad(300))) fhMiDiffSectorDCALPHOSMod[2+etaside]->Fill(pt, m, GetEventWeight());
2763  else if ( (phi1 > DegToRad(260) && phi2 > DegToRad(300) && phi1 < DegToRad(280) && phi2 < DegToRad(320))
2764  || (phi1 > DegToRad(300) && phi2 > DegToRad(260) && phi1 < DegToRad(320) && phi2 < DegToRad(280))) fhMiDiffSectorDCALPHOSMod[4+etaside]->Fill(pt, m, GetEventWeight());
2765  else fhMiDiffSectorDCALPHOSMod[6+etaside]->Fill(pt, m, GetEventWeight());
2766  }
2767  }
2768 
2769  Bool_t ok = kTRUE;
2770  if(fSameSM)
2771  {
2773  {
2774  if(module1!=module2) ok=kFALSE;
2775  }
2776  else // PHOS and DCal in same sector
2777  {
2778  Float_t phi1 = GetPhi(fPhotonMom1.Phi());
2779  Float_t phi2 = GetPhi(fPhotonMom2.Phi());
2780  ok=kFALSE;
2781  if ( phi1 > DegToRad(260) && phi2 > DegToRad(260) && phi1 < DegToRad(280) && phi2 < DegToRad(280)) ok = kTRUE;
2782  else if ( phi1 > DegToRad(280) && phi2 > DegToRad(280) && phi1 < DegToRad(300) && phi2 < DegToRad(300)) ok = kTRUE;
2783  else if ( phi1 > DegToRad(300) && phi2 > DegToRad(300) && phi1 < DegToRad(320) && phi2 < DegToRad(320)) ok = kTRUE;
2784  }
2785  } // Pair only in same SM
2786 
2787  if(!ok) continue ;
2788 
2789  //
2790  // Do the final histograms with the selected clusters
2791  //
2792 
2793  // Check if one of the clusters comes from a conversion
2794  if(fCheckConversion)
2795  {
2796  if (p1->IsTagged() && p2->IsTagged()) fhMiConv2->Fill(pt, m, GetEventWeight());
2797  else if(p1->IsTagged() || p2->IsTagged()) fhMiConv ->Fill(pt, m, GetEventWeight());
2798  }
2799 
2800  // Fill histograms for different bad channel distance, centrality, assymmetry cut and pid bit
2801  for(Int_t ipid=0; ipid<fNPIDBits; ipid++)
2802  {
2803  if((p1->IsPIDOK(ipid,AliCaloPID::kPhoton)) && (p2->IsPIDOK(ipid,AliCaloPID::kPhoton)))
2804  {
2805  for(Int_t iasym=0; iasym < fNAsymCuts; iasym++)
2806  {
2807  if(a < fAsymCuts[iasym])
2808  {
2809  Int_t index = ((curCentrBin*fNPIDBits)+ipid)*fNAsymCuts + iasym;
2810 
2811  if(index < 0 || index >= ncentr*fNPIDBits*fNAsymCuts) continue ;
2812 
2813  fhMi1[index]->Fill(pt, m, GetEventWeight()) ;
2814 
2815  if(fMakeInvPtPlots)fhMiInvPt1[index]->Fill(pt, m, 1./pt * GetEventWeight()) ;
2816 
2817  if(fFillBadDistHisto)
2818  {
2819  if(p1->DistToBad()>0 && p2->DistToBad()>0)
2820  {
2821  fhMi2[index]->Fill(pt, m, GetEventWeight()) ;
2822  if(fMakeInvPtPlots)fhMiInvPt2[index]->Fill(pt, m, 1./pt * GetEventWeight()) ;
2823 
2824  if(p1->DistToBad()>1 && p2->DistToBad()>1)
2825  {
2826  fhMi3[index]->Fill(pt, m, GetEventWeight()) ;
2827  if(fMakeInvPtPlots)fhMiInvPt3[index]->Fill(pt, m, 1./pt * GetEventWeight()) ;
2828  }
2829  }
2830  }// Fill bad dist histo
2831  }//Asymmetry cut
2832  }// Asymmetry loop
2833  }//PID cut
2834  }//loop for histograms
2835 
2836  //-----------------------
2837  // Multi cuts analysis
2838  //-----------------------
2839  if(fMultiCutAna)
2840  {
2841  // Several pt,ncell and asymmetry cuts
2842 
2843  for(Int_t ipt=0; ipt<fNPtCuts; ipt++)
2844  {
2845  for(Int_t icell=0; icell<fNCellNCuts; icell++)
2846  {
2847  for(Int_t iasym=0; iasym<fNAsymCuts; iasym++)
2848  {
2849  Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
2850  if(p1->Pt() > fPtCuts[ipt] && p2->Pt() > fPtCuts[ipt] &&
2851  a < fAsymCuts[iasym] //&&
2852  //p1->GetBtag() >= fCellNCuts[icell] && p2->GetBtag() >= fCellNCuts[icell] // trick, correct it.
2853  )
2854  {
2855  fhMiPtNCellAsymCuts[index]->Fill(pt, m, GetEventWeight()) ;
2856  //printf("ipt %d, icell%d, iasym %d, name %s\n",ipt, icell, iasym, fhRePtNCellAsymCuts[((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym]->GetName());
2857  }
2858  }// pid bit cut loop
2859  }// icell loop
2860  }// pt cut loop
2861  } // Multi cut ana
2862 
2863  // Fill histograms with opening angle
2864  if(fFillAngleHisto)
2865  {
2866  fhMixedOpeningAngle ->Fill(pt, angle, GetEventWeight());
2867  fhMixedCosOpeningAngle->Fill(pt, TMath::Cos(angle), GetEventWeight());
2868  }
2869  }// second cluster loop
2870  }//first cluster loop
2871  }//loop on mixed events
2872 
2873  //--------------------------------------------------------
2874  // Add the current event to the list of events for mixing
2875  //--------------------------------------------------------
2876 
2877  //TClonesArray *currentEvent = new TClonesArray(*GetInputAODBranch());
2878  TClonesArray *currentEvent = new TClonesArray(*secondLoopInputData);
2879 
2880  // Add current event to buffer and Remove redundant events
2881  if( currentEvent->GetEntriesFast() > 0 )
2882  {
2883  evMixList->AddFirst(currentEvent) ;
2884  currentEvent=0 ; //Now list of particles belongs to buffer and it will be deleted with buffer
2885  if( evMixList->GetSize() >= GetNMaxEvMix() )
2886  {
2887  TClonesArray * tmp = (TClonesArray*) (evMixList->Last()) ;
2888  evMixList->RemoveLast() ;
2889  delete tmp ;
2890  }
2891  }
2892  else
2893  { // empty event
2894  delete currentEvent ;
2895  currentEvent=0 ;
2896  }
2897  }// DoOwnMix
2898 
2899  AliDebug(1,"End fill histograms");
2900 }
2901 
2902 //________________________________________________________________________
2906 //________________________________________________________________________
2907 Int_t AliAnaPi0::GetEventIndex(AliAODPWG4Particle * part, Double_t * vert)
2908 {
2909  Int_t evtIndex = -1 ;
2910 
2911  if(GetReader()->GetDataType()!=AliCaloTrackReader::kMC)
2912  {
2913  if (GetMixedEvent())
2914  {
2915  evtIndex = GetMixedEvent()->EventIndexForCaloCluster(part->GetCaloLabel(0)) ;
2916  GetVertex(vert,evtIndex);
2917 
2918  if(TMath::Abs(vert[2])> GetZvertexCut())
2919  evtIndex = -2 ; //Event can not be used (vertex, centrality,... cuts not fulfilled)
2920  }
2921  else
2922  {
2923  // Single event
2924  GetVertex(vert);
2925 
2926  if(TMath::Abs(vert[2])> GetZvertexCut())
2927  evtIndex = -1 ; //Event can not be used (vertex, centrality,... cuts not fulfilled)
2928  else
2929  evtIndex = 0 ;
2930  }
2931  } // No MC reader
2932  else
2933  {
2934  evtIndex = 0;
2935  vert[0] = 0. ;
2936  vert[1] = 0. ;
2937  vert[2] = 0. ;
2938  }
2939 
2940  return evtIndex ;
2941 }
2942 
TH2F * fhPrimEtaY
! Rapidity distribution of primary particles vs pT
Definition: AliAnaPi0.h:387
Float_t GetHistoPtMax() const
TH2F ** fhRe3
[GetNCentrBin()*fNPIDBits*fNAsymCuts]
Definition: AliAnaPi0.h:285
TH2F * fhPrimEtaAccPtCentrality
! primary eta with accepted daughters reconstructed centrality vs pT
Definition: AliAnaPi0.h:400
Int_t pdg
ClassImp(AliAnalysisTaskTriggerRates) AliAnalysisTaskTriggerRates
TH2F * fhPrimPi0Y
! Rapidity distribution of primary particles vs pT
Definition: AliAnaPi0.h:364
TLorentzVector fPhotonMom1Boost
! Photon cluster momentum, temporary array
Definition: AliAnaPi0.h:208
Int_t fNModules
[GetNCentrBin()*GetNZvertBin()*GetNRPBin()]
Definition: AliAnaPi0.h:170
TH2F * fhPrimEtaYetaYcut
! PseudoRapidity distribution of primary particles vs pT, Y<1
Definition: AliAnaPi0.h:390
TH1F * fhPrimPi0AccPt
! Spectrum of primary with accepted daughters
Definition: AliAnaPi0.h:362
Float_t GetHistoPtMin() const
TH2F * fhPrimEtaAccPhi
! Azimutal distribution of primary with accepted daughters vs pT
Definition: AliAnaPi0.h:393
TLorentzVector fPi0Mom
! Pi0 cluster momentum, temporary array
Definition: AliAnaPi0.h:210
TH2F * fhReMCFromMixConversion
! Invariant mass of 2 clusters one from conversion and the other not
Definition: AliAnaPi0.h:451
Int_t fPIDBits[10]
Array with different PID bits.
Definition: AliAnaPi0.h:187
virtual void AddToHistogramsName(TString add)
TH2F * fhReConv
[8]
Definition: AliAnaPi0.h:264
Bool_t fFillArmenterosThetaStar
Fill armenteros histograms.
Definition: AliAnaPi0.h:199
Int_t fNPIDBits
Number of possible PID bit combinations.
Definition: AliAnaPi0.h:186
TLorentzVector fPhotonMom2
! Photon cluster momentum, temporary array
Definition: AliAnaPi0.h:209
const char * title
Definition: MakeQAPdf.C:26
TH2F ** fhMiDiffSectorDCALPHOSMod
[6]
Definition: AliAnaPi0.h:260
Int_t fCellNCuts[10]
Array with different cell number cluster cuts.
Definition: AliAnaPi0.h:185
TH2F * fhMCOrgDeltaEta[13]
! Delta Eta vs pt of real pairs, check common origin of pair
Definition: AliAnaPi0.h:416
TH2F * fhReSS[3]
[fNAsymCuts]
Definition: AliAnaPi0.h:332
TH2F * fhPrimEtaCosOpeningAngle
! Cosinus of opening angle of pair version pair energy, eta primaries
Definition: AliAnaPi0.h:397
TH2F * fhPrimPi0PtStatus
! Spectrum of generated pi0 vs pi0 status
Definition: AliAnaPi0.h:408
Bool_t fCheckConversion
Fill histograms with tagged photons as conversion.
Definition: AliAnaPi0.h:193
virtual Float_t GetPairTimeCut() const
Time cut in ns.
virtual void GetVertex(Double_t vertex[3]) const
Int_t GetHistoTrackMultiplicityMax() const
TH1F * fhPrimPi0E
! Spectrum of Primary
Definition: AliAnaPi0.h:359
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:173
TH2F ** fhMiSameSectorEMCALMod
[fNModules-2]
Definition: AliAnaPi0.h:251
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:344
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:204
Float_t DegToRad(Float_t deg) const
TH2F * fhPrimPi0PtCentrality
! primary pi0 reconstructed centrality vs pT
Definition: AliAnaPi0.h:375
TH2F * fhPrimEtaPtCentrality
! primary eta reconstructed centrality vs pT
Definition: AliAnaPi0.h:398
virtual Float_t GetZvertexCut() const
Maximal number of events for mixin.
virtual Double_t GetEventPlaneAngle() const
TH2F * fhReMCFromConversion
! Invariant mass of 2 clusters originated in conversions
Definition: AliAnaPi0.h:449
TH2F ** fhReSameSideEMCALMod
[fNModules]
Definition: AliAnaPi0.h:230
TH1F * fhPrimPi0AccPtPhotonCuts
! Spectrum of primary with accepted daughters, photon pt or angle cuts
Definition: AliAnaPi0.h:363
TH2F * fhPrimPi0OpeningAngle
! Opening angle of pair versus pair energy, primaries
Definition: AliAnaPi0.h:371
virtual AliNeutralMesonSelection * GetNeutralMesonSelection()
TH2F * fhRealOpeningAngle
! Opening angle of pair versus pair energy
Definition: AliAnaPi0.h:351
TH2F * fhMCOrgMass[13]
! Mass vs pt of real pairs, check common origin of pair
Definition: AliAnaPi0.h:414
Int_t GetHistoMassBins() const
Int_t GetHistoPhiBins() const
TH2F * fhRePtAsymPi0
! REAL two-photon pt vs asymmetry, close to pi0 mass
Definition: AliAnaPi0.h:337
TLorentzVector fPhotonMom1
! Photon cluster momentum, temporary array
Definition: AliAnaPi0.h:207
TH2F * fhRealCosOpeningAngle
! Cosinus of opening angle of pair version pair energy
Definition: AliAnaPi0.h:352
TH2F ** fhRe1
REAL two-photon invariant mass distribution for different centralities and Asymmetry.
Definition: AliAnaPi0.h:270
Float_t GetHistoMassMin() const
TH2F * fhPrimPi0ProdVertex
! Spectrum of primary pi0 vs production vertex
Definition: AliAnaPi0.h:446
TH2F ** fhReDiffPHOSMod
[fNModules/2]
Definition: AliAnaPi0.h:236
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:265
Int_t GetHistoTrackMultiplicityMin() const
Float_t fAsymCuts[10]
Array with different assymetry cuts.
Definition: AliAnaPi0.h:183
Float_t fPtCuts[10]
Array with different pt cuts.
Definition: AliAnaPi0.h:181
TH2F * fhMCOrgAsym[13]
! Asymmetry vs pt of real pairs, check common origin of pair
Definition: AliAnaPi0.h:415
TH2F ** fhRePtNCellAsymCuts
[GetNCentrBin()*fNPIDBits*fNAsymCuts]
Definition: AliAnaPi0.h:318
TH2F * fhPrimPi0PtOrigin
! Spectrum of generated pi0 vs mother
Definition: AliAnaPi0.h:405
Float_t GetHistoDiffTimeMin() const
Bool_t IsMCParticleInCalorimeterAcceptance(Int_t calo, TParticle *particle)
TH2F ** fhMCPi0MassPtRec
Real pi0 pairs, reconstructed mass vs reconstructed pt of original pair.
Definition: AliAnaPi0.h:422
TH2F * fhPrimPi0AccPtCentrality
! primary pi0 with accepted daughters reconstructed centrality vs pT
Definition: AliAnaPi0.h:377
Float_t GetHistoPhiMin() const
Int_t fNAsymCuts
Number of assymmetry cuts.
Definition: AliAnaPi0.h:182
Float_t GetHistoDiffTimeMax() const
TH2F * fhMixedCosOpeningAngle
! Cosinus of opening angle of pair version pair energy
Definition: AliAnaPi0.h:354
TVector3 fProdVertex
! Production vertex, temporary array
Definition: AliAnaPi0.h:211
TH2F ** fhReInvPt2
[GetNCentrBin()*fNPIDBits*fNAsymCuts]
Definition: AliAnaPi0.h:301
TH2F ** fhMiDiffPHOSMod
[fNModules/2]
Definition: AliAnaPi0.h:254
TH2F ** fhMi2
[GetNCentrBin()*fNPIDBits*fNAsymCuts]
Definition: AliAnaPi0.h:281
TList * GetCreateOutputObjects()
Definition: AliAnaPi0.cxx:240
TH2F ** fhReInvPt3
[GetNCentrBin()*fNPIDBits*fNAsymCuts]
Definition: AliAnaPi0.h:309
Bool_t fFillOriginHisto
Fill histograms depending on their origin.
Definition: AliAnaPi0.h:198
TH2F * fhPrimPi0Phi
! Azimutal distribution of primary particles vs pT
Definition: AliAnaPi0.h:369
TH2F * fhPrimEtaOpeningAnglePhotonCuts
! Opening angle of pair versus pair energy, eta primaries, photon pT cuts
Definition: AliAnaPi0.h:395
Bool_t IsAngleInWindow(Float_t angle, Float_t e) const
TH2F ** fhReSameSectorDCALPHOSMod
[fNModules]
Definition: AliAnaPi0.h:239
const Double_t etamin
Float_t GetHistoAsymmetryMax() const
Float_t GetHistoMassMax() const
Base class for CaloTrackCorr analysis algorithms.
Bool_t fFillAngleHisto
Fill histograms with pair opening angle.
Definition: AliAnaPi0.h:196
virtual TString GetCalorimeterString() const
Bool_t fFillBadDistHisto
Do plots for different distances to bad channels.
Definition: AliAnaPi0.h:194
virtual Bool_t IsRealCaloAcceptanceOn() const
virtual AliFiducialCut * GetFiducialCut()
virtual TClonesArray * GetInputAODBranch() const
TH2F ** fhMCEtaPtTruePtRec
[fNPtCuts*fNAsymCuts*fNCellNCuts]
Definition: AliAnaPi0.h:437
Bool_t fMultiCutAna
Do analysis with several or fixed cut.
Definition: AliAnaPi0.h:178
virtual TClonesArray * GetAODMCParticles() const
virtual AliHistogramRanges * GetHistogramRanges()
TH2F ** fhReMod
REAL two-photon invariant mass distribution for different calorimeter modules.
Definition: AliAnaPi0.h:227
Bool_t ReadStack() const
TH1F * fhCentralityNoPair
! Histogram with centrality bins with no pair
Definition: AliAnaPi0.h:345
Bool_t fFillAsymmetryHisto
Fill histograms with asymmetry vs pt.
Definition: AliAnaPi0.h:197
Int_t GetHistoDiffTimeBins() const
virtual TList * GetAODBranchList() const
Int_t GetEventIndex(AliAODPWG4Particle *part, Double_t *vert)
Definition: AliAnaPi0.cxx:2907
void FillAcceptanceHistograms()
Fill acceptance histograms if MC data is available.
Definition: AliAnaPi0.cxx:1393
TH2F ** fhReSameSectorEMCALMod
[fNModules-2]
Definition: AliAnaPi0.h:233
TH2F ** fhMCEtaMassPtTrue
[fNPtCuts*fNAsymCuts*fNCellNCuts]
Definition: AliAnaPi0.h:434
TH2F ** fhReInvPt1
[GetNCentrBin()*fNPIDBits*fNAsymCuts]
Definition: AliAnaPi0.h:294
Int_t fNCellNCuts
Number of cuts with number of cells in cluster.
Definition: AliAnaPi0.h:184
TH2F ** fhMiInvPt1
[GetNCentrBin()*fNPIDBits*fNAsymCuts]
Definition: AliAnaPi0.h:297
Int_t GetHistoAsymmetryBins() const
const Double_t ptmax
TH2F * fhPrimEtaAccYeta
! PseudoRapidity distribution of primary with accepted daughters vs pT
Definition: AliAnaPi0.h:391
Float_t fAngleMaxCut
Select pairs with opening angle smaller than a threshold.
Definition: AliAnaPi0.h:175
TH1I * fhEventMixBin
! Number of mixed pairs in a particular bin (cen,vz,rp)
Definition: AliAnaPi0.h:343
void FillArmenterosThetaStar(Int_t pdg)
Fill armenteros plots.
Definition: AliAnaPi0.cxx:1809
virtual AliEMCALGeometry * GetEMCALGeometry() const
TH2F * fhPrimPi0CosOpeningAngle
! Cosinus of opening angle of pair version pair energy, pi0 primaries
Definition: AliAnaPi0.h:374
TH2F * fhMCPi0PtOrigin
[fNPtCuts*fNAsymCuts*fNCellNCuts]
Definition: AliAnaPi0.h:439
TH2F ** fhMiSameSideEMCALMod
[fNModules]
Definition: AliAnaPi0.h:248
Bool_t fFillOnlyMCAcceptanceHisto
Do analysis only of MC kinematics input.
Definition: AliAnaPi0.h:200
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:190
virtual AliAODEvent * GetOutputEvent() const
TH2F * fhPrimEtaPtEventPlane
! primary eta reconstructed event plane vs pT
Definition: AliAnaPi0.h:399
TH1F * fhPrimPi0AccE
! Spectrum of primary with accepted daughters
Definition: AliAnaPi0.h:361
Int_t fNPtCuts
Number of pt cuts.
Definition: AliAnaPi0.h:180
TH1F * fhPrimEtaE
! Spectrum of Primary
Definition: AliAnaPi0.h:382
TH2F * fhPrimEtaAccPtEventPlane
! primary eta with accepted daughters reconstructed event plane vs pT
Definition: AliAnaPi0.h:401
virtual AliCalorimeterUtils * GetCaloUtils() const
TH2F * fhPrimNotResonancePi0PtOrigin
! Spectrum of generated pi0 vs mother
Definition: AliAnaPi0.h:407
TH2F * fhPrimEtaAccY
! Rapidity distribution of primary with accepted daughters vs pT
Definition: AliAnaPi0.h:388
TH1F * fhPrimEtaPt
! Spectrum of Primary
Definition: AliAnaPi0.h:383
Bool_t fMultiCutAnaSim
Do analysis with several or fixed cut, in the simulation related part.
Definition: AliAnaPi0.h:179
TH1F * fhPrimEtaAccPt
! Spectrum of primary with accepted daughters
Definition: AliAnaPi0.h:385
TH2F * fhArmPrimPi0[4]
! Armenteros plots for primary pi0 in 6 energy bins
Definition: AliAnaPi0.h:453
void Print(const Option_t *opt) const
Print some relevant parameters set for the analysis.
Definition: AliAnaPi0.cxx:1350
TH2F * fhPrimPi0AccPtEventPlane
! primary pi0 with accepted daughters reconstructed event plane vs pT
Definition: AliAnaPi0.h:378
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:394
TH2F ** fhMiPtNCellAsymCuts
[fNPtCuts*fNAsymCuts*fNCellNCuts]
Definition: AliAnaPi0.h:321
TH2F * fhPrimPi0OpeningAngleAsym
! Opening angle of pair versus pair E asymmetry, pi0 primaries
Definition: AliAnaPi0.h:373
TH2F * fhPrimPi0AccY
! Rapidity distribution of primary with accepted daughters vs pT
Definition: AliAnaPi0.h:365
TH2F * fhPrimPi0AccPhi
! Azimutal distribution of primary with accepted daughters vs pT
Definition: AliAnaPi0.h:370
Float_t fAngleCut
Select pairs with opening angle larger than a threshold.
Definition: AliAnaPi0.h:174
virtual Double_t GetEventWeight() const
TString fOtherDetectorInputName
String with name of extra detector data.
Definition: AliAnaPi0.h:205
TH2F * fhPrimPi0AccYeta
! PseudoRapidity distribution of primary with accepted daughters vs pT
Definition: AliAnaPi0.h:368
virtual Int_t GetNCentrBin() const
Number of bins in reaction plain.
Int_t GetHistoTrackMultiplicityBins() const
TH2F ** fhMi1
[GetNCentrBin()*fNPIDBits*fNAsymCuts]
Definition: AliAnaPi0.h:273
TH1F * fhPrimEtaAccPtPhotonCuts
! Spectrum of primary with accepted daughters, photon pt or angle cuts
Definition: AliAnaPi0.h:386
TH2F * fhCosThStarPrimEta
! cos(theta*) plots vs E for primary eta, same as asymmetry ...
Definition: AliAnaPi0.h:456
TObjString * GetAnalysisCuts()
Save parameters used for analysis.
Definition: AliAnaPi0.cxx:192
void MakeAnalysisFillHistograms()
Definition: AliAnaPi0.cxx:2145
Float_t GetHistoEtaMin() const
Bool_t fFillSMCombinations
Fill histograms with different cluster pairs in SM combinations.
Definition: AliAnaPi0.h:192
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:450
TH2F * fhPrimPi0YetaYcut
! PseudoRapidity distribution of primary particles vs pT, Y<1
Definition: AliAnaPi0.h:367
virtual Int_t GetNMaxEvMix() const
Number of bins in track multiplicity.
TH2F * fhPrimPi0Yeta
! PseudoRapidity distribution of primary particles vs pT
Definition: AliAnaPi0.h:366
TH3F ** fhRePtMult
[fNPIDBits]
Definition: AliAnaPi0.h:330
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:441
TH2F * fhRePtAsym
! REAL two-photon pt vs asymmetry
Definition: AliAnaPi0.h:336
Float_t GetHistoEtaMax() const
TH2F * fhEPairDiffTime
! E pair vs Pair of clusters time difference vs E
Definition: AliAnaPi0.h:458
void InitParameters()
Definition: AliAnaPi0.cxx:158
TH2F * fhMCEtaPtOrigin
! Mass of reconstructed eta pairs in calorimeter vs mother origin ID.
Definition: AliAnaPi0.h:440
TH2F * fhMCPi0PtStatus
! Mass of reconstructed pi0 pairs in calorimeter vs mother status.
Definition: AliAnaPi0.h:442
Int_t GetHistoPtBins() const
TH2F * fhMiConv2
! MIXED two-photon invariant mass distribution both pair photons recombined from 2 clusters with smal...
Definition: AliAnaPi0.h:267
virtual ~AliAnaPi0()
Destructor.
Definition: AliAnaPi0.cxx:132
TH2F * fhPrimEtaProdVertex
! Spectrum of primary eta vs production vertex
Definition: AliAnaPi0.h:447
Bool_t fSameSM
Select only pairs in same SM;.
Definition: AliAnaPi0.h:191
Bool_t fCheckAccInSector
Check that the decay pi0 falls in the same SM or sector.
Definition: AliAnaPi0.h:202
TH1I * fhEventBin
! Number of real pairs in a particular bin (cen,vz,rp)
Definition: AliAnaPi0.h:342
TH1F * fhPrimPi0Pt
! Spectrum of Primary
Definition: AliAnaPi0.h:360
TH2F ** fhRe2
[GetNCentrBin()*fNPIDBits*fNAsymCuts]
Definition: AliAnaPi0.h:277
TH1F * fhPrimEtaAccE
! Spectrum of primary with accepted daughters
Definition: AliAnaPi0.h:384
TH2F ** fhMCEtaMassPtRec
[fNPtCuts*fNAsymCuts*fNCellNCuts]
Definition: AliAnaPi0.h:431
Bool_t fUseAngleCut
Select pairs depending on their opening angle.
Definition: AliAnaPi0.h:172
TH2F * fhRePtAsymEta
! REAL two-photon pt vs asymmetry, close to eta mass
Definition: AliAnaPi0.h:338
virtual Int_t GetNRPBin() const
Number of bins in vertex.
TH2F * fhEventPlaneResolution
! Histogram with Event plane resolution vs centrality
Definition: AliAnaPi0.h:347
TH2F * fhMixedOpeningAngle
! Opening angle of pair versus pair energy
Definition: AliAnaPi0.h:353
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:1872
const Double_t etamax
TH2F ** fhMCPi0PtTruePtRec
[fNPtCuts*fNAsymCuts*fNCellNCuts]
Definition: AliAnaPi0.h:428
TH2F * fhPrimEtaYeta
! PseudoRapidity distribution of primary particles vs pT
Definition: AliAnaPi0.h:389
virtual AliMCAnalysisUtils * GetMCAnalysisUtils()
Float_t GetHistoAsymmetryMin() const
TH2F * fhCosThStarPrimPi0
! cos(theta*) plots vs E for primary pi0, same as asymmetry ...
Definition: AliAnaPi0.h:455
virtual void Print(const Option_t *) const
Print some relevant parameters set for the analysis.
TH2F ** fhRePIDBits
[fNPtCuts*fNAsymCuts*fNCellNCutsfNModules]
Definition: AliAnaPi0.h:327
TH2F ** fhReDiffSectorDCALPHOSMod
[6]
Definition: AliAnaPi0.h:242
TH2F * fhMCPi0ProdVertex
! Spectrum of selected pi0 vs production vertex
Definition: AliAnaPi0.h:444
TH2F * fhPrimPi0PtEventPlane
! primary pi0 reconstructed event plane vs pT
Definition: AliAnaPi0.h:376
TH2F ** fhRePtNCellAsymCutsSM[12]
[fNPtCuts*fNAsymCuts*fNCellNCuts]
Definition: AliAnaPi0.h:324
TList ** fEventsList
Containers for photons in stored events.
Definition: AliAnaPi0.h:168
TH2F * fhMCEtaProdVertex
! Spectrum of selected eta vs production vertex
Definition: AliAnaPi0.h:445
TH2F * fhPrimEtaOpeningAngleAsym
! Opening angle of pair versus pair E asymmetry, eta primaries
Definition: AliAnaPi0.h:396
TH2F ** fhMiSameSectorDCALPHOSMod
[fNModules-1]
Definition: AliAnaPi0.h:257
Float_t GetHistoPhiMax() const
virtual AliCaloTrackReader * GetReader() const
TH2F ** fhMiInvPt3
[GetNCentrBin()*fNPIDBits*fNAsymCuts]
Definition: AliAnaPi0.h:313
AliAnaPi0()
Default Constructor. Initialized parameters with default values.
Definition: AliAnaPi0.cxx:53
Int_t GetHistoEtaBins() const
TH2F ** fhMi3
[GetNCentrBin()*fNPIDBits*fNAsymCuts]
Definition: AliAnaPi0.h:289
TH2F * fhArmPrimEta[4]
! Armenteros plots for primary eta in 6 energy bins
Definition: AliAnaPi0.h:454
Bool_t fFillSSCombinations
Do invariant mass for different combination of shower shape clusters.
Definition: AliAnaPi0.h:195
TH2F * fhPrimEtaPtOrigin
! Spectrum of generated eta vs mother
Definition: AliAnaPi0.h:406
TH2F ** fhMiInvPt2
[GetNCentrBin()*fNPIDBits*fNAsymCuts]
Definition: AliAnaPi0.h:305
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
TH2F * fhMCOrgDeltaPhi[13]
! Delta Phi vs pt of real pairs, check common origin of pair
Definition: AliAnaPi0.h:417
Int_t nptbins
TH2F ** fhMiMod
[8]
Definition: AliAnaPi0.h:245
virtual AliMixedEvent * GetMixedEvent() const
TH2F * fhReConv2
! REAL two-photon invariant mass distribution both pair photons recombined from 2 clusters with small...
Definition: AliAnaPi0.h:266
const Double_t phimin
TH2F * fhPrimEtaPhi
! Azimutal distribution of primary particles vs pT
Definition: AliAnaPi0.h:392
TH2F * fhPrimPi0OpeningAnglePhotonCuts
! Opening angle of pair versus pair energy, primaries, photon pt cuts
Definition: AliAnaPi0.h:372
TH2F ** fhMCPi0MassPtTrue
[fNPtCuts*fNAsymCuts*fNCellNCuts]
Definition: AliAnaPi0.h:425