1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 //_________________________________________________________________________
17 // Class to collect two-photon invariant mass distributions for
18 // extracting raw pi0 yield.
19 // Input is produced by AliAnaPhoton (or any other analysis producing output AliAODPWG4Particles),
20 // it will do nothing if executed alone
22 //-- Author: Dmitri Peressounko (RRC "KI")
23 //-- Adapted to CaloTrackCorr frame by Lamia Benhabib (SUBATECH)
24 //-- and Gustavo Conesa (INFN-Frascati)
25 //_________________________________________________________________________
28 // --- ROOT system ---
31 //#include "Riostream.h"
35 #include "TClonesArray.h"
36 #include "TObjString.h"
37 #include "TDatabasePDG.h"
39 //---- AliRoot system ----
40 #include "AliAnaPi0.h"
41 #include "AliCaloTrackReader.h"
42 #include "AliCaloPID.h"
44 #include "AliFiducialCut.h"
45 #include "TParticle.h"
46 #include "AliVEvent.h"
47 #include "AliESDCaloCluster.h"
48 #include "AliESDEvent.h"
49 #include "AliAODEvent.h"
50 #include "AliNeutralMesonSelection.h"
51 #include "AliMixedEvent.h"
52 #include "AliAODMCParticle.h"
55 #include "AliPHOSGeoUtils.h"
56 #include "AliEMCALGeometry.h"
60 //______________________________________________________
61 AliAnaPi0::AliAnaPi0() : AliAnaCaloTrackCorrBaseClass(),
63 fCalorimeter(""), fNModules(12),
64 fUseAngleCut(kFALSE), fUseAngleEDepCut(kFALSE), fAngleCut(0), fAngleMaxCut(7.),
65 fMultiCutAna(kFALSE), fMultiCutAnaSim(kFALSE),
66 fNPtCuts(0), fNAsymCuts(0), fNCellNCuts(0), fNPIDBits(0),
67 fMakeInvPtPlots(kFALSE), fSameSM(kFALSE),
68 fFillSMCombinations(kFALSE), fCheckConversion(kFALSE),
69 fFillBadDistHisto(kFALSE), fFillSSCombinations(kFALSE),
70 fFillAngleHisto(kFALSE), fFillAsymmetryHisto(kFALSE), fFillOriginHisto(0), fFillArmenterosThetaStar(0),
72 fhAverTotECluster(0), fhAverTotECell(0), fhAverTotECellvsCluster(0),
73 fhEDensityCluster(0), fhEDensityCell(0), fhEDensityCellvsCluster(0),
74 fhReMod(0x0), fhReSameSideEMCALMod(0x0), fhReSameSectorEMCALMod(0x0), fhReDiffPHOSMod(0x0),
75 fhMiMod(0x0), fhMiSameSideEMCALMod(0x0), fhMiSameSectorEMCALMod(0x0), fhMiDiffPHOSMod(0x0),
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),
88 fhPrimPi0E(0x0), fhPrimPi0Pt(0x0),
89 fhPrimPi0AccE(0x0), fhPrimPi0AccPt(0x0),
90 fhPrimPi0Y(0x0), fhPrimPi0AccY(0x0),
91 fhPrimPi0Phi(0x0), fhPrimPi0AccPhi(0x0),
92 fhPrimPi0OpeningAngle(0x0), fhPrimPi0OpeningAngleAsym(0x0),fhPrimPi0CosOpeningAngle(0x0),
93 fhPrimPi0PtCentrality(0), fhPrimPi0PtEventPlane(0),
94 fhPrimPi0AccPtCentrality(0), fhPrimPi0AccPtEventPlane(0),
95 fhPrimEtaE(0x0), fhPrimEtaPt(0x0),
96 fhPrimEtaAccE(0x0), fhPrimEtaAccPt(0x0),
97 fhPrimEtaY(0x0), fhPrimEtaAccY(0x0),
98 fhPrimEtaPhi(0x0), fhPrimEtaAccPhi(0x0),
99 fhPrimEtaOpeningAngle(0x0), fhPrimEtaOpeningAngleAsym(0x0),fhPrimEtaCosOpeningAngle(0x0),
100 fhPrimEtaPtCentrality(0), fhPrimEtaPtEventPlane(0),
101 fhPrimEtaAccPtCentrality(0), fhPrimEtaAccPtEventPlane(0),
102 fhPrimPi0PtOrigin(0x0), fhPrimEtaPtOrigin(0x0),
103 fhMCOrgMass(), fhMCOrgAsym(), fhMCOrgDeltaEta(), fhMCOrgDeltaPhi(),
104 fhMCPi0MassPtRec(), fhMCPi0MassPtTrue(), fhMCPi0PtTruePtRec(),
105 fhMCEtaMassPtRec(), fhMCEtaMassPtTrue(), fhMCEtaPtTruePtRec(),
106 fhMCPi0PtOrigin(0x0), fhMCEtaPtOrigin(0x0),
107 fhReMCFromConversion(0), fhReMCFromNotConversion(0), fhReMCFromMixConversion(0),
108 fhCosThStarPrimPi0(0), fhCosThStarPrimEta(0)//,
114 for(Int_t i = 0; i < 4; i++)
121 //_____________________
122 AliAnaPi0::~AliAnaPi0()
124 // Remove event containers
126 if(DoOwnMix() && fEventsList)
128 for(Int_t ic=0; ic<GetNCentrBin(); ic++)
130 for(Int_t iz=0; iz<GetNZvertBin(); iz++)
132 for(Int_t irp=0; irp<GetNRPBin(); irp++)
134 Int_t bin = GetEventMixBin(ic,iz,irp);
135 fEventsList[bin]->Delete() ;
136 delete fEventsList[bin] ;
140 delete[] fEventsList;
145 //______________________________
146 void AliAnaPi0::InitParameters()
148 //Init parameters when first called the analysis
149 //Set default parameters
150 SetInputAODName("PWG4Particle");
152 AddToHistogramsName("AnaPi0_");
153 fNModules = 12; // set maximum to maximum number of EMCAL modules
155 fCalorimeter = "PHOS";
156 fUseAngleCut = kFALSE;
157 fUseAngleEDepCut = kFALSE;
159 fAngleMaxCut = TMath::Pi();
161 fMultiCutAna = kFALSE;
164 fPtCuts[0] = 0.; fPtCuts[1] = 0.3; fPtCuts[2] = 0.5;
165 for(Int_t i = fNPtCuts; i < 10; i++)fPtCuts[i] = 0.;
168 fAsymCuts[0] = 1.; fAsymCuts[1] = 0.7; //fAsymCuts[2] = 0.6; // fAsymCuts[3] = 0.1;
169 for(Int_t i = fNAsymCuts; i < 10; i++)fAsymCuts[i] = 0.;
172 fCellNCuts[0] = 0; fCellNCuts[1] = 1; fCellNCuts[2] = 2;
173 for(Int_t i = fNCellNCuts; i < 10; i++)fCellNCuts[i] = 0;
176 fPIDBits[0] = 0; fPIDBits[1] = 2; // fPIDBits[2] = 4; fPIDBits[3] = 6;// check, no cut, dispersion, neutral, dispersion&&neutral
177 for(Int_t i = fNPIDBits; i < 10; i++)fPIDBits[i] = 0;
182 //_______________________________________
183 TObjString * AliAnaPi0::GetAnalysisCuts()
185 //Save parameters used for analysis
186 TString parList ; //this will be list of parameters used for this analysis.
187 const Int_t buffersize = 255;
188 char onePar[buffersize] ;
189 snprintf(onePar,buffersize,"--- AliAnaPi0 ---\n") ;
191 snprintf(onePar,buffersize,"Number of bins in Centrality: %d \n",GetNCentrBin()) ;
193 snprintf(onePar,buffersize,"Number of bins in Z vert. pos: %d \n",GetNZvertBin()) ;
195 snprintf(onePar,buffersize,"Number of bins in Reac. Plain: %d \n",GetNRPBin()) ;
197 snprintf(onePar,buffersize,"Depth of event buffer: %d \n",GetNMaxEvMix()) ;
199 snprintf(onePar,buffersize,"Select pairs with their angle: %d, edep %d, min angle %2.3f, max angle %2.3f,\n",fUseAngleCut, fUseAngleEDepCut,fAngleCut,fAngleMaxCut) ;
201 snprintf(onePar,buffersize," Asymmetry cuts: n = %d, asymmetry < ",fNAsymCuts) ;
202 for(Int_t i = 0; i < fNAsymCuts; i++) snprintf(onePar,buffersize,"%s %2.2f;",onePar,fAsymCuts[i]);
204 snprintf(onePar,buffersize," PID selection bits: n = %d, PID bit =\n",fNPIDBits) ;
205 for(Int_t i = 0; i < fNPIDBits; i++) snprintf(onePar,buffersize,"%s %d;",onePar,fPIDBits[i]);
207 snprintf(onePar,buffersize,"Cuts: \n") ;
209 snprintf(onePar,buffersize,"Z vertex position: -%f < z < %f \n",GetZvertexCut(),GetZvertexCut()) ;
211 snprintf(onePar,buffersize,"Calorimeter: %s \n",fCalorimeter.Data()) ;
213 snprintf(onePar,buffersize,"Number of modules: %d \n",fNModules) ;
216 snprintf(onePar, buffersize," pT cuts: n = %d, pt > ",fNPtCuts) ;
217 for(Int_t i = 0; i < fNPtCuts; i++) snprintf(onePar,buffersize,"%s %2.2f;",onePar,fPtCuts[i]);
219 snprintf(onePar,buffersize, " N cell in cluster cuts: n = %d, nCell > ",fNCellNCuts) ;
220 for(Int_t i = 0; i < fNCellNCuts; i++) snprintf(onePar,buffersize,"%s %d;",onePar,fCellNCuts[i]);
224 return new TObjString(parList) ;
227 //_________________________________________
228 TList * AliAnaPi0::GetCreateOutputObjects()
230 // Create histograms to be saved in output file and
231 // store them in fOutputContainer
233 //create event containers
234 fEventsList = new TList*[GetNCentrBin()*GetNZvertBin()*GetNRPBin()] ;
236 for(Int_t ic=0; ic<GetNCentrBin(); ic++)
238 for(Int_t iz=0; iz<GetNZvertBin(); iz++)
240 for(Int_t irp=0; irp<GetNRPBin(); irp++)
242 Int_t bin = GetEventMixBin(ic,iz,irp);
243 fEventsList[bin] = new TList() ;
244 fEventsList[bin]->SetOwner(kFALSE);
249 TList * outputContainer = new TList() ;
250 outputContainer->SetName(GetName());
252 fhReMod = new TH2F*[fNModules] ;
253 fhMiMod = new TH2F*[fNModules] ;
255 if(fCalorimeter == "PHOS")
257 fhReDiffPHOSMod = new TH2F*[fNModules] ;
258 fhMiDiffPHOSMod = new TH2F*[fNModules] ;
262 fhReSameSectorEMCALMod = new TH2F*[fNModules/2] ;
263 fhReSameSideEMCALMod = new TH2F*[fNModules-2] ;
264 fhMiSameSectorEMCALMod = new TH2F*[fNModules/2] ;
265 fhMiSameSideEMCALMod = new TH2F*[fNModules-2] ;
269 fhRe1 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
270 fhMi1 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
271 if(fFillBadDistHisto)
273 fhRe2 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
274 fhRe3 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
275 fhMi2 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
276 fhMi3 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
280 fhReInvPt1 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
281 fhMiInvPt1 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
282 if(fFillBadDistHisto){
283 fhReInvPt2 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
284 fhReInvPt3 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
285 fhMiInvPt2 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
286 fhMiInvPt3 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
290 const Int_t buffersize = 255;
291 char key[buffersize] ;
292 char title[buffersize] ;
294 Int_t nptbins = GetHistogramRanges()->GetHistoPtBins();
295 Int_t nphibins = GetHistogramRanges()->GetHistoPhiBins();
296 Int_t netabins = GetHistogramRanges()->GetHistoEtaBins();
297 Float_t ptmax = GetHistogramRanges()->GetHistoPtMax();
298 Float_t phimax = GetHistogramRanges()->GetHistoPhiMax();
299 Float_t etamax = GetHistogramRanges()->GetHistoEtaMax();
300 Float_t ptmin = GetHistogramRanges()->GetHistoPtMin();
301 Float_t phimin = GetHistogramRanges()->GetHistoPhiMin();
302 Float_t etamin = GetHistogramRanges()->GetHistoEtaMin();
304 Int_t nmassbins = GetHistogramRanges()->GetHistoMassBins();
305 Int_t nasymbins = GetHistogramRanges()->GetHistoAsymmetryBins();
306 Float_t massmax = GetHistogramRanges()->GetHistoMassMax();
307 Float_t asymmax = GetHistogramRanges()->GetHistoAsymmetryMax();
308 Float_t massmin = GetHistogramRanges()->GetHistoMassMin();
309 Float_t asymmin = GetHistogramRanges()->GetHistoAsymmetryMin();
310 Int_t ntrmbins = GetHistogramRanges()->GetHistoTrackMultiplicityBins();
311 Int_t ntrmmax = GetHistogramRanges()->GetHistoTrackMultiplicityMax();
312 Int_t ntrmmin = GetHistogramRanges()->GetHistoTrackMultiplicityMin();
316 fhReConv = new TH2F("hReConv","Real Pair with one recombined conversion ",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
317 fhReConv->SetXTitle("p_{T} (GeV/c)");
318 fhReConv->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
319 outputContainer->Add(fhReConv) ;
321 fhReConv2 = new TH2F("hReConv2","Real Pair with 2 recombined conversion ",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
322 fhReConv2->SetXTitle("p_{T} (GeV/c)");
323 fhReConv2->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
324 outputContainer->Add(fhReConv2) ;
328 fhMiConv = new TH2F("hMiConv","Mixed Pair with one recombined conversion ",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
329 fhMiConv->SetXTitle("p_{T} (GeV/c)");
330 fhMiConv->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
331 outputContainer->Add(fhMiConv) ;
333 fhMiConv2 = new TH2F("hMiConv2","Mixed Pair with 2 recombined conversion ",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
334 fhMiConv2->SetXTitle("p_{T} (GeV/c)");
335 fhMiConv2->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
336 outputContainer->Add(fhMiConv2) ;
340 for(Int_t ic=0; ic<GetNCentrBin(); ic++)
342 for(Int_t ipid=0; ipid<fNPIDBits; ipid++)
344 for(Int_t iasym=0; iasym<fNAsymCuts; iasym++)
346 Int_t index = ((ic*fNPIDBits)+ipid)*fNAsymCuts + iasym;
347 //printf("cen %d, pid %d, asy %d, Index %d\n",ic,ipid,iasym,index);
348 //Distance to bad module 1
349 snprintf(key, buffersize,"hRe_cen%d_pidbit%d_asy%d_dist1",ic,ipid,iasym) ;
350 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 1",
351 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
352 fhRe1[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
353 fhRe1[index]->SetXTitle("p_{T} (GeV/c)");
354 fhRe1[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
355 //printf("name: %s\n ",fhRe1[index]->GetName());
356 outputContainer->Add(fhRe1[index]) ;
358 if(fFillBadDistHisto)
360 //Distance to bad module 2
361 snprintf(key, buffersize,"hRe_cen%d_pidbit%d_asy%d_dist2",ic,ipid,iasym) ;
362 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 2",
363 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
364 fhRe2[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
365 fhRe2[index]->SetXTitle("p_{T} (GeV/c)");
366 fhRe2[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
367 outputContainer->Add(fhRe2[index]) ;
369 //Distance to bad module 3
370 snprintf(key, buffersize,"hRe_cen%d_pidbit%d_asy%d_dist3",ic,ipid,iasym) ;
371 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 3",
372 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
373 fhRe3[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
374 fhRe3[index]->SetXTitle("p_{T} (GeV/c)");
375 fhRe3[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
376 outputContainer->Add(fhRe3[index]) ;
382 //Distance to bad module 1
383 snprintf(key, buffersize,"hReInvPt_cen%d_pidbit%d_asy%d_dist1",ic,ipid,iasym) ;
384 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 1",
385 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
386 fhReInvPt1[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
387 fhReInvPt1[index]->SetXTitle("p_{T} (GeV/c)");
388 fhReInvPt1[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
389 outputContainer->Add(fhReInvPt1[index]) ;
391 if(fFillBadDistHisto){
392 //Distance to bad module 2
393 snprintf(key, buffersize,"hReInvPt_cen%d_pidbit%d_asy%d_dist2",ic,ipid,iasym) ;
394 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 2",
395 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
396 fhReInvPt2[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
397 fhReInvPt2[index]->SetXTitle("p_{T} (GeV/c)");
398 fhReInvPt2[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
399 outputContainer->Add(fhReInvPt2[index]) ;
401 //Distance to bad module 3
402 snprintf(key, buffersize,"hReInvPt_cen%d_pidbit%d_asy%d_dist3",ic,ipid,iasym) ;
403 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 3",
404 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
405 fhReInvPt3[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
406 fhReInvPt3[index]->SetXTitle("p_{T} (GeV/c)");
407 fhReInvPt3[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
408 outputContainer->Add(fhReInvPt3[index]) ;
414 //Distance to bad module 1
415 snprintf(key, buffersize,"hMi_cen%d_pidbit%d_asy%d_dist1",ic,ipid,iasym) ;
416 snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 1",
417 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
418 fhMi1[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
419 fhMi1[index]->SetXTitle("p_{T} (GeV/c)");
420 fhMi1[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
421 outputContainer->Add(fhMi1[index]) ;
422 if(fFillBadDistHisto){
423 //Distance to bad module 2
424 snprintf(key, buffersize,"hMi_cen%d_pidbit%d_asy%d_dist2",ic,ipid,iasym) ;
425 snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 2",
426 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
427 fhMi2[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
428 fhMi2[index]->SetXTitle("p_{T} (GeV/c)");
429 fhMi2[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
430 outputContainer->Add(fhMi2[index]) ;
432 //Distance to bad module 3
433 snprintf(key, buffersize,"hMi_cen%d_pidbit%d_asy%d_dist3",ic,ipid,iasym) ;
434 snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 3",
435 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
436 fhMi3[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
437 fhMi3[index]->SetXTitle("p_{T} (GeV/c)");
438 fhMi3[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
439 outputContainer->Add(fhMi3[index]) ;
445 //Distance to bad module 1
446 snprintf(key, buffersize,"hMiInvPt_cen%d_pidbit%d_asy%d_dist1",ic,ipid,iasym) ;
447 snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 1",
448 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
449 fhMiInvPt1[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
450 fhMiInvPt1[index]->SetXTitle("p_{T} (GeV/c)");
451 fhMiInvPt1[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
452 outputContainer->Add(fhMiInvPt1[index]) ;
453 if(fFillBadDistHisto){
454 //Distance to bad module 2
455 snprintf(key, buffersize,"hMiInvPt_cen%d_pidbit%d_asy%d_dist2",ic,ipid,iasym) ;
456 snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 2",
457 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
458 fhMiInvPt2[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
459 fhMiInvPt2[index]->SetXTitle("p_{T} (GeV/c)");
460 fhMiInvPt2[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
461 outputContainer->Add(fhMiInvPt2[index]) ;
463 //Distance to bad module 3
464 snprintf(key, buffersize,"hMiInvPt_cen%d_pidbit%d_asy%d_dist3",ic,ipid,iasym) ;
465 snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f,dist bad 3",
466 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
467 fhMiInvPt3[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
468 fhMiInvPt3[index]->SetXTitle("p_{T} (GeV/c)");
469 fhMiInvPt3[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
470 outputContainer->Add(fhMiInvPt3[index]) ;
478 if(fFillAsymmetryHisto)
480 fhRePtAsym = new TH2F("hRePtAsym","Asymmetry vs pt, for pairs",nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax) ;
481 fhRePtAsym->SetXTitle("p_{T} (GeV/c)");
482 fhRePtAsym->SetYTitle("Asymmetry");
483 outputContainer->Add(fhRePtAsym);
485 fhRePtAsymPi0 = new TH2F("hRePtAsymPi0","Asymmetry vs pt, for pairs close to #pi^{0} mass",nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax) ;
486 fhRePtAsymPi0->SetXTitle("p_{T} (GeV/c)");
487 fhRePtAsymPi0->SetYTitle("Asymmetry");
488 outputContainer->Add(fhRePtAsymPi0);
490 fhRePtAsymEta = new TH2F("hRePtAsymEta","Asymmetry vs pt, for pairs close to #eta mass",nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax) ;
491 fhRePtAsymEta->SetXTitle("p_{T} (GeV/c)");
492 fhRePtAsymEta->SetYTitle("Asymmetry");
493 outputContainer->Add(fhRePtAsymEta);
498 fhRePIDBits = new TH2F*[fNPIDBits];
499 for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
500 snprintf(key, buffersize,"hRe_pidbit%d",ipid) ;
501 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for PIDBit=%d",fPIDBits[ipid]) ;
502 fhRePIDBits[ipid] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
503 fhRePIDBits[ipid]->SetXTitle("p_{T} (GeV/c)");
504 fhRePIDBits[ipid]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
505 outputContainer->Add(fhRePIDBits[ipid]) ;
508 fhRePtNCellAsymCuts = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
509 fhMiPtNCellAsymCuts = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
511 if(fFillSMCombinations)
512 for(Int_t iSM = 0; iSM < fNModules; iSM++) fhRePtNCellAsymCutsSM[iSM] = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
515 for(Int_t ipt=0; ipt<fNPtCuts; ipt++)
517 for(Int_t icell=0; icell<fNCellNCuts; icell++)
519 for(Int_t iasym=0; iasym<fNAsymCuts; iasym++)
521 snprintf(key, buffersize,"hRe_pt%d_cell%d_asym%d",ipt,icell,iasym) ;
522 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for pt >%2.2f, ncell>%d and asym >%1.2f ",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]) ;
523 Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
524 //printf("ipt %d, icell %d, iassym %d, index %d\n",ipt, icell, iasym, index);
525 fhRePtNCellAsymCuts[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
526 fhRePtNCellAsymCuts[index]->SetXTitle("p_{T} (GeV/c)");
527 fhRePtNCellAsymCuts[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
528 outputContainer->Add(fhRePtNCellAsymCuts[index]) ;
530 snprintf(key, buffersize,"hMi_pt%d_cell%d_asym%d",ipt,icell,iasym) ;
531 snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for pt >%2.2f, ncell>%d and asym >%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]) ;
532 fhMiPtNCellAsymCuts[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
533 fhMiPtNCellAsymCuts[index]->SetXTitle("p_{T} (GeV/c)");
534 fhMiPtNCellAsymCuts[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
535 outputContainer->Add(fhMiPtNCellAsymCuts[index]) ;
537 if(fFillSMCombinations)
539 for(Int_t iSM = 0; iSM < fNModules; iSM++)
541 snprintf(key, buffersize,"hRe_pt%d_cell%d_asym%d_SM%d",ipt,icell,iasym,iSM) ;
542 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for pt >%2.2f, ncell>%d and asym >%1.2f, SM %d ",
543 fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym],iSM) ;
544 fhRePtNCellAsymCutsSM[iSM][index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
545 fhRePtNCellAsymCutsSM[iSM][index]->SetXTitle("p_{T} (GeV/c)");
546 fhRePtNCellAsymCutsSM[iSM][index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
547 outputContainer->Add(fhRePtNCellAsymCutsSM[iSM][index]) ;
557 fhRePtMult = new TH3F*[fNAsymCuts] ;
558 for(Int_t iasym = 0; iasym<fNAsymCuts; iasym++)
560 fhRePtMult[iasym] = new TH3F(Form("hRePtMult_asym%d",iasym),Form("(p_{T},C,M)_{#gamma#gamma}, A<%1.2f",fAsymCuts[iasym]),
561 nptbins,ptmin,ptmax,ntrmbins,ntrmmin,ntrmmax,nmassbins,massmin,massmax);
562 fhRePtMult[iasym]->SetXTitle("p_{T} (GeV/c)");
563 fhRePtMult[iasym]->SetYTitle("Track multiplicity");
564 fhRePtMult[iasym]->SetZTitle("m_{#gamma,#gamma} (GeV/c^{2})");
565 outputContainer->Add(fhRePtMult[iasym]) ;
568 }// multi cuts analysis
570 if(fFillSSCombinations)
573 fhReSS[0] = new TH2F("hRe_SS_Tight"," 0.01 < #lambda_{0}^{2} < 0.4",
574 nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
575 fhReSS[0]->SetXTitle("p_{T} (GeV/c)");
576 fhReSS[0]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
577 outputContainer->Add(fhReSS[0]) ;
580 fhReSS[1] = new TH2F("hRe_SS_Loose"," #lambda_{0}^{2} > 0.4",
581 nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
582 fhReSS[1]->SetXTitle("p_{T} (GeV/c)");
583 fhReSS[1]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
584 outputContainer->Add(fhReSS[1]) ;
587 fhReSS[2] = new TH2F("hRe_SS_Both"," cluster_{1} #lambda_{0}^{2} > 0.4; cluster_{2} 0.01 < #lambda_{0}^{2} < 0.4",
588 nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
589 fhReSS[2]->SetXTitle("p_{T} (GeV/c)");
590 fhReSS[2]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
591 outputContainer->Add(fhReSS[2]) ;
594 fhEventBin=new TH1I("hEventBin","Number of real pairs per bin(cen,vz,rp)",
595 GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1,0,
596 GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1) ;
597 fhEventBin->SetXTitle("bin");
598 outputContainer->Add(fhEventBin) ;
600 fhEventMixBin=new TH1I("hEventMixBin","Number of mixed pairs per bin(cen,vz,rp)",
601 GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1,0,
602 GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1) ;
603 fhEventMixBin->SetXTitle("bin");
604 outputContainer->Add(fhEventMixBin) ;
608 fhCentrality=new TH1F("hCentralityBin","Number of events in centrality bin",GetNCentrBin(),0.,1.*GetNCentrBin()) ;
609 fhCentrality->SetXTitle("Centrality bin");
610 outputContainer->Add(fhCentrality) ;
612 fhCentralityNoPair=new TH1F("hCentralityBinNoPair","Number of events in centrality bin, with no cluster pairs",GetNCentrBin(),0.,1.*GetNCentrBin()) ;
613 fhCentralityNoPair->SetXTitle("Centrality bin");
614 outputContainer->Add(fhCentralityNoPair) ;
617 if(GetNRPBin() > 1 && GetNCentrBin()>1 )
619 fhEventPlaneResolution=new TH2F("hEventPlaneResolution","Event plane resolution",GetNCentrBin(),0,GetNCentrBin(),100,0.,TMath::TwoPi()) ;
620 fhEventPlaneResolution->SetYTitle("Resolution");
621 fhEventPlaneResolution->SetXTitle("Centrality Bin");
622 outputContainer->Add(fhEventPlaneResolution) ;
627 fhRealOpeningAngle = new TH2F
628 ("hRealOpeningAngle","Angle between all #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,300,0,TMath::Pi());
629 fhRealOpeningAngle->SetYTitle("#theta(rad)");
630 fhRealOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
631 outputContainer->Add(fhRealOpeningAngle) ;
633 fhRealCosOpeningAngle = new TH2F
634 ("hRealCosOpeningAngle","Cosinus of angle between all #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,100,0,1);
635 fhRealCosOpeningAngle->SetYTitle("cos (#theta) ");
636 fhRealCosOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
637 outputContainer->Add(fhRealCosOpeningAngle) ;
641 fhMixedOpeningAngle = new TH2F
642 ("hMixedOpeningAngle","Angle between all #gamma pair vs E_{#pi^{0}}, Mixed pairs",nptbins,ptmin,ptmax,300,0,TMath::Pi());
643 fhMixedOpeningAngle->SetYTitle("#theta(rad)");
644 fhMixedOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
645 outputContainer->Add(fhMixedOpeningAngle) ;
647 fhMixedCosOpeningAngle = new TH2F
648 ("hMixedCosOpeningAngle","Cosinus of angle between all #gamma pair vs E_{#pi^{0}}, Mixed pairs",nptbins,ptmin,ptmax,100,0,1);
649 fhMixedCosOpeningAngle->SetYTitle("cos (#theta) ");
650 fhMixedCosOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
651 outputContainer->Add(fhMixedCosOpeningAngle) ;
656 //Histograms filled only if MC data is requested
659 fhReMCFromConversion = new TH2F("hReMCFromConversion","Invariant mass of 2 clusters originated in conversions",
660 nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
661 fhReMCFromConversion->SetXTitle("p_{T} (GeV/c)");
662 fhReMCFromConversion->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
663 outputContainer->Add(fhReMCFromConversion) ;
665 fhReMCFromNotConversion = new TH2F("hReMCNotFromConversion","Invariant mass of 2 clusters not originated in conversions",
666 nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
667 fhReMCFromNotConversion->SetXTitle("p_{T} (GeV/c)");
668 fhReMCFromNotConversion->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
669 outputContainer->Add(fhReMCFromNotConversion) ;
671 fhReMCFromMixConversion = new TH2F("hReMCFromMixConversion","Invariant mass of 2 clusters one from conversion and the other not",
672 nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
673 fhReMCFromMixConversion->SetXTitle("p_{T} (GeV/c)");
674 fhReMCFromMixConversion->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
675 outputContainer->Add(fhReMCFromMixConversion) ;
679 fhPrimPi0E = new TH1F("hPrimPi0E","Primary pi0 E, Y<1",nptbins,ptmin,ptmax) ;
680 fhPrimPi0AccE = new TH1F("hPrimPi0AccE","Primary pi0 E with both photons in acceptance",nptbins,ptmin,ptmax) ;
681 fhPrimPi0E ->SetXTitle("E (GeV)");
682 fhPrimPi0AccE->SetXTitle("E (GeV)");
683 outputContainer->Add(fhPrimPi0E) ;
684 outputContainer->Add(fhPrimPi0AccE) ;
686 fhPrimPi0Pt = new TH1F("hPrimPi0Pt","Primary pi0 pt, Y<1",nptbins,ptmin,ptmax) ;
687 fhPrimPi0AccPt = new TH1F("hPrimPi0AccPt","Primary pi0 pt with both photons in acceptance",nptbins,ptmin,ptmax) ;
688 fhPrimPi0Pt ->SetXTitle("p_{T} (GeV/c)");
689 fhPrimPi0AccPt->SetXTitle("p_{T} (GeV/c)");
690 outputContainer->Add(fhPrimPi0Pt) ;
691 outputContainer->Add(fhPrimPi0AccPt) ;
693 Int_t netabinsopen = TMath::Nint(netabins*4/(etamax-etamin));
694 fhPrimPi0Y = new TH2F("hPrimPi0Rapidity","Rapidity of primary pi0",nptbins,ptmin,ptmax,netabinsopen,-2, 2) ;
695 fhPrimPi0Y ->SetYTitle("Rapidity");
696 fhPrimPi0Y ->SetXTitle("p_{T} (GeV/c)");
697 outputContainer->Add(fhPrimPi0Y) ;
699 fhPrimPi0AccY = new TH2F("hPrimPi0AccRapidity","Rapidity of primary pi0",nptbins,ptmin,ptmax,netabins,etamin,etamax) ;
700 fhPrimPi0AccY->SetYTitle("Rapidity");
701 fhPrimPi0AccY->SetXTitle("p_{T} (GeV/c)");
702 outputContainer->Add(fhPrimPi0AccY) ;
704 Int_t nphibinsopen = TMath::Nint(nphibins*TMath::TwoPi()/(phimax-phimin));
705 fhPrimPi0Phi = new TH2F("hPrimPi0Phi","Azimuthal of primary pi0, Y<1",nptbins,ptmin,ptmax,nphibinsopen,0,360) ;
706 fhPrimPi0Phi->SetYTitle("#phi (deg)");
707 fhPrimPi0Phi->SetXTitle("p_{T} (GeV/c)");
708 outputContainer->Add(fhPrimPi0Phi) ;
710 fhPrimPi0AccPhi = new TH2F("hPrimPi0AccPhi","Azimuthal of primary pi0 with accepted daughters",nptbins,ptmin,ptmax,
711 nphibins,phimin*TMath::RadToDeg(),phimax*TMath::RadToDeg()) ;
712 fhPrimPi0AccPhi->SetYTitle("#phi (deg)");
713 fhPrimPi0AccPhi->SetXTitle("p_{T} (GeV/c)");
714 outputContainer->Add(fhPrimPi0AccPhi) ;
716 fhPrimPi0PtCentrality = new TH2F("hPrimPi0PtCentrality","Primary pi0 pt vs reco centrality, Y<1",nptbins,ptmin,ptmax, 100, 0, 100) ;
717 fhPrimPi0AccPtCentrality = new TH2F("hPrimPi0AccPtCentrality","Primary pi0 with both photons in acceptance pt vs reco centrality",nptbins,ptmin,ptmax, 100, 0, 100) ;
718 fhPrimPi0PtCentrality ->SetXTitle("p_{T} (GeV/c)");
719 fhPrimPi0AccPtCentrality->SetXTitle("p_{T} (GeV/c)");
720 fhPrimPi0PtCentrality ->SetYTitle("Centrality");
721 fhPrimPi0AccPtCentrality->SetYTitle("Centrality");
722 outputContainer->Add(fhPrimPi0PtCentrality) ;
723 outputContainer->Add(fhPrimPi0AccPtCentrality) ;
725 fhPrimPi0PtEventPlane = new TH2F("hPrimPi0PtEventPlane","Primary pi0 pt vs reco event plane angle, Y<1",nptbins,ptmin,ptmax, 100, 0, TMath::Pi()) ;
726 fhPrimPi0AccPtEventPlane = new TH2F("hPrimPi0AccPtEventPlane","Primary pi0 with both photons in acceptance pt vs reco event plane angle",nptbins,ptmin,ptmax, 100, 0, TMath::Pi()) ;
727 fhPrimPi0PtEventPlane ->SetXTitle("p_{T} (GeV/c)");
728 fhPrimPi0AccPtEventPlane->SetXTitle("p_{T} (GeV/c)");
729 fhPrimPi0PtEventPlane ->SetYTitle("Event Plane Angle (rad)");
730 fhPrimPi0AccPtEventPlane->SetYTitle("Event Plane Angle (rad)");
731 outputContainer->Add(fhPrimPi0PtEventPlane) ;
732 outputContainer->Add(fhPrimPi0AccPtEventPlane) ;
736 fhPrimEtaE = new TH1F("hPrimEtaE","Primary eta E",nptbins,ptmin,ptmax) ;
737 fhPrimEtaAccE = new TH1F("hPrimEtaAccE","Primary eta E with both photons in acceptance",nptbins,ptmin,ptmax) ;
738 fhPrimEtaE ->SetXTitle("E (GeV)");
739 fhPrimEtaAccE->SetXTitle("E (GeV)");
740 outputContainer->Add(fhPrimEtaE) ;
741 outputContainer->Add(fhPrimEtaAccE) ;
743 fhPrimEtaPt = new TH1F("hPrimEtaPt","Primary eta pt",nptbins,ptmin,ptmax) ;
744 fhPrimEtaAccPt = new TH1F("hPrimEtaAccPt","Primary eta pt with both photons in acceptance",nptbins,ptmin,ptmax) ;
745 fhPrimEtaPt ->SetXTitle("p_{T} (GeV/c)");
746 fhPrimEtaAccPt->SetXTitle("p_{T} (GeV/c)");
747 outputContainer->Add(fhPrimEtaPt) ;
748 outputContainer->Add(fhPrimEtaAccPt) ;
750 fhPrimEtaY = new TH2F("hPrimEtaRapidity","Rapidity of primary eta",nptbins,ptmin,ptmax,netabins,etamin,etamax) ;
751 fhPrimEtaY->SetYTitle("Rapidity");
752 fhPrimEtaY->SetXTitle("p_{T} (GeV/c)");
753 outputContainer->Add(fhPrimEtaY) ;
755 fhPrimEtaAccY = new TH2F("hPrimEtaAccRapidity","Rapidity of primary eta",nptbins,ptmin,ptmax, netabins,etamin,etamax) ;
756 fhPrimEtaAccY->SetYTitle("Rapidity");
757 fhPrimEtaAccY->SetXTitle("p_{T} (GeV/c)");
758 outputContainer->Add(fhPrimEtaAccY) ;
760 fhPrimEtaPhi = new TH2F("hPrimEtaPhi","Azimuthal of primary eta",nptbins,ptmin,ptmax, nphibins,phimin*TMath::RadToDeg(),phimax*TMath::RadToDeg()) ;
761 fhPrimEtaPhi->SetYTitle("#phi (deg)");
762 fhPrimEtaPhi->SetXTitle("p_{T} (GeV/c)");
763 outputContainer->Add(fhPrimEtaPhi) ;
765 fhPrimEtaAccPhi = new TH2F("hPrimEtaAccPhi","Azimuthal of primary eta with accepted daughters",nptbins,ptmin,ptmax, nphibins,phimin*TMath::RadToDeg(),phimax*TMath::RadToDeg()) ;
766 fhPrimEtaAccPhi->SetYTitle("#phi (deg)");
767 fhPrimEtaAccPhi->SetXTitle("p_{T} (GeV/c)");
768 outputContainer->Add(fhPrimEtaAccPhi) ;
770 fhPrimEtaPtCentrality = new TH2F("hPrimEtaPtCentrality","Primary eta pt vs reco centrality, Y<1",nptbins,ptmin,ptmax, 100, 0, 100) ;
771 fhPrimEtaAccPtCentrality = new TH2F("hPrimEtaAccPtCentrality","Primary eta with both photons in acceptance pt vs reco centrality",nptbins,ptmin,ptmax, 100, 0, 100) ;
772 fhPrimEtaPtCentrality ->SetXTitle("p_{T} (GeV/c)");
773 fhPrimEtaAccPtCentrality->SetXTitle("p_{T} (GeV/c)");
774 fhPrimEtaPtCentrality ->SetYTitle("Centrality");
775 fhPrimEtaAccPtCentrality->SetYTitle("Centrality");
776 outputContainer->Add(fhPrimEtaPtCentrality) ;
777 outputContainer->Add(fhPrimEtaAccPtCentrality) ;
779 fhPrimEtaPtEventPlane = new TH2F("hPrimEtaPtEventPlane","Primary eta pt vs reco event plane angle, Y<1",nptbins,ptmin,ptmax, 100, 0, TMath::Pi()) ;
780 fhPrimEtaAccPtEventPlane = new TH2F("hPrimEtaAccPtEventPlane","Primary eta with both photons in acceptance pt vs reco event plane angle",nptbins,ptmin,ptmax, 100, 0, TMath::Pi()) ;
781 fhPrimEtaPtEventPlane ->SetXTitle("p_{T} (GeV/c)");
782 fhPrimEtaAccPtEventPlane->SetXTitle("p_{T} (GeV/c)");
783 fhPrimEtaPtEventPlane ->SetYTitle("Event Plane Angle (rad)");
784 fhPrimEtaAccPtEventPlane->SetYTitle("Event Plane Angle (rad)");
785 outputContainer->Add(fhPrimEtaPtEventPlane) ;
786 outputContainer->Add(fhPrimEtaAccPtEventPlane) ;
790 fhPrimPi0OpeningAngle = new TH2F
791 ("hPrimPi0OpeningAngle","Angle between all primary #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,100,0,0.5);
792 fhPrimPi0OpeningAngle->SetYTitle("#theta(rad)");
793 fhPrimPi0OpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
794 outputContainer->Add(fhPrimPi0OpeningAngle) ;
796 fhPrimPi0OpeningAngleAsym = new TH2F
797 ("hPrimPi0OpeningAngleAsym","Angle between all primary #gamma pair vs Asymmetry, p_{T}>5 GeV/c",100,0,1,100,0,0.5);
798 fhPrimPi0OpeningAngleAsym->SetXTitle("|A|=| (E_{1}-E_{2}) / (E_{1}+E_{2}) |");
799 fhPrimPi0OpeningAngleAsym->SetYTitle("#theta(rad)");
800 outputContainer->Add(fhPrimPi0OpeningAngleAsym) ;
802 fhPrimPi0CosOpeningAngle = new TH2F
803 ("hPrimPi0CosOpeningAngle","Cosinus of angle between all primary #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,100,-1,1);
804 fhPrimPi0CosOpeningAngle->SetYTitle("cos (#theta) ");
805 fhPrimPi0CosOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
806 outputContainer->Add(fhPrimPi0CosOpeningAngle) ;
808 fhPrimEtaOpeningAngle = new TH2F
809 ("hPrimEtaOpeningAngle","Angle between all primary #gamma pair vs E_{#eta}",nptbins,ptmin,ptmax,100,0,0.5);
810 fhPrimEtaOpeningAngle->SetYTitle("#theta(rad)");
811 fhPrimEtaOpeningAngle->SetXTitle("E_{#eta} (GeV)");
812 outputContainer->Add(fhPrimEtaOpeningAngle) ;
814 fhPrimEtaOpeningAngleAsym = new TH2F
815 ("hPrimEtaOpeningAngleAsym","Angle between all primary #gamma pair vs Asymmetry, p_{T}>5 GeV/c",100,0,1,100,0,0.5);
816 fhPrimEtaOpeningAngleAsym->SetXTitle("|A|=| (E_{1}-E_{2}) / (E_{1}+E_{2}) |");
817 fhPrimEtaOpeningAngleAsym->SetYTitle("#theta(rad)");
818 outputContainer->Add(fhPrimEtaOpeningAngleAsym) ;
821 fhPrimEtaCosOpeningAngle = new TH2F
822 ("hPrimEtaCosOpeningAngle","Cosinus of angle between all primary #gamma pair vs E_{#eta}",nptbins,ptmin,ptmax,100,-1,1);
823 fhPrimEtaCosOpeningAngle->SetYTitle("cos (#theta) ");
824 fhPrimEtaCosOpeningAngle->SetXTitle("E_{ #eta} (GeV)");
825 outputContainer->Add(fhPrimEtaCosOpeningAngle) ;
834 fhPrimPi0PtOrigin = new TH2F("hPrimPi0PtOrigin","Primary pi0 pt vs origin",nptbins,ptmin,ptmax,11,0,11) ;
835 fhPrimPi0PtOrigin->SetXTitle("p_{T} (GeV/c)");
836 fhPrimPi0PtOrigin->SetYTitle("Origin");
837 fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(1 ,"Status 21");
838 fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(2 ,"Quark");
839 fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(3 ,"qq Resonances ");
840 fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(4 ,"Resonances");
841 fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(5 ,"#rho");
842 fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(6 ,"#omega");
843 fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(7 ,"K");
844 fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(8 ,"Other");
845 fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(9 ,"#eta");
846 fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(10 ,"#eta prime");
847 outputContainer->Add(fhPrimPi0PtOrigin) ;
849 fhMCPi0PtOrigin = new TH2F("hMCPi0PtOrigin","Reconstructed pair from generated pi0 pt vs origin",nptbins,ptmin,ptmax,11,0,11) ;
850 fhMCPi0PtOrigin->SetXTitle("p_{T} (GeV/c)");
851 fhMCPi0PtOrigin->SetYTitle("Origin");
852 fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(1 ,"Status 21");
853 fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(2 ,"Quark");
854 fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(3 ,"qq Resonances");
855 fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(4 ,"Resonances");
856 fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(5 ,"#rho");
857 fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(6 ,"#omega");
858 fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(7 ,"K");
859 fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(8 ,"Other");
860 fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(9 ,"#eta");
861 fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(10 ,"#eta prime");
862 outputContainer->Add(fhMCPi0PtOrigin) ;
865 fhPrimEtaPtOrigin = new TH2F("hPrimEtaPtOrigin","Primary pi0 pt vs origin",nptbins,ptmin,ptmax,7,0,7) ;
866 fhPrimEtaPtOrigin->SetXTitle("p_{T} (GeV/c)");
867 fhPrimEtaPtOrigin->SetYTitle("Origin");
868 fhPrimEtaPtOrigin->GetYaxis()->SetBinLabel(1 ,"Status 21");
869 fhPrimEtaPtOrigin->GetYaxis()->SetBinLabel(2 ,"Quark");
870 fhPrimEtaPtOrigin->GetYaxis()->SetBinLabel(3 ,"qq Resonances");
871 fhPrimEtaPtOrigin->GetYaxis()->SetBinLabel(4 ,"Resonances");
872 fhPrimEtaPtOrigin->GetYaxis()->SetBinLabel(5 ,"Other");
873 fhPrimEtaPtOrigin->GetYaxis()->SetBinLabel(6 ,"#eta prime ");
875 outputContainer->Add(fhPrimEtaPtOrigin) ;
877 fhMCEtaPtOrigin = new TH2F("hMCEtaPtOrigin","Reconstructed pair from generated pi0 pt vs origin",nptbins,ptmin,ptmax,7,0,7) ;
878 fhMCEtaPtOrigin->SetXTitle("p_{T} (GeV/c)");
879 fhMCEtaPtOrigin->SetYTitle("Origin");
880 fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(1 ,"Status 21");
881 fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(2 ,"Quark");
882 fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(3 ,"qq Resonances");
883 fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(4 ,"Resonances");
884 fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(5 ,"Other");
885 fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(6 ,"#eta prime");
887 outputContainer->Add(fhMCEtaPtOrigin) ;
889 for(Int_t i = 0; i<13; i++)
891 fhMCOrgMass[i] = new TH2F(Form("hMCOrgMass_%d",i),Form("mass vs pt, origin %d",i),nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
892 fhMCOrgMass[i]->SetXTitle("p_{T} (GeV/c)");
893 fhMCOrgMass[i]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
894 outputContainer->Add(fhMCOrgMass[i]) ;
896 fhMCOrgAsym[i]= new TH2F(Form("hMCOrgAsym_%d",i),Form("asymmetry vs pt, origin %d",i),nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax) ;
897 fhMCOrgAsym[i]->SetXTitle("p_{T} (GeV/c)");
898 fhMCOrgAsym[i]->SetYTitle("A");
899 outputContainer->Add(fhMCOrgAsym[i]) ;
901 fhMCOrgDeltaEta[i] = new TH2F(Form("hMCOrgDeltaEta_%d",i),Form("#Delta #eta of pair vs pt, origin %d",i),nptbins,ptmin,ptmax,netabins,-1.4,1.4) ;
902 fhMCOrgDeltaEta[i]->SetXTitle("p_{T} (GeV/c)");
903 fhMCOrgDeltaEta[i]->SetYTitle("#Delta #eta");
904 outputContainer->Add(fhMCOrgDeltaEta[i]) ;
906 fhMCOrgDeltaPhi[i]= new TH2F(Form("hMCOrgDeltaPhi_%d",i),Form("#Delta #phi of pair vs p_{T}, origin %d",i),nptbins,ptmin,ptmax,nphibins,-0.7,0.7) ;
907 fhMCOrgDeltaPhi[i]->SetXTitle("p_{T} (GeV/c)");
908 fhMCOrgDeltaPhi[i]->SetYTitle("#Delta #phi (rad)");
909 outputContainer->Add(fhMCOrgDeltaPhi[i]) ;
915 fhMCPi0MassPtTrue = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
916 fhMCPi0MassPtRec = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
917 fhMCPi0PtTruePtRec = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
918 fhMCEtaMassPtRec = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
919 fhMCEtaMassPtTrue = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
920 fhMCEtaPtTruePtRec = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
921 for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
922 for(Int_t icell=0; icell<fNCellNCuts; icell++){
923 for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
924 Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
926 fhMCPi0MassPtRec[index] = new TH2F(Form("hMCPi0MassPtRec_pt%d_cell%d_asym%d",ipt,icell,iasym),
927 Form("Reconstructed Mass vs reconstructed p_T of true #pi^{0} cluster pairs for pt >%2.2f, ncell>%d and asym >%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]),
928 nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
929 fhMCPi0MassPtRec[index]->SetXTitle("p_{T, reconstructed} (GeV/c)");
930 fhMCPi0MassPtRec[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
931 outputContainer->Add(fhMCPi0MassPtRec[index]) ;
933 fhMCPi0MassPtTrue[index] = new TH2F(Form("hMCPi0MassPtTrue_pt%d_cell%d_asym%d",ipt,icell,iasym),
934 Form("Reconstructed Mass vs generated p_T of true #pi^{0} cluster pairs for pt >%2.2f, ncell>%d and asym >%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]),
935 nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
936 fhMCPi0MassPtTrue[index]->SetXTitle("p_{T, generated} (GeV/c)");
937 fhMCPi0MassPtTrue[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
938 outputContainer->Add(fhMCPi0MassPtTrue[index]) ;
940 fhMCPi0PtTruePtRec[index] = new TH2F(Form("hMCPi0PtTruePtRec_pt%d_cell%d_asym%d",ipt,icell,iasym),
941 Form("Generated vs reconstructed p_T of true #pi^{0} cluster pairs, 0.01 < rec. mass < 0.17 MeV/c^{2} for pt >%2.2f, ncell>%d and asym >%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]),
942 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax) ;
943 fhMCPi0PtTruePtRec[index]->SetXTitle("p_{T, generated} (GeV/c)");
944 fhMCPi0PtTruePtRec[index]->SetYTitle("p_{T, reconstructed} (GeV/c)");
945 outputContainer->Add(fhMCPi0PtTruePtRec[index]) ;
947 fhMCEtaMassPtRec[index] = new TH2F(Form("hMCEtaMassPtRec_pt%d_cell%d_asym%d",ipt,icell,iasym),
948 Form("Reconstructed Mass vs reconstructed p_T of true #eta cluster pairs for pt >%2.2f, ncell>%d and asym >%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]),
949 nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
950 fhMCEtaMassPtRec[index]->SetXTitle("p_{T, generated} (GeV/c)");
951 fhMCEtaMassPtRec[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
952 outputContainer->Add(fhMCEtaMassPtRec[index]) ;
954 fhMCEtaMassPtTrue[index] = new TH2F(Form("hMCEtaMassPtTrue_pt%d_cell%d_asym%d",ipt,icell,iasym),
955 Form("Reconstructed Mass vs generated p_T of true #eta cluster pairs for pt >%2.2f, ncell>%d and asym >%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]),
956 nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
957 fhMCEtaMassPtTrue[index]->SetXTitle("p_{T, generated} (GeV/c)");
958 fhMCEtaMassPtTrue[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
959 outputContainer->Add(fhMCEtaMassPtTrue[index]) ;
961 fhMCEtaPtTruePtRec[index] = new TH2F(Form("hMCEtaPtTruePtRec_pt%d_cell%d_asym%d",ipt,icell,iasym),
962 Form("Generated vs reconstructed p_T of true #eta cluster pairs, 0.01 < rec. mass < 0.17 MeV/c^{2} for pt >%2.2f, ncell>%d and asym >%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]),
963 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax) ;
964 fhMCEtaPtTruePtRec[index]->SetXTitle("p_{T, generated} (GeV/c)");
965 fhMCEtaPtTruePtRec[index]->SetYTitle("p_{T, reconstructed} (GeV/c)");
966 outputContainer->Add(fhMCEtaPtTruePtRec[index]) ;
973 fhMCPi0MassPtTrue = new TH2F*[1];
974 fhMCPi0PtTruePtRec = new TH2F*[1];
975 fhMCEtaMassPtTrue = new TH2F*[1];
976 fhMCEtaPtTruePtRec = new TH2F*[1];
978 fhMCPi0MassPtTrue[0] = new TH2F("hMCPi0MassPtTrue","Reconstructed Mass vs generated p_T of true #pi^{0} cluster pairs",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
979 fhMCPi0MassPtTrue[0]->SetXTitle("p_{T, generated} (GeV/c)");
980 fhMCPi0MassPtTrue[0]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
981 outputContainer->Add(fhMCPi0MassPtTrue[0]) ;
983 fhMCPi0PtTruePtRec[0]= new TH2F("hMCPi0PtTruePtRec","Generated vs reconstructed p_T of true #pi^{0} cluster pairs, 0.01 < rec. mass < 0.17 MeV/c^{2}",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax) ;
984 fhMCPi0PtTruePtRec[0]->SetXTitle("p_{T, generated} (GeV/c)");
985 fhMCPi0PtTruePtRec[0]->SetYTitle("p_{T, reconstructed} (GeV/c)");
986 outputContainer->Add(fhMCPi0PtTruePtRec[0]) ;
988 fhMCEtaMassPtTrue[0] = new TH2F("hMCEtaMassPtTrue","Reconstructed Mass vs generated p_T of true #eta cluster pairs",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
989 fhMCEtaMassPtTrue[0]->SetXTitle("p_{T, generated} (GeV/c)");
990 fhMCEtaMassPtTrue[0]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
991 outputContainer->Add(fhMCEtaMassPtTrue[0]) ;
993 fhMCEtaPtTruePtRec[0]= new TH2F("hMCEtaPtTruePtRec","Generated vs reconstructed p_T of true #eta cluster pairs, 0.01 < rec. mass < 0.17 MeV/c^{2}",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax) ;
994 fhMCEtaPtTruePtRec[0]->SetXTitle("p_{T, generated} (GeV/c)");
995 fhMCEtaPtTruePtRec[0]->SetYTitle("p_{T, reconstructed} (GeV/c)");
996 outputContainer->Add(fhMCEtaPtTruePtRec[0]) ;
1001 if(fFillSMCombinations)
1003 TString pairnamePHOS[] = {"(0-1)","(0-2)","(1-2)","(0-3)","(0-4)","(1-3)","(1-4)","(2-3)","(2-4)","(3-4)"};
1004 for(Int_t imod=0; imod<fNModules; imod++)
1006 //Module dependent invariant mass
1007 snprintf(key, buffersize,"hReMod_%d",imod) ;
1008 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for Module %d",imod) ;
1009 fhReMod[imod] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1010 fhReMod[imod]->SetXTitle("p_{T} (GeV/c)");
1011 fhReMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
1012 outputContainer->Add(fhReMod[imod]) ;
1013 if(fCalorimeter=="PHOS")
1015 snprintf(key, buffersize,"hReDiffPHOSMod_%d",imod) ;
1016 snprintf(title, buffersize,"Real pairs PHOS, clusters in different Modules: %s",(pairnamePHOS[imod]).Data()) ;
1017 fhReDiffPHOSMod[imod] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1018 fhReDiffPHOSMod[imod]->SetXTitle("p_{T} (GeV/c)");
1019 fhReDiffPHOSMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
1020 outputContainer->Add(fhReDiffPHOSMod[imod]) ;
1023 if(imod<fNModules/2)
1025 snprintf(key, buffersize,"hReSameSectorEMCAL_%d",imod) ;
1026 snprintf(title, buffersize,"Real pairs EMCAL, clusters in same sector, SM(%d,%d)",imod*2,imod*2+1) ;
1027 fhReSameSectorEMCALMod[imod] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1028 fhReSameSectorEMCALMod[imod]->SetXTitle("p_{T} (GeV/c)");
1029 fhReSameSectorEMCALMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
1030 outputContainer->Add(fhReSameSectorEMCALMod[imod]) ;
1032 if(imod<fNModules-2)
1034 snprintf(key, buffersize,"hReSameSideEMCAL_%d",imod) ;
1035 snprintf(title, buffersize,"Real pairs EMCAL, clusters in same side SM(%d,%d)",imod, imod+2) ;
1036 fhReSameSideEMCALMod[imod] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1037 fhReSameSideEMCALMod[imod]->SetXTitle("p_{T} (GeV/c)");
1038 fhReSameSideEMCALMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
1039 outputContainer->Add(fhReSameSideEMCALMod[imod]) ;
1045 snprintf(key, buffersize,"hMiMod_%d",imod) ;
1046 snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for Module %d",imod) ;
1047 fhMiMod[imod] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1048 fhMiMod[imod]->SetXTitle("p_{T} (GeV/c)");
1049 fhMiMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
1050 outputContainer->Add(fhMiMod[imod]) ;
1052 if(fCalorimeter=="PHOS"){
1053 snprintf(key, buffersize,"hMiDiffPHOSMod_%d",imod) ;
1054 snprintf(title, buffersize,"Mixed pairs PHOS, clusters in different Modules: %s",(pairnamePHOS[imod]).Data()) ;
1055 fhMiDiffPHOSMod[imod] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1056 fhMiDiffPHOSMod[imod]->SetXTitle("p_{T} (GeV/c)");
1057 fhMiDiffPHOSMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
1058 outputContainer->Add(fhMiDiffPHOSMod[imod]) ;
1061 if(imod<fNModules/2)
1063 snprintf(key, buffersize,"hMiSameSectorEMCALMod_%d",imod) ;
1064 snprintf(title, buffersize,"Mixed pairs EMCAL, clusters in same sector, SM(%d,%d)",imod*2,imod*2+1) ;
1065 fhMiSameSectorEMCALMod[imod] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1066 fhMiSameSectorEMCALMod[imod]->SetXTitle("p_{T} (GeV/c)");
1067 fhMiSameSectorEMCALMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
1068 outputContainer->Add(fhMiSameSectorEMCALMod[imod]) ;
1070 if(imod<fNModules-2){
1072 snprintf(key, buffersize,"hMiSameSideEMCALMod_%d",imod) ;
1073 snprintf(title, buffersize,"Mixed pairs EMCAL, clusters in same side SM(%d,%d)",imod, imod+2) ;
1074 fhMiSameSideEMCALMod[imod] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1075 fhMiSameSideEMCALMod[imod]->SetXTitle("p_{T} (GeV/c)");
1076 fhMiSameSideEMCALMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
1077 outputContainer->Add(fhMiSameSideEMCALMod[imod]) ;
1081 }//loop combinations
1082 } // SM combinations
1084 if(fFillArmenterosThetaStar && IsDataMC())
1086 TString ebin[] = {"8 < E < 12 GeV","12 < E < 16 GeV", "16 < E < 20 GeV", "E > 20 GeV" };
1087 Int_t narmbins = 400;
1089 Float_t armmax = 0.4;
1091 for(Int_t i = 0; i < 4; i++)
1093 fhArmPrimPi0[i] = new TH2F(Form("hArmenterosPrimPi0EBin%d",i),
1094 Form("Armenteros of primary #pi^{0}, %s",ebin[i].Data()),
1095 200, -1, 1, narmbins,armmin,armmax);
1096 fhArmPrimPi0[i]->SetYTitle("p_{T}^{Arm}");
1097 fhArmPrimPi0[i]->SetXTitle("#alpha^{Arm}");
1098 outputContainer->Add(fhArmPrimPi0[i]) ;
1100 fhArmPrimEta[i] = new TH2F(Form("hArmenterosPrimEtaEBin%d",i),
1101 Form("Armenteros of primary #eta, %s",ebin[i].Data()),
1102 200, -1, 1, narmbins,armmin,armmax);
1103 fhArmPrimEta[i]->SetYTitle("p_{T}^{Arm}");
1104 fhArmPrimEta[i]->SetXTitle("#alpha^{Arm}");
1105 outputContainer->Add(fhArmPrimEta[i]) ;
1109 // Same as asymmetry ...
1110 fhCosThStarPrimPi0 = new TH2F
1111 ("hCosThStarPrimPi0","cos(#theta *) for primary #pi^{0}",nptbins,ptmin,ptmax,200,-1,1);
1112 fhCosThStarPrimPi0->SetYTitle("cos(#theta *)");
1113 fhCosThStarPrimPi0->SetXTitle("E_{ #pi^{0}} (GeV)");
1114 outputContainer->Add(fhCosThStarPrimPi0) ;
1116 fhCosThStarPrimEta = new TH2F
1117 ("hCosThStarPrimEta","cos(#theta *) for primary #eta",nptbins,ptmin,ptmax,200,-1,1);
1118 fhCosThStarPrimEta->SetYTitle("cos(#theta *)");
1119 fhCosThStarPrimEta->SetXTitle("E_{ #eta} (GeV)");
1120 outputContainer->Add(fhCosThStarPrimEta) ;
1124 // for(Int_t i = 0; i < outputContainer->GetEntries() ; i++){
1126 // printf("Histogram %d, name: %s\n ",i, outputContainer->At(i)->GetName());
1130 return outputContainer;
1133 //___________________________________________________
1134 void AliAnaPi0::Print(const Option_t * /*opt*/) const
1136 //Print some relevant parameters set for the analysis
1137 printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
1138 AliAnaCaloTrackCorrBaseClass::Print(" ");
1140 printf("Number of bins in Centrality: %d \n",GetNCentrBin()) ;
1141 printf("Number of bins in Z vert. pos: %d \n",GetNZvertBin()) ;
1142 printf("Number of bins in Reac. Plain: %d \n",GetNRPBin()) ;
1143 printf("Depth of event buffer: %d \n",GetNMaxEvMix()) ;
1144 printf("Pair in same Module: %d \n",fSameSM) ;
1145 printf("Cuts: \n") ;
1146 // printf("Z vertex position: -%2.3f < z < %2.3f \n",GetZvertexCut(),GetZvertexCut()) ; //It crashes here, why?
1147 printf("Number of modules: %d \n",fNModules) ;
1148 printf("Select pairs with their angle: %d, edep %d, min angle %2.3f, max angle %2.3f \n",fUseAngleCut, fUseAngleEDepCut, fAngleCut, fAngleMaxCut) ;
1149 printf("Asymmetry cuts: n = %d, \n",fNAsymCuts) ;
1150 printf("\tasymmetry < ");
1151 for(Int_t i = 0; i < fNAsymCuts; i++) printf("%2.2f ",fAsymCuts[i]);
1154 printf("PID selection bits: n = %d, \n",fNPIDBits) ;
1155 printf("\tPID bit = ");
1156 for(Int_t i = 0; i < fNPIDBits; i++) printf("%d ",fPIDBits[i]);
1160 printf("pT cuts: n = %d, \n",fNPtCuts) ;
1162 for(Int_t i = 0; i < fNPtCuts; i++) printf("%2.2f ",fPtCuts[i]);
1165 printf("N cell in cluster cuts: n = %d, \n",fNCellNCuts) ;
1166 printf("\tnCell > ");
1167 for(Int_t i = 0; i < fNCellNCuts; i++) printf("%d ",fCellNCuts[i]);
1171 printf("------------------------------------------------------\n") ;
1174 //________________________________________
1175 void AliAnaPi0::FillAcceptanceHistograms()
1177 //Fill acceptance histograms if MC data is available
1179 Float_t cen = GetEventCentrality();
1180 Float_t ep = GetEventPlaneAngle();
1182 if(GetReader()->ReadStack())
1184 AliStack * stack = GetMCStack();
1187 for(Int_t i=0 ; i<stack->GetNtrack(); i++)
1189 if(GetReader()->AcceptOnlyHIJINGLabels() && !GetReader()->IsHIJINGLabel(i)) continue ;
1191 TParticle * prim = stack->Particle(i) ;
1192 Int_t pdg = prim->GetPdgCode();
1193 //printf("i %d, %s %d %s %d \n",i, stack->Particle(i)->GetName(), stack->Particle(i)->GetPdgCode(),
1194 // prim->GetName(), prim->GetPdgCode());
1196 if( pdg == 111 || pdg == 221)
1198 Double_t pi0Pt = prim->Pt() ;
1199 Double_t pi0E = prim->Energy() ;
1200 if(pi0E == TMath::Abs(prim->Pz())) continue ; //Protection against floating point exception
1201 Double_t pi0Y = 0.5*TMath::Log((pi0E-prim->Pz())/(pi0E+prim->Pz())) ;
1202 Double_t phi = TMath::RadToDeg()*prim->Phi() ;
1205 if(TMath::Abs(pi0Y) < 1.0)
1207 fhPrimPi0E ->Fill(pi0E ) ;
1208 fhPrimPi0Pt ->Fill(pi0Pt) ;
1209 fhPrimPi0Phi->Fill(pi0Pt, phi) ;
1210 fhPrimPi0PtCentrality->Fill(pi0Pt,cen) ;
1211 fhPrimPi0PtEventPlane->Fill(pi0Pt,ep ) ;
1213 fhPrimPi0Y ->Fill(pi0Pt, pi0Y) ;
1217 if(TMath::Abs(pi0Y) < 1.0)
1219 fhPrimEtaE ->Fill(pi0E ) ;
1220 fhPrimEtaPt ->Fill(pi0Pt) ;
1221 fhPrimEtaPhi->Fill(pi0Pt, phi) ;
1222 fhPrimEtaPtCentrality->Fill(pi0Pt,cen) ;
1223 fhPrimEtaPtEventPlane->Fill(pi0Pt,ep ) ;
1225 fhPrimEtaY ->Fill(pi0Pt, pi0Y) ;
1229 if(fFillOriginHisto)
1231 Int_t momindex = prim->GetFirstMother();
1234 TParticle* mother = stack->Particle(momindex);
1235 Int_t mompdg = TMath::Abs(mother->GetPdgCode());
1236 Int_t momstatus = mother->GetStatusCode();
1239 if (momstatus == 21)fhPrimPi0PtOrigin->Fill(pi0Pt,0.5);//parton
1240 else if(mompdg < 22 ) fhPrimPi0PtOrigin->Fill(pi0Pt,1.5);//quark
1241 else if(mompdg > 2100 && mompdg < 2210) fhPrimPi0PtOrigin->Fill(pi0Pt,2.5);// resonances
1242 else if(mompdg == 221) fhPrimPi0PtOrigin->Fill(pi0Pt,8.5);//eta
1243 else if(mompdg == 331) fhPrimPi0PtOrigin->Fill(pi0Pt,9.5);//eta prime
1244 else if(mompdg == 213) fhPrimPi0PtOrigin->Fill(pi0Pt,4.5);//rho
1245 else if(mompdg == 223) fhPrimPi0PtOrigin->Fill(pi0Pt,5.5);//omega
1246 else if(mompdg >= 310 && mompdg <= 323) fhPrimPi0PtOrigin->Fill(pi0Pt,6.5);//k0S, k+-,k*
1247 else if(mompdg == 130) fhPrimPi0PtOrigin->Fill(pi0Pt,6.5);//k0L
1248 else if(momstatus == 11 || momstatus == 12 ) fhPrimPi0PtOrigin->Fill(pi0Pt,3.5);//resonances
1249 else fhPrimPi0PtOrigin->Fill(pi0Pt,7.5);//other?
1253 if (momstatus == 21 ) fhPrimEtaPtOrigin->Fill(pi0Pt,0.5);//parton
1254 else if(mompdg < 22 ) fhPrimEtaPtOrigin->Fill(pi0Pt,1.5);//quark
1255 else if(mompdg > 2100 && mompdg < 2210) fhPrimEtaPtOrigin->Fill(pi0Pt,2.5);//qq resonances
1256 else if(mompdg == 331) fhPrimEtaPtOrigin->Fill(pi0Pt,5.5);//eta prime
1257 else if(momstatus == 11 || momstatus == 12 ) fhPrimEtaPtOrigin->Fill(pi0Pt,3.5);//resonances
1258 else fhPrimEtaPtOrigin->Fill(pi0Pt,4.5);//stable, conversions?
1259 //printf("Other Meson pdg %d, Mother %s, pdg %d, status %d\n",pdg, TDatabasePDG::Instance()->GetParticle(mompdg)->GetName(),mompdg, momstatus );
1264 //Check if both photons hit Calorimeter
1265 if(prim->GetNDaughters()!=2) continue; //Only interested in 2 gamma decay
1266 Int_t iphot1=prim->GetFirstDaughter() ;
1267 Int_t iphot2=prim->GetLastDaughter() ;
1268 if(iphot1>-1 && iphot1<stack->GetNtrack() && iphot2>-1 && iphot2<stack->GetNtrack())
1270 TParticle * phot1 = stack->Particle(iphot1) ;
1271 TParticle * phot2 = stack->Particle(iphot2) ;
1272 if(phot1 && phot2 && phot1->GetPdgCode()==22 && phot2->GetPdgCode()==22)
1274 //printf("2 photons: photon 1: pt %2.2f, phi %3.2f, eta %1.2f; photon 2: pt %2.2f, phi %3.2f, eta %1.2f\n",
1275 // phot1->Pt(), phot1->Phi()*180./3.1415, phot1->Eta(), phot2->Pt(), phot2->Phi()*180./3.1415, phot2->Eta());
1277 TLorentzVector lv1, lv2,lvmeson;
1278 phot1->Momentum(lv1);
1279 phot2->Momentum(lv2);
1280 prim ->Momentum(lvmeson);
1282 if(fFillArmenterosThetaStar) FillArmenterosThetaStar(pdg,lvmeson,lv1,lv2);
1284 Bool_t inacceptance = kFALSE;
1285 if(fCalorimeter == "PHOS")
1287 if(GetPHOSGeometry() && GetCaloUtils()->IsPHOSGeoMatrixSet())
1291 if(GetPHOSGeometry()->ImpactOnEmc(phot1,mod,z,x) && GetPHOSGeometry()->ImpactOnEmc(phot2,mod,z,x))
1292 inacceptance = kTRUE;
1293 if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
1297 if(GetFiducialCut()->IsInFiducialCut(lv1,fCalorimeter) && GetFiducialCut()->IsInFiducialCut(lv2,fCalorimeter))
1298 inacceptance = kTRUE ;
1299 if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
1302 else if(fCalorimeter == "EMCAL" && GetCaloUtils()->IsEMCALGeoMatrixSet())
1304 if(GetEMCALGeometry())
1309 GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(phot1->Eta(),phot1->Phi(),absID1);
1310 GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(phot2->Eta(),phot2->Phi(),absID2);
1312 if( absID1 >= 0 && absID2 >= 0)
1313 inacceptance = kTRUE;
1315 // if(GetEMCALGeometry()->Impact(phot1) && GetEMCALGeometry()->Impact(phot2))
1316 // inacceptance = kTRUE;
1317 if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
1321 if(GetFiducialCut()->IsInFiducialCut(lv1,fCalorimeter) && GetFiducialCut()->IsInFiducialCut(lv2,fCalorimeter))
1322 inacceptance = kTRUE ;
1323 if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
1329 Float_t asym = TMath::Abs((lv1.E()-lv2.E()) / (lv1.E()+lv2.E()));
1330 Double_t angle = lv1.Angle(lv2.Vect());
1334 fhPrimPi0AccE ->Fill(pi0E) ;
1335 fhPrimPi0AccPt ->Fill(pi0Pt) ;
1336 fhPrimPi0AccPhi->Fill(pi0Pt, phi) ;
1337 fhPrimPi0AccY ->Fill(pi0Pt, pi0Y) ;
1338 fhPrimPi0AccPtCentrality->Fill(pi0Pt,cen) ;
1339 fhPrimPi0AccPtEventPlane->Fill(pi0Pt,ep ) ;
1343 fhPrimPi0OpeningAngle ->Fill(pi0Pt,angle);
1344 if(pi0Pt > 5)fhPrimPi0OpeningAngleAsym->Fill(asym,angle);
1345 fhPrimPi0CosOpeningAngle ->Fill(pi0Pt,TMath::Cos(angle));
1350 fhPrimEtaAccE ->Fill(pi0E ) ;
1351 fhPrimEtaAccPt ->Fill(pi0Pt) ;
1352 fhPrimEtaAccPhi->Fill(pi0Pt, phi) ;
1353 fhPrimEtaAccY ->Fill(pi0Pt, pi0Y) ;
1354 fhPrimEtaAccPtCentrality->Fill(pi0Pt,cen) ;
1355 fhPrimEtaAccPtEventPlane->Fill(pi0Pt,ep ) ;
1359 fhPrimEtaOpeningAngle ->Fill(pi0Pt,angle);
1360 if(pi0Pt > 5)fhPrimEtaOpeningAngleAsym->Fill(asym,angle);
1361 fhPrimEtaCosOpeningAngle ->Fill(pi0Pt,TMath::Cos(angle));
1366 }//Check daughters exist
1367 }// Primary pi0 or eta
1368 }//loop on primaries
1369 }//stack exists and data is MC
1371 else if(GetReader()->ReadAODMCParticles())
1373 TClonesArray * mcparticles = GetReader()->GetAODMCParticles();
1376 Int_t nprim = mcparticles->GetEntriesFast();
1378 for(Int_t i=0; i < nprim; i++)
1380 if(GetReader()->AcceptOnlyHIJINGLabels() && !GetReader()->IsHIJINGLabel(i)) continue ;
1382 AliAODMCParticle * prim = (AliAODMCParticle *) mcparticles->At(i);
1384 // Only generator particles, when they come from PYTHIA, PHOJET, HERWIG ...
1385 //if( prim->GetStatus() == 0 && (GetMCAnalysisUtils()->GetMCGenerator()).Length()!=0) break;
1387 Int_t pdg = prim->GetPdgCode();
1388 if( pdg == 111 || pdg == 221)
1390 Double_t pi0Pt = prim->Pt() ;
1391 Double_t pi0E = prim->E() ;
1392 //printf("pi0, pt %2.2f, eta %f, phi %f\n",pi0Pt, prim->Eta(), prim->Phi());
1393 if(pi0E == TMath::Abs(prim->Pz())) continue ; //Protection against floating point exception
1395 Double_t pi0Y = 0.5*TMath::Log((prim->E()-prim->Pz())/(prim->E()+prim->Pz())) ;
1396 Double_t phi = TMath::RadToDeg()*prim->Phi() ;
1399 if(TMath::Abs(pi0Y) < 1)
1401 fhPrimPi0E ->Fill(pi0E ) ;
1402 fhPrimPi0Pt ->Fill(pi0Pt) ;
1403 fhPrimPi0Phi->Fill(pi0Pt, phi) ;
1404 fhPrimPi0PtCentrality->Fill(pi0Pt,cen) ;
1405 fhPrimPi0PtEventPlane->Fill(pi0Pt,ep ) ;
1407 fhPrimPi0Y ->Fill(pi0Pt, pi0Y) ;
1411 if(TMath::Abs(pi0Y) < 1)
1413 fhPrimEtaE ->Fill(pi0E ) ;
1414 fhPrimEtaPt ->Fill(pi0Pt) ;
1415 fhPrimEtaPhi->Fill(pi0Pt, phi) ;
1416 fhPrimEtaPtCentrality->Fill(pi0Pt,cen) ;
1417 fhPrimEtaPtEventPlane->Fill(pi0Pt,ep ) ;
1419 fhPrimEtaY ->Fill(pi0Pt, pi0Y) ;
1423 Int_t momindex = prim->GetMother();
1426 AliAODMCParticle* mother = (AliAODMCParticle *) mcparticles->At(momindex);
1427 Int_t mompdg = TMath::Abs(mother->GetPdgCode());
1428 Int_t momstatus = mother->GetStatus();
1429 if(fFillOriginHisto)
1433 if (momstatus == 21) fhPrimPi0PtOrigin->Fill(pi0Pt,0.5);//parton
1434 else if(mompdg < 22 ) fhPrimPi0PtOrigin->Fill(pi0Pt,1.5);//quark
1435 else if(mompdg > 2100 && mompdg < 2210) fhPrimPi0PtOrigin->Fill(pi0Pt,2.5);// resonances
1436 else if(mompdg == 221) fhPrimPi0PtOrigin->Fill(pi0Pt,8.5);//eta
1437 else if(mompdg == 331) fhPrimPi0PtOrigin->Fill(pi0Pt,9.5);//eta prime
1438 else if(mompdg == 213) fhPrimPi0PtOrigin->Fill(pi0Pt,4.5);//rho
1439 else if(mompdg == 223) fhPrimPi0PtOrigin->Fill(pi0Pt,5.5);//omega
1440 else if(mompdg >= 310 && mompdg <= 323) fhPrimPi0PtOrigin->Fill(pi0Pt,6.5);//k0S, k+-,k*
1441 else if(mompdg == 130) fhPrimPi0PtOrigin->Fill(pi0Pt,6.5);//k0L
1442 else if(momstatus == 11 || momstatus == 12 ) fhPrimPi0PtOrigin->Fill(pi0Pt,3.5);//resonances
1443 else fhPrimPi0PtOrigin->Fill(pi0Pt,7.5);//other?
1447 if (momstatus == 21 ) fhPrimEtaPtOrigin->Fill(pi0Pt,0.5);//parton
1448 else if(mompdg < 22 ) fhPrimEtaPtOrigin->Fill(pi0Pt,1.5);//quark
1449 else if(mompdg > 2100 && mompdg < 2210) fhPrimEtaPtOrigin->Fill(pi0Pt,2.5);//qq resonances
1450 else if(mompdg == 331) fhPrimEtaPtOrigin->Fill(pi0Pt,5.5);//eta prime
1451 else if(momstatus == 11 || momstatus == 12 ) fhPrimEtaPtOrigin->Fill(pi0Pt,3.5);//resonances
1452 else fhPrimEtaPtOrigin->Fill(pi0Pt,4.5);//stable, conversions?
1453 //printf("Other Meson pdg %d, Mother %s, pdg %d, status %d\n",pdg, TDatabasePDG::Instance()->GetParticle(mompdg)->GetName(),mompdg, momstatus );
1458 //Check if both photons hit Calorimeter
1459 if(prim->GetNDaughters()!=2) continue; //Only interested in 2 gamma decay
1460 Int_t iphot1=prim->GetDaughter(0) ;
1461 Int_t iphot2=prim->GetDaughter(1) ;
1462 if(iphot1>-1 && iphot1<nprim && iphot2>-1 && iphot2<nprim)
1464 AliAODMCParticle * phot1 = (AliAODMCParticle *) mcparticles->At(iphot1);
1465 AliAODMCParticle * phot2 = (AliAODMCParticle *) mcparticles->At(iphot2);
1466 if(phot1 && phot2 && phot1->GetPdgCode()==22 && phot2->GetPdgCode()==22)
1468 TLorentzVector lv1, lv2,lvmeson;
1469 lv1.SetPxPyPzE(phot1->Px(),phot1->Py(),phot1->Pz(),phot1->E());
1470 lv2.SetPxPyPzE(phot2->Px(),phot2->Py(),phot2->Pz(),phot2->E());
1471 lvmeson.SetPxPyPzE(prim->Px(),prim->Py(),prim->Pz(),prim->E());
1473 if(fFillArmenterosThetaStar) FillArmenterosThetaStar(pdg,lvmeson,lv1,lv2);
1475 Bool_t inacceptance = kFALSE;
1476 if(fCalorimeter == "PHOS")
1478 if(GetPHOSGeometry() && GetCaloUtils()->IsPHOSGeoMatrixSet())
1482 Double_t vtx []={phot1->Xv(),phot1->Yv(),phot1->Zv()};
1483 Double_t vtx2[]={phot2->Xv(),phot2->Yv(),phot2->Zv()};
1484 if(GetPHOSGeometry()->ImpactOnEmc(vtx, phot1->Theta(),phot1->Phi(),mod,z,x) &&
1485 GetPHOSGeometry()->ImpactOnEmc(vtx2,phot2->Theta(),phot2->Phi(),mod,z,x))
1486 inacceptance = kTRUE;
1487 if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
1491 if(GetFiducialCut()->IsInFiducialCut(lv1,fCalorimeter) && GetFiducialCut()->IsInFiducialCut(lv2,fCalorimeter))
1492 inacceptance = kTRUE ;
1493 if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
1496 else if(fCalorimeter == "EMCAL" && GetCaloUtils()->IsEMCALGeoMatrixSet())
1498 if(GetEMCALGeometry())
1503 GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(phot1->Eta(),phot1->Phi(),absID1);
1504 GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(phot2->Eta(),phot2->Phi(),absID2);
1506 if( absID1 >= 0 && absID2 >= 0)
1507 inacceptance = kTRUE;
1509 if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
1513 if(GetFiducialCut()->IsInFiducialCut(lv1,fCalorimeter) && GetFiducialCut()->IsInFiducialCut(lv2,fCalorimeter))
1514 inacceptance = kTRUE ;
1515 if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
1521 Float_t asym = TMath::Abs((lv1.E()-lv2.E()) / (lv1.E()+lv2.E()));
1522 Double_t angle = lv1.Angle(lv2.Vect());
1526 // printf("ACCEPTED pi0: pt %2.2f, phi %3.2f, eta %1.2f\n",pi0Pt,phi,pi0Y);
1527 fhPrimPi0AccE ->Fill(pi0E ) ;
1528 fhPrimPi0AccPt ->Fill(pi0Pt) ;
1529 fhPrimPi0AccPhi->Fill(pi0Pt, phi) ;
1530 fhPrimPi0AccY ->Fill(pi0Pt, pi0Y) ;
1531 fhPrimPi0AccPtCentrality->Fill(pi0Pt,cen) ;
1532 fhPrimPi0AccPtEventPlane->Fill(pi0Pt,ep ) ;
1536 fhPrimPi0OpeningAngle ->Fill(pi0Pt,angle);
1537 if(pi0Pt > 5)fhPrimPi0OpeningAngleAsym->Fill(asym,angle);
1538 fhPrimPi0CosOpeningAngle ->Fill(pi0Pt,TMath::Cos(angle));
1543 fhPrimEtaAccE ->Fill(pi0E ) ;
1544 fhPrimEtaAccPt ->Fill(pi0Pt) ;
1545 fhPrimEtaAccPhi->Fill(pi0Pt, phi) ;
1546 fhPrimEtaAccY ->Fill(pi0Pt, pi0Y) ;
1547 fhPrimEtaAccPtCentrality->Fill(pi0Pt,cen) ;
1548 fhPrimEtaAccPtEventPlane->Fill(pi0Pt,ep ) ;
1552 fhPrimEtaOpeningAngle ->Fill(pi0Pt,angle);
1553 if(pi0Pt > 5)fhPrimEtaOpeningAngleAsym->Fill(asym,angle);
1554 fhPrimEtaCosOpeningAngle ->Fill(pi0Pt,TMath::Cos(angle));
1559 }//Check daughters exist
1560 }// Primary pi0 or eta
1561 }//loop on primaries
1562 }//stack exists and data is MC
1568 //______________________________________________________________________________________
1569 void AliAnaPi0::FillArmenterosThetaStar(const Int_t pdg, const TLorentzVector meson,
1570 const TLorentzVector daugh1, const TLorentzVector daugh2)
1572 // Fill armenteros plots
1574 // Get pTArm and AlphaArm
1575 Float_t momentumSquaredMother = meson.P()*meson.P();
1576 Float_t momentumDaughter1AlongMother = 0.;
1577 Float_t momentumDaughter2AlongMother = 0.;
1579 if (momentumSquaredMother > 0.)
1581 momentumDaughter1AlongMother = (daugh1.Px()*meson.Px() + daugh1.Py()*meson.Py()+ daugh1.Pz()*meson.Pz()) / sqrt(momentumSquaredMother);
1582 momentumDaughter2AlongMother = (daugh2.Px()*meson.Px() + daugh2.Py()*meson.Py()+ daugh2.Pz()*meson.Pz()) / sqrt(momentumSquaredMother);
1585 Float_t momentumSquaredDaughter1 = daugh1.P()*daugh1.P();
1586 Float_t ptArmSquared = momentumSquaredDaughter1 - momentumDaughter1AlongMother*momentumDaughter1AlongMother;
1589 if (ptArmSquared > 0.)
1590 pTArm = sqrt(ptArmSquared);
1592 Float_t alphaArm = 0.;
1593 if(momentumDaughter1AlongMother +momentumDaughter2AlongMother > 0)
1594 alphaArm = (momentumDaughter1AlongMother -momentumDaughter2AlongMother) / (momentumDaughter1AlongMother + momentumDaughter2AlongMother);
1596 TLorentzVector daugh1Boost = daugh1;
1597 daugh1Boost.Boost(-meson.BoostVector());
1598 Float_t cosThStar=TMath::Cos(daugh1Boost.Vect().Angle(meson.Vect()));
1600 Float_t en = meson.Energy();
1602 if(en > 8 && en <= 12) ebin = 0;
1603 if(en > 12 && en <= 16) ebin = 1;
1604 if(en > 16 && en <= 20) ebin = 2;
1605 if(en > 20) ebin = 3;
1606 if(ebin < 0 || ebin > 3) return ;
1609 if(daugh1.E()+daugh2.E() > 0) asym = TMath::Abs(daugh1.E()-daugh2.E())/(daugh1.E()+daugh2.E());
1613 fhCosThStarPrimPi0->Fill(en,cosThStar);
1614 fhArmPrimPi0[ebin]->Fill(alphaArm,pTArm);
1618 fhCosThStarPrimEta->Fill(en,cosThStar);
1619 fhArmPrimEta[ebin]->Fill(alphaArm,pTArm);
1622 if(GetDebug() > 2 ) printf("AliAnaPi0::FillArmenterosThetaStar() - E %f, alphaArm %f, pTArm %f, cos(theta*) %f\n",en,alphaArm,pTArm,cosThStar);
1627 //_____________________________________________________________
1628 void AliAnaPi0::FillMCVersusRecDataHistograms(const Int_t index1, const Int_t index2,
1629 const Float_t pt1, const Float_t pt2,
1630 const Int_t ncell1, const Int_t ncell2,
1631 const Double_t mass, const Double_t pt, const Double_t asym,
1632 const Double_t deta, const Double_t dphi){
1633 //Do some MC checks on the origin of the pair, is there any common ancestor and if there is one, who?
1634 //Adjusted for Pythia, need to see what to do for other generators.
1635 //Array of histograms ordered as follows: 0-Photon, 1-electron, 2-pi0, 3-eta, 4-a-proton, 5-a-neutron, 6-stable particles,
1636 // 7-other decays, 8-string, 9-final parton, 10-initial parton, intermediate, 11-colliding proton, 12-unrelated
1638 if(!fFillOriginHisto) return;
1641 Int_t ancStatus = 0;
1642 TLorentzVector ancMomentum;
1643 TVector3 prodVertex;
1644 Int_t ancLabel = GetMCAnalysisUtils()->CheckCommonAncestor(index1, index2,
1645 GetReader(), ancPDG, ancStatus,ancMomentum, prodVertex);
1647 Int_t momindex = -1;
1649 Int_t momstatus = -1;
1650 if(GetDebug() > 1) printf("AliAnaPi0::FillMCVersusRecDataHistograms() - Common ancestor label %d, pdg %d, name %s, status %d; \n",
1651 ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus);
1654 if(ancPDG==22){//gamma
1655 fhMCOrgMass[0]->Fill(pt,mass);
1656 fhMCOrgAsym[0]->Fill(pt,asym);
1657 fhMCOrgDeltaEta[0]->Fill(pt,deta);
1658 fhMCOrgDeltaPhi[0]->Fill(pt,dphi);
1660 else if(TMath::Abs(ancPDG)==11){//e
1661 fhMCOrgMass[1]->Fill(pt,mass);
1662 fhMCOrgAsym[1]->Fill(pt,asym);
1663 fhMCOrgDeltaEta[1]->Fill(pt,deta);
1664 fhMCOrgDeltaPhi[1]->Fill(pt,dphi);
1666 else if(ancPDG==111){//Pi0
1667 fhMCOrgMass[2]->Fill(pt,mass);
1668 fhMCOrgAsym[2]->Fill(pt,asym);
1669 fhMCOrgDeltaEta[2]->Fill(pt,deta);
1670 fhMCOrgDeltaPhi[2]->Fill(pt,dphi);
1671 if(fMultiCutAnaSim){
1672 for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
1673 for(Int_t icell=0; icell<fNCellNCuts; icell++){
1674 for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
1675 Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
1676 if(pt1 > fPtCuts[ipt] && pt2 > fPtCuts[ipt] &&
1677 asym < fAsymCuts[iasym] &&
1678 ncell1 >= fCellNCuts[icell] && ncell2 >= fCellNCuts[icell]){
1679 fhMCPi0MassPtRec [index]->Fill(pt,mass);
1680 fhMCPi0MassPtTrue[index]->Fill(ancMomentum.Pt(),mass);
1681 if(mass < 0.17 && mass > 0.1) fhMCPi0PtTruePtRec[index]->Fill(ancMomentum.Pt(),pt);
1682 }//pass the different cuts
1683 }// pid bit cut loop
1686 }//Multi cut ana sim
1688 fhMCPi0MassPtTrue[0]->Fill(ancMomentum.Pt(),mass);
1689 if(mass < 0.17 && mass > 0.1) {
1690 fhMCPi0PtTruePtRec[0]->Fill(ancMomentum.Pt(),pt);
1691 if(fFillOriginHisto)
1693 if(GetReader()->ReadStack())
1695 TParticle* ancestor = GetMCStack()->Particle(ancLabel);
1696 momindex = ancestor->GetFirstMother();
1697 if(momindex < 0) return;
1698 TParticle* mother = GetMCStack()->Particle(momindex);
1699 mompdg = TMath::Abs(mother->GetPdgCode());
1700 momstatus = mother->GetStatusCode();
1704 TClonesArray * mcparticles = GetReader()->GetAODMCParticles();
1705 AliAODMCParticle* ancestor = (AliAODMCParticle *) mcparticles->At(ancLabel);
1706 momindex = ancestor->GetMother();
1707 if(momindex < 0) return;
1708 AliAODMCParticle* mother = (AliAODMCParticle *) mcparticles->At(momindex);
1709 mompdg = TMath::Abs(mother->GetPdgCode());
1710 momstatus = mother->GetStatus();
1713 if (momstatus == 21) fhMCPi0PtOrigin->Fill(pt,0.5);//parton
1714 else if(mompdg < 22 ) fhMCPi0PtOrigin->Fill(pt,1.5);//quark
1715 else if(mompdg > 2100 && mompdg < 2210) fhMCPi0PtOrigin->Fill(pt,2.5);// resonances
1716 else if(mompdg == 221) fhMCPi0PtOrigin->Fill(pt,8.5);//eta
1717 else if(mompdg == 331) fhMCPi0PtOrigin->Fill(pt,9.5);//eta prime
1718 else if(mompdg == 213) fhMCPi0PtOrigin->Fill(pt,4.5);//rho
1719 else if(mompdg == 223) fhMCPi0PtOrigin->Fill(pt,5.5);//omega
1720 else if(mompdg >= 310 && mompdg <= 323) fhMCPi0PtOrigin->Fill(pt,6.5);//k0S, k+-,k*
1721 else if(mompdg == 130) fhMCPi0PtOrigin->Fill(pt,6.5);//k0L
1722 else if(momstatus == 11 || momstatus == 12 ) fhMCPi0PtOrigin->Fill(pt,3.5);//resonances
1723 else fhMCPi0PtOrigin->Fill(pt,7.5);//other?
1729 else if(ancPDG==221){//Eta
1730 fhMCOrgMass[3]->Fill(pt,mass);
1731 fhMCOrgAsym[3]->Fill(pt,asym);
1732 fhMCOrgDeltaEta[3]->Fill(pt,deta);
1733 fhMCOrgDeltaPhi[3]->Fill(pt,dphi);
1734 if(fMultiCutAnaSim){
1735 for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
1736 for(Int_t icell=0; icell<fNCellNCuts; icell++){
1737 for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
1738 Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
1739 if(pt1 > fPtCuts[ipt] && pt2 > fPtCuts[ipt] &&
1740 asym < fAsymCuts[iasym] &&
1741 ncell1 >= fCellNCuts[icell] && ncell2 >= fCellNCuts[icell]){
1742 fhMCEtaMassPtRec [index]->Fill(pt,mass);
1743 fhMCEtaMassPtTrue[index]->Fill(ancMomentum.Pt(),mass);
1744 if(mass < 0.65 && mass > 0.45) fhMCEtaPtTruePtRec[index]->Fill(ancMomentum.Pt(),pt);
1745 }//pass the different cuts
1746 }// pid bit cut loop
1749 } //Multi cut ana sim
1751 fhMCEtaMassPtTrue[0]->Fill(ancMomentum.Pt(),mass);
1752 if(mass < 0.65 && mass > 0.45) fhMCEtaPtTruePtRec[0]->Fill(ancMomentum.Pt(),pt);
1754 if(fFillOriginHisto)
1756 if(GetReader()->ReadStack())
1758 TParticle* ancestor = GetMCStack()->Particle(ancLabel);
1759 momindex = ancestor->GetFirstMother();
1760 if(momindex < 0) return;
1761 TParticle* mother = GetMCStack()->Particle(momindex);
1762 mompdg = TMath::Abs(mother->GetPdgCode());
1763 momstatus = mother->GetStatusCode();
1767 TClonesArray * mcparticles = GetReader()->GetAODMCParticles();
1768 AliAODMCParticle* ancestor = (AliAODMCParticle *) mcparticles->At(ancLabel);
1769 momindex = ancestor->GetMother();
1770 if(momindex < 0) return;
1771 AliAODMCParticle* mother = (AliAODMCParticle *) mcparticles->At(momindex);
1772 mompdg = TMath::Abs(mother->GetPdgCode());
1773 momstatus = mother->GetStatus();
1776 if (momstatus == 21 ) fhMCEtaPtOrigin->Fill(pt,0.5);//parton
1777 else if(mompdg < 22 ) fhMCEtaPtOrigin->Fill(pt,1.5);//quark
1778 else if(mompdg > 2100 && mompdg < 2210) fhMCEtaPtOrigin->Fill(pt,2.5);//qq resonances
1779 else if(mompdg == 331) fhMCEtaPtOrigin->Fill(pt,5.5);//eta prime
1780 else if(momstatus == 11 || momstatus == 12 ) fhMCEtaPtOrigin->Fill(pt,3.5);//resonances
1781 else fhMCEtaPtOrigin->Fill(pt,4.5);//stable, conversions?
1782 //printf("Other Meson pdg %d, Mother %s, pdg %d, status %d\n",pdg, TDatabasePDG::Instance()->GetParticle(mompdg)->GetName(),mompdg, momstatus );
1786 else if(ancPDG==-2212){//AProton
1787 fhMCOrgMass[4]->Fill(pt,mass);
1788 fhMCOrgAsym[4]->Fill(pt,asym);
1789 fhMCOrgDeltaEta[4]->Fill(pt,deta);
1790 fhMCOrgDeltaPhi[4]->Fill(pt,dphi);
1792 else if(ancPDG==-2112){//ANeutron
1793 fhMCOrgMass[5]->Fill(pt,mass);
1794 fhMCOrgAsym[5]->Fill(pt,asym);
1795 fhMCOrgDeltaEta[5]->Fill(pt,deta);
1796 fhMCOrgDeltaPhi[5]->Fill(pt,dphi);
1798 else if(TMath::Abs(ancPDG)==13){//muons
1799 fhMCOrgMass[6]->Fill(pt,mass);
1800 fhMCOrgAsym[6]->Fill(pt,asym);
1801 fhMCOrgDeltaEta[6]->Fill(pt,deta);
1802 fhMCOrgDeltaPhi[6]->Fill(pt,dphi);
1804 else if (TMath::Abs(ancPDG) > 100 && ancLabel > 7) {
1805 if(ancStatus==1){//Stable particles, converted? not decayed resonances
1806 fhMCOrgMass[6]->Fill(pt,mass);
1807 fhMCOrgAsym[6]->Fill(pt,asym);
1808 fhMCOrgDeltaEta[6]->Fill(pt,deta);
1809 fhMCOrgDeltaPhi[6]->Fill(pt,dphi);
1811 else{//resonances and other decays, more hadron conversions?
1812 fhMCOrgMass[7]->Fill(pt,mass);
1813 fhMCOrgAsym[7]->Fill(pt,asym);
1814 fhMCOrgDeltaEta[7]->Fill(pt,deta);
1815 fhMCOrgDeltaPhi[7]->Fill(pt,dphi);
1818 else {//Partons, colliding protons, strings, intermediate corrections
1819 if(ancStatus==11 || ancStatus==12){//String fragmentation
1820 fhMCOrgMass[8]->Fill(pt,mass);
1821 fhMCOrgAsym[8]->Fill(pt,asym);
1822 fhMCOrgDeltaEta[8]->Fill(pt,deta);
1823 fhMCOrgDeltaPhi[8]->Fill(pt,dphi);
1825 else if (ancStatus==21){
1826 if(ancLabel < 2) {//Colliding protons
1827 fhMCOrgMass[11]->Fill(pt,mass);
1828 fhMCOrgAsym[11]->Fill(pt,asym);
1829 fhMCOrgDeltaEta[11]->Fill(pt,deta);
1830 fhMCOrgDeltaPhi[11]->Fill(pt,dphi);
1831 }//colliding protons
1832 else if(ancLabel < 6){//partonic initial states interactions
1833 fhMCOrgMass[9]->Fill(pt,mass);
1834 fhMCOrgAsym[9]->Fill(pt,asym);
1835 fhMCOrgDeltaEta[9]->Fill(pt,deta);
1836 fhMCOrgDeltaPhi[9]->Fill(pt,dphi);
1838 else if(ancLabel < 8){//Final state partons radiations?
1839 fhMCOrgMass[10]->Fill(pt,mass);
1840 fhMCOrgAsym[10]->Fill(pt,asym);
1841 fhMCOrgDeltaEta[10]->Fill(pt,deta);
1842 fhMCOrgDeltaPhi[10]->Fill(pt,dphi);
1845 // printf("AliAnaPi0::FillMCVersusRecDataHistograms() - Check ** Common ancestor label %d, pdg %d, name %s, status %d; \n",
1846 // ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus);
1850 // printf("AliAnaPi0::FillMCVersusRecDataHistograms() - Check *** Common ancestor label %d, pdg %d, name %s, status %d; \n",
1851 // ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus);
1853 }////Partons, colliding protons, strings, intermediate corrections
1855 else { //ancLabel <= -1
1856 //printf("Not related at all label = %d\n",ancLabel);
1857 fhMCOrgMass[12]->Fill(pt,mass);
1858 fhMCOrgAsym[12]->Fill(pt,asym);
1859 fhMCOrgDeltaEta[12]->Fill(pt,deta);
1860 fhMCOrgDeltaPhi[12]->Fill(pt,dphi);
1864 //__________________________________________
1865 void AliAnaPi0::MakeAnalysisFillHistograms()
1867 //Process one event and extract photons from AOD branch
1868 // filled with AliAnaPhoton and fill histos with invariant mass
1870 //In case of simulated data, fill acceptance histograms
1871 if(IsDataMC())FillAcceptanceHistograms();
1873 //if (GetReader()->GetEventNumber()%10000 == 0)
1874 // printf("--- Event %d ---\n",GetReader()->GetEventNumber());
1876 if(!GetInputAODBranch())
1878 printf("AliAnaPi0::MakeAnalysisFillHistograms() - No input aod photons in AOD with name branch < %s >, STOP \n",GetInputAODName().Data());
1882 //Init some variables
1883 Int_t nPhot = GetInputAODBranch()->GetEntriesFast() ;
1886 printf("AliAnaPi0::MakeAnalysisFillHistograms() - Photon entries %d\n", nPhot);
1888 //If less than photon 2 entries in the list, skip this event
1892 printf("AliAnaPi0::MakeAnalysisFillHistograms() - nPhotons %d, cent bin %d continue to next event\n",nPhot, GetEventCentrality());
1894 if(GetNCentrBin() > 1) fhCentralityNoPair->Fill(GetEventCentrality() * GetNCentrBin() / GetReader()->GetCentralityOpt());
1899 Int_t ncentr = GetNCentrBin();
1900 if(GetNCentrBin()==0) ncentr = 1;
1905 Double_t vert[] = {0.0, 0.0, 0.0} ; //vertex
1906 Int_t evtIndex1 = 0 ;
1907 Int_t currentEvtIndex = -1;
1908 Int_t curCentrBin = GetEventCentralityBin();
1909 //Int_t curVzBin = GetEventVzBin();
1910 //Int_t curRPBin = GetEventRPBin();
1911 Int_t eventbin = GetEventMixBin();
1913 //Get shower shape information of clusters
1914 TObjArray *clusters = 0;
1915 if (fCalorimeter=="EMCAL") clusters = GetEMCALClusters();
1916 else if(fCalorimeter=="PHOS" ) clusters = GetPHOSClusters() ;
1918 //---------------------------------
1919 //First loop on photons/clusters
1920 //---------------------------------
1921 for(Int_t i1=0; i1<nPhot-1; i1++)
1923 AliAODPWG4Particle * p1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i1)) ;
1924 //printf("AliAnaPi0::MakeAnalysisFillHistograms() : cluster1 id %d\n",p1->GetCaloLabel(0));
1926 // get the event index in the mixed buffer where the photon comes from
1927 // in case of mixing with analysis frame, not own mixing
1928 evtIndex1 = GetEventIndex(p1, vert) ;
1929 //printf("charge = %d\n", track->Charge());
1930 if ( evtIndex1 == -1 )
1932 if ( evtIndex1 == -2 )
1935 //printf("z vertex %f < %f\n",vert[2],GetZvertexCut());
1936 if(TMath::Abs(vert[2]) > GetZvertexCut()) continue ; //vertex cut
1939 if (evtIndex1 != currentEvtIndex)
1941 //Fill event bin info
1942 fhEventBin->Fill(eventbin) ;
1943 if(GetNCentrBin() > 1)
1945 fhCentrality->Fill(curCentrBin);
1946 if(GetNRPBin() > 1 && GetEventPlane()) fhEventPlaneResolution->Fill(curCentrBin,TMath::Cos(2.*GetEventPlane()->GetQsubRes()));
1948 currentEvtIndex = evtIndex1 ;
1951 //printf("AliAnaPi0::MakeAnalysisFillHistograms(): Photon 1 Evt %d Vertex : %f,%f,%f\n",evtIndex1, GetVertex(evtIndex1)[0] ,GetVertex(evtIndex1)[1],GetVertex(evtIndex1)[2]);
1953 //Get the momentum of this cluster
1954 TLorentzVector photon1(p1->Px(),p1->Py(),p1->Pz(),p1->E());
1956 //Get (Super)Module number of this cluster
1957 module1 = GetModuleNumber(p1);
1959 //------------------------------------------
1960 //Get index in VCaloCluster array
1961 AliVCluster *cluster1 = 0;
1962 Bool_t bFound1 = kFALSE;
1963 Int_t caloLabel1 = p1->GetCaloLabel(0);
1967 for(Int_t iclus = 0; iclus < clusters->GetEntriesFast(); iclus++){
1968 AliVCluster *cluster= dynamic_cast<AliVCluster*> (clusters->At(iclus));
1971 if (cluster->GetID()==caloLabel1)
1980 }// calorimeter clusters loop
1982 //---------------------------------
1983 //Second loop on photons/clusters
1984 //---------------------------------
1985 for(Int_t i2=i1+1; i2<nPhot; i2++)
1987 AliAODPWG4Particle * p2 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i2)) ;
1989 //In case of mixing frame, check we are not in the same event as the first cluster
1990 Int_t evtIndex2 = GetEventIndex(p2, vert) ;
1991 if ( evtIndex2 == -1 )
1993 if ( evtIndex2 == -2 )
1995 if (GetMixedEvent() && (evtIndex1 == evtIndex2))
1998 //------------------------------------------
1999 //Get index in VCaloCluster array
2000 AliVCluster *cluster2 = 0;
2001 Bool_t bFound2 = kFALSE;
2002 Int_t caloLabel2 = p2->GetCaloLabel(0);
2004 for(Int_t iclus = iclus1+1; iclus < clusters->GetEntriesFast(); iclus++){
2005 AliVCluster *cluster= dynamic_cast<AliVCluster*> (clusters->At(iclus));
2007 if(cluster->GetID()==caloLabel2) {
2013 }// calorimeter clusters loop
2018 if(cluster1 && bFound1){
2019 tof1 = cluster1->GetTOF()*1e9;
2020 l01 = cluster1->GetM02();
2022 // else printf("cluster1 not available: calo label %d / %d, cluster ID %d\n",
2023 // p1->GetCaloLabel(0),(GetReader()->GetInputEvent())->GetNumberOfCaloClusters()-1,cluster1->GetID());
2027 if(cluster2 && bFound2)
2029 tof2 = cluster2->GetTOF()*1e9;
2030 l02 = cluster2->GetM02();
2033 // else printf("cluster2 not available: calo label %d / %d, cluster ID %d\n",
2034 // p2->GetCaloLabel(0),(GetReader()->GetInputEvent())->GetNumberOfCaloClusters()-1,cluster2->GetID());
2038 Double_t t12diff = tof1-tof2;
2039 if(TMath::Abs(t12diff) > GetPairTimeCut()) continue;
2041 //------------------------------------------
2043 //printf("AliAnaPi0::MakeAnalysisFillHistograms(): Photon 2 Evt %d Vertex : %f,%f,%f\n",evtIndex2, GetVertex(evtIndex2)[0] ,GetVertex(evtIndex2)[1],GetVertex(evtIndex2)[2]);
2045 //Get the momentum of this cluster
2046 TLorentzVector photon2(p2->Px(),p2->Py(),p2->Pz(),p2->E());
2048 module2 = GetModuleNumber(p2);
2050 //---------------------------------
2051 // Get pair kinematics
2052 //---------------------------------
2053 Double_t m = (photon1 + photon2).M() ;
2054 Double_t pt = (photon1 + photon2).Pt();
2055 Double_t deta = photon1.Eta() - photon2.Eta();
2056 Double_t dphi = photon1.Phi() - photon2.Phi();
2057 Double_t a = TMath::Abs(p1->E()-p2->E())/(p1->E()+p2->E()) ;
2060 printf(" E: photon1 %f, photon2 %f; Pair: pT %f, mass %f, a %f\n", p1->E(), p2->E(), (photon1 + photon2).E(),m,a);
2062 //--------------------------------
2063 // Opening angle selection
2064 //--------------------------------
2065 //Check if opening angle is too large or too small compared to what is expected
2066 Double_t angle = photon1.Angle(photon2.Vect());
2067 if(fUseAngleEDepCut && !GetNeutralMesonSelection()->IsAngleInWindow((photon1+photon2).E(),angle+0.05)) {
2069 printf("AliAnaPi0::MakeAnalysisFillHistograms() -Real pair angle %f not in E %f window\n",angle, (photon1+photon2).E());
2073 if(fUseAngleCut && (angle < fAngleCut || angle > fAngleMaxCut)) {
2075 printf("AliAnaPi0::MakeAnalysisFillHistograms() - Real pair cut %f < angle %f < cut %f\n",fAngleCut, angle, fAngleMaxCut);
2079 //-------------------------------------------------------------------------------------------------
2080 //Fill module dependent histograms, put a cut on assymmetry on the first available cut in the array
2081 //-------------------------------------------------------------------------------------------------
2082 if(a < fAsymCuts[0] && fFillSMCombinations)
2084 if(module1==module2 && module1 >=0 && module1<fNModules)
2085 fhReMod[module1]->Fill(pt,m) ;
2087 if(fCalorimeter=="EMCAL")
2091 for(Int_t i = 0; i < fNModules/2; i++)
2094 if((module1==j && module2==j+1) || (module1==j+1 && module2==j)) fhReSameSectorEMCALMod[i]->Fill(pt,m) ;
2098 for(Int_t i = 0; i < fNModules-2; i++){
2099 if((module1==i && module2==i+2) || (module1==i+2 && module2==i)) fhReSameSideEMCALMod[i]->Fill(pt,m);
2103 if((module1==0 && module2==1) || (module1==1 && module2==0)) fhReDiffPHOSMod[0]->Fill(pt,m) ;
2104 if((module1==0 && module2==2) || (module1==2 && module2==0)) fhReDiffPHOSMod[1]->Fill(pt,m) ;
2105 if((module1==1 && module2==2) || (module1==2 && module2==1)) fhReDiffPHOSMod[2]->Fill(pt,m) ;
2109 //In case we want only pairs in same (super) module, check their origin.
2111 if(fSameSM && module1!=module2) ok=kFALSE;
2114 //Check if one of the clusters comes from a conversion
2115 if(fCheckConversion)
2117 if (p1->IsTagged() && p2->IsTagged()) fhReConv2->Fill(pt,m);
2118 else if(p1->IsTagged() || p2->IsTagged()) fhReConv ->Fill(pt,m);
2121 // Fill shower shape cut histograms
2122 if(fFillSSCombinations)
2124 if ( l01 > 0.01 && l01 < 0.4 &&
2125 l02 > 0.01 && l02 < 0.4 ) fhReSS[0]->Fill(pt,m); // Tight
2126 else if( l01 > 0.4 && l02 > 0.4 ) fhReSS[1]->Fill(pt,m); // Loose
2127 else if( l01 > 0.01 && l01 < 0.4 && l02 > 0.4 ) fhReSS[2]->Fill(pt,m); // Both
2128 else if( l02 > 0.01 && l02 < 0.4 && l01 > 0.4 ) fhReSS[2]->Fill(pt,m); // Both
2131 //Fill histograms for different bad channel distance, centrality, assymmetry cut and pid bit
2132 for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
2133 if((p1->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton)) && (p2->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton))){
2134 for(Int_t iasym=0; iasym < fNAsymCuts; iasym++){
2135 if(a < fAsymCuts[iasym])
2137 Int_t index = ((curCentrBin*fNPIDBits)+ipid)*fNAsymCuts + iasym;
2138 //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);
2140 if(index < 0 || index >= ncentr*fNPIDBits*fNAsymCuts) continue ;
2142 fhRe1 [index]->Fill(pt,m);
2143 if(fMakeInvPtPlots)fhReInvPt1[index]->Fill(pt,m,1./pt) ;
2144 if(fFillBadDistHisto){
2145 if(p1->DistToBad()>0 && p2->DistToBad()>0){
2146 fhRe2 [index]->Fill(pt,m) ;
2147 if(fMakeInvPtPlots)fhReInvPt2[index]->Fill(pt,m,1./pt) ;
2148 if(p1->DistToBad()>1 && p2->DistToBad()>1){
2149 fhRe3 [index]->Fill(pt,m) ;
2150 if(fMakeInvPtPlots)fhReInvPt3[index]->Fill(pt,m,1./pt) ;
2153 }// Fill bad dist histos
2155 }// asymmetry cut loop
2159 //Fill histograms with opening angle
2162 fhRealOpeningAngle ->Fill(pt,angle);
2163 fhRealCosOpeningAngle->Fill(pt,TMath::Cos(angle));
2166 //Fill histograms with pair assymmetry
2167 if(fFillAsymmetryHisto)
2169 fhRePtAsym->Fill(pt,a);
2170 if(m > 0.10 && m < 0.17) fhRePtAsymPi0->Fill(pt,a);
2171 if(m > 0.45 && m < 0.65) fhRePtAsymEta->Fill(pt,a);
2174 //-------------------------------------------------------
2175 //Get the number of cells needed for multi cut analysis.
2176 //-------------------------------------------------------
2179 if(fMultiCutAna || (IsDataMC() && fMultiCutAnaSim))
2181 AliVEvent * event = GetReader()->GetInputEvent();
2183 for(Int_t iclus = 0; iclus < event->GetNumberOfCaloClusters(); iclus++)
2185 AliVCluster *cluster = event->GetCaloCluster(iclus);
2188 if (fCalorimeter == "EMCAL" && GetReader()->IsEMCALCluster(cluster)) is = kTRUE;
2189 else if(fCalorimeter == "PHOS" && GetReader()->IsPHOSCluster (cluster)) is = kTRUE;
2192 if (p1->GetCaloLabel(0) == cluster->GetID()) ncell1 = cluster->GetNCells();
2193 else if (p2->GetCaloLabel(0) == cluster->GetID()) ncell2 = cluster->GetNCells();
2194 } // PHOS or EMCAL cluster as requested in analysis
2196 if(ncell2 > 0 && ncell1 > 0) break; // No need to continue the iteration
2199 //printf("e 1: %2.2f, e 2: %2.2f, ncells: n1 %d, n2 %d\n", p1->E(), p2->E(),ncell1,ncell2);
2206 //Do some MC checks on the origin of the pair, is there any common ancestor and if there is one, who?
2209 if(GetMCAnalysisUtils()->CheckTagBit(p1->GetTag(),AliMCAnalysisUtils::kMCConversion) &&
2210 GetMCAnalysisUtils()->CheckTagBit(p2->GetTag(),AliMCAnalysisUtils::kMCConversion))
2212 fhReMCFromConversion->Fill(pt,m);
2214 else if(!GetMCAnalysisUtils()->CheckTagBit(p1->GetTag(),AliMCAnalysisUtils::kMCConversion) &&
2215 !GetMCAnalysisUtils()->CheckTagBit(p2->GetTag(),AliMCAnalysisUtils::kMCConversion))
2217 fhReMCFromNotConversion->Fill(pt,m);
2221 fhReMCFromMixConversion->Fill(pt,m);
2224 FillMCVersusRecDataHistograms(p1->GetLabel(), p2->GetLabel(),p1->Pt(), p2->Pt(),ncell1, ncell2, m, pt, a,deta, dphi);
2227 //-----------------------
2228 //Multi cuts analysis
2229 //-----------------------
2232 //Histograms for different PID bits selection
2233 for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
2235 if(p1->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton) &&
2236 p2->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton)) fhRePIDBits[ipid]->Fill(pt,m) ;
2238 //printf("ipt %d, ipid%d, name %s\n",ipt, ipid, fhRePtPIDCuts[ipt*fNPIDBitsBits+ipid]->GetName());
2239 } // pid bit cut loop
2241 //Several pt,ncell and asymmetry cuts
2242 for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
2243 for(Int_t icell=0; icell<fNCellNCuts; icell++){
2244 for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
2245 Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
2246 if(p1->E() > fPtCuts[ipt] && p2->E() > fPtCuts[ipt] &&
2247 a < fAsymCuts[iasym] &&
2248 ncell1 >= fCellNCuts[icell] && ncell2 >= fCellNCuts[icell]){
2249 fhRePtNCellAsymCuts[index]->Fill(pt,m) ;
2250 //printf("ipt %d, icell%d, iasym %d, name %s\n",ipt, icell, iasym, fhRePtNCellAsymCuts[((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym]->GetName());
2251 if(fFillSMCombinations && module1==module2){
2252 fhRePtNCellAsymCutsSM[module1][index]->Fill(pt,m) ;
2255 }// pid bit cut loop
2258 if(GetHistogramRanges()->GetHistoTrackMultiplicityBins()){
2259 for(Int_t iasym = 0; iasym < fNAsymCuts; iasym++){
2260 if(a < fAsymCuts[iasym])fhRePtMult[iasym]->Fill(pt,GetTrackMultiplicity(),m) ;
2263 }// multiple cuts analysis
2265 }// second same event particle
2268 //-------------------------------------------------------------
2270 //-------------------------------------------------------------
2273 //Recover events in with same characteristics as the current event
2275 //Check that the bin exists, if not (bad determination of RP, centrality or vz bin) do nothing
2276 if(eventbin < 0) return ;
2278 TList * evMixList=fEventsList[eventbin] ;
2282 printf("AliAnaPi0::MakeAnalysisFillHistograms() - Mix event list not available, bin %d \n",eventbin);
2286 Int_t nMixed = evMixList->GetSize() ;
2287 for(Int_t ii=0; ii<nMixed; ii++)
2289 TClonesArray* ev2= (TClonesArray*) (evMixList->At(ii));
2290 Int_t nPhot2=ev2->GetEntriesFast() ;
2293 printf("AliAnaPi0::MakeAnalysisFillHistograms() - Mixed event %d photon entries %d, centrality bin %d\n", ii, nPhot2, GetEventCentralityBin());
2295 fhEventMixBin->Fill(eventbin) ;
2297 //---------------------------------
2298 //First loop on photons/clusters
2299 //---------------------------------
2300 for(Int_t i1=0; i1<nPhot; i1++){
2301 AliAODPWG4Particle * p1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i1)) ;
2302 if(fSameSM && GetModuleNumber(p1)!=module1) continue;
2304 //Get kinematics of cluster and (super) module of this cluster
2305 TLorentzVector photon1(p1->Px(),p1->Py(),p1->Pz(),p1->E());
2306 module1 = GetModuleNumber(p1);
2308 //---------------------------------
2309 //First loop on photons/clusters
2310 //---------------------------------
2311 for(Int_t i2=0; i2<nPhot2; i2++){
2312 AliAODPWG4Particle * p2 = (AliAODPWG4Particle*) (ev2->At(i2)) ;
2314 //Get kinematics of second cluster and calculate those of the pair
2315 TLorentzVector photon2(p2->Px(),p2->Py(),p2->Pz(),p2->E());
2316 m = (photon1+photon2).M() ;
2317 Double_t pt = (photon1 + photon2).Pt();
2318 Double_t a = TMath::Abs(p1->E()-p2->E())/(p1->E()+p2->E()) ;
2320 //Check if opening angle is too large or too small compared to what is expected
2321 Double_t angle = photon1.Angle(photon2.Vect());
2322 if(fUseAngleEDepCut && !GetNeutralMesonSelection()->IsAngleInWindow((photon1+photon2).E(),angle+0.05)){
2324 printf("AliAnaPi0::MakeAnalysisFillHistograms() -Mix pair angle %f not in E %f window\n",angle, (photon1+photon2).E());
2327 if(fUseAngleCut && (angle < fAngleCut || angle > fAngleMaxCut)) {
2329 printf("AliAnaPi0::MakeAnalysisFillHistograms() -Mix pair angle %f < cut %f\n",angle,fAngleCut);
2335 printf("AliAnaPi0::MakeAnalysisFillHistograms() - Mixed Event: pT: photon1 %2.2f, photon2 %2.2f; Pair: pT %2.2f, mass %2.3f, a %f2.3\n",
2336 p1->Pt(), p2->Pt(), pt,m,a);
2338 //In case we want only pairs in same (super) module, check their origin.
2339 module2 = GetModuleNumber(p2);
2341 //-------------------------------------------------------------------------------------------------
2342 //Fill module dependent histograms, put a cut on assymmetry on the first available cut in the array
2343 //-------------------------------------------------------------------------------------------------
2344 if(a < fAsymCuts[0] && fFillSMCombinations){
2345 if(module1==module2 && module1 >=0 && module1<fNModules)
2346 fhMiMod[module1]->Fill(pt,m) ;
2348 if(fCalorimeter=="EMCAL"){
2352 for(Int_t i = 0; i < fNModules/2; i++){
2354 if((module1==j && module2==j+1) || (module1==j+1 && module2==j)) fhMiSameSectorEMCALMod[i]->Fill(pt,m) ;
2358 for(Int_t i = 0; i < fNModules-2; i++){
2359 if((module1==i && module2==i+2) || (module1==i+2 && module2==i)) fhMiSameSideEMCALMod[i]->Fill(pt,m);
2363 if((module1==0 && module2==1) || (module1==1 && module2==0)) fhMiDiffPHOSMod[0]->Fill(pt,m) ;
2364 if((module1==0 && module2==2) || (module1==2 && module2==0)) fhMiDiffPHOSMod[1]->Fill(pt,m) ;
2365 if((module1==1 && module2==2) || (module1==2 && module2==1)) fhMiDiffPHOSMod[2]->Fill(pt,m) ;
2372 if(fSameSM && module1!=module2) ok=kFALSE;
2375 //Check if one of the clusters comes from a conversion
2376 if(fCheckConversion){
2377 if (p1->IsTagged() && p2->IsTagged()) fhMiConv2->Fill(pt,m);
2378 else if(p1->IsTagged() || p2->IsTagged()) fhMiConv ->Fill(pt,m);
2380 //Fill histograms for different bad channel distance, centrality, assymmetry cut and pid bit
2381 for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
2382 if((p1->IsPIDOK(ipid,AliCaloPID::kPhoton)) && (p2->IsPIDOK(ipid,AliCaloPID::kPhoton)))
2384 for(Int_t iasym=0; iasym < fNAsymCuts; iasym++)
2386 if(a < fAsymCuts[iasym])
2388 Int_t index = ((curCentrBin*fNPIDBits)+ipid)*fNAsymCuts + iasym;
2390 if(index < 0 || index >= ncentr*fNPIDBits*fNAsymCuts) continue ;
2392 fhMi1 [index]->Fill(pt,m) ;
2393 if(fMakeInvPtPlots)fhMiInvPt1[index]->Fill(pt,m,1./pt) ;
2394 if(fFillBadDistHisto)
2396 if(p1->DistToBad()>0 && p2->DistToBad()>0)
2398 fhMi2 [index]->Fill(pt,m) ;
2399 if(fMakeInvPtPlots)fhMiInvPt2[index]->Fill(pt,m,1./pt) ;
2400 if(p1->DistToBad()>1 && p2->DistToBad()>1)
2402 fhMi3 [index]->Fill(pt,m) ;
2403 if(fMakeInvPtPlots)fhMiInvPt3[index]->Fill(pt,m,1./pt) ;
2406 }// Fill bad dist histo
2410 }//loop for histograms
2412 //-----------------------
2413 //Multi cuts analysis
2414 //-----------------------
2416 //Several pt,ncell and asymmetry cuts
2418 for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
2419 for(Int_t icell=0; icell<fNCellNCuts; icell++){
2420 for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
2421 Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
2422 if(p1->Pt() > fPtCuts[ipt] && p2->Pt() > fPtCuts[ipt] &&
2423 a < fAsymCuts[iasym] //&&
2424 //p1->GetBtag() >= fCellNCuts[icell] && p2->GetBtag() >= fCellNCuts[icell] // trick, correct it.
2426 fhMiPtNCellAsymCuts[index]->Fill(pt,m) ;
2427 //printf("ipt %d, icell%d, iasym %d, name %s\n",ipt, icell, iasym, fhRePtNCellAsymCuts[((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym]->GetName());
2429 }// pid bit cut loop
2434 //Fill histograms with opening angle
2435 if(fFillAngleHisto){
2436 fhMixedOpeningAngle ->Fill(pt,angle);
2437 fhMixedCosOpeningAngle->Fill(pt,TMath::Cos(angle));
2441 }// second cluster loop
2442 }//first cluster loop
2443 }//loop on mixed events
2445 //--------------------------------------------------------
2446 //Add the current event to the list of events for mixing
2447 //--------------------------------------------------------
2448 TClonesArray *currentEvent = new TClonesArray(*GetInputAODBranch());
2449 //Add current event to buffer and Remove redundant events
2450 if(currentEvent->GetEntriesFast()>0){
2451 evMixList->AddFirst(currentEvent) ;
2452 currentEvent=0 ; //Now list of particles belongs to buffer and it will be deleted with buffer
2453 if(evMixList->GetSize() >= GetNMaxEvMix())
2455 TClonesArray * tmp = (TClonesArray*) (evMixList->Last()) ;
2456 evMixList->RemoveLast() ;
2461 delete currentEvent ;
2468 //________________________________________________________________________
2469 Int_t AliAnaPi0::GetEventIndex(AliAODPWG4Particle * part, Double_t * vert)
2471 // retieves the event index and checks the vertex
2472 // in the mixed buffer returns -2 if vertex NOK
2473 // for normal events returns 0 if vertex OK and -1 if vertex NOK
2475 Int_t evtIndex = -1 ;
2476 if(GetReader()->GetDataType()!=AliCaloTrackReader::kMC){
2478 if (GetMixedEvent()){
2480 evtIndex = GetMixedEvent()->EventIndexForCaloCluster(part->GetCaloLabel(0)) ;
2481 GetVertex(vert,evtIndex);
2483 if(TMath::Abs(vert[2])> GetZvertexCut())
2484 evtIndex = -2 ; //Event can not be used (vertex, centrality,... cuts not fulfilled)
2485 } else {// Single event
2489 if(TMath::Abs(vert[2])> GetZvertexCut())
2490 evtIndex = -1 ; //Event can not be used (vertex, centrality,... cuts not fulfilled)