AliAnaPhoton: add histograms for conversion photons identification, MC and real data...
[u/mrichter/AliRoot.git] / PWG4 / PartCorrDep / AliAnaPi0.cxx
CommitLineData
1c5acb87 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15/* $Id: $ */
16
17//_________________________________________________________________________
18// Class to collect two-photon invariant mass distributions for
6175da48 19// extracting raw pi0 yield.
20// Input is produced by AliAnaPhoton (or any other analysis producing output AliAODPWG4Particles),
21// it will do nothing if executed alone
1c5acb87 22//
23//-- Author: Dmitri Peressounko (RRC "KI")
24//-- Adapted to PartCorr frame by Lamia Benhabib (SUBATECH)
25//-- and Gustavo Conesa (INFN-Frascati)
26//_________________________________________________________________________
27
28
29// --- ROOT system ---
30#include "TH3.h"
50f39b97 31#include "TH2D.h"
1c5acb87 32//#include "Riostream.h"
6639984f 33#include "TCanvas.h"
34#include "TPad.h"
35#include "TROOT.h"
477d6cee 36#include "TClonesArray.h"
0c1383b5 37#include "TObjString.h"
6175da48 38#include "TDatabasePDG.h"
1c5acb87 39
40//---- AliRoot system ----
41#include "AliAnaPi0.h"
42#include "AliCaloTrackReader.h"
43#include "AliCaloPID.h"
6639984f 44#include "AliStack.h"
ff45398a 45#include "AliFiducialCut.h"
477d6cee 46#include "TParticle.h"
477d6cee 47#include "AliVEvent.h"
6921fa00 48#include "AliESDCaloCluster.h"
49#include "AliESDEvent.h"
50#include "AliAODEvent.h"
50f39b97 51#include "AliNeutralMesonSelection.h"
c8fe2783 52#include "AliMixedEvent.h"
6175da48 53#include "AliAODMCParticle.h"
591cc579 54
1c5acb87 55ClassImp(AliAnaPi0)
56
57//________________________________________________________________________________________________________________________________________________
58AliAnaPi0::AliAnaPi0() : AliAnaPartCorrBaseClass(),
5025c139 59fDoOwnMix(kFALSE),fNCentrBin(0),//fNZvertBin(0),fNrpBin(0),
af7b3903 60fNmaxMixEv(0), fCalorimeter(""),
6175da48 61fNModules(12), fUseAngleCut(kFALSE), fUseAngleEDepCut(kFALSE),fAngleCut(0), fAngleMaxCut(7.),fEventsList(0x0), fMultiCutAna(kFALSE), fMultiCutAnaSim(kFALSE),
62fNPtCuts(0),fNAsymCuts(0), fNCellNCuts(0),fNPIDBits(0), fMakeInvPtPlots(kFALSE), fSameSM(kFALSE),
63fUseTrackMultBins(kFALSE),fUsePhotonMultBins(kFALSE),fUseAverClusterEBins(kFALSE),fUseAverCellEBins(kFALSE), fFillBadDistHisto(kFALSE),
64fhAverTotECluster(0),fhAverTotECell(0),
65fhReMod(0x0), fhReDiffMod(0x0), fhMiMod(0x0), fhMiDiffMod(0x0),
66fhReConv(0x0), fhMiConv(0x0), fhReConv2(0x0), fhMiConv2(0x0),
5ae09196 67fhRe1(0x0), fhMi1(0x0), fhRe2(0x0), fhMi2(0x0), fhRe3(0x0), fhMi3(0x0),
68fhReInvPt1(0x0), fhMiInvPt1(0x0), fhReInvPt2(0x0), fhMiInvPt2(0x0), fhReInvPt3(0x0), fhMiInvPt3(0x0),
6175da48 69fhRePtNCellAsymCuts(0x0), fhRePtNCellAsymCutsSM0(0x0), fhRePtNCellAsymCutsSM1(0x0), fhRePtNCellAsymCutsSM2(0x0), fhRePtNCellAsymCutsSM3(0x0), fhMiPtNCellAsymCuts(0x0),
70fhRePIDBits(0x0),fhRePtMult(0x0), fhRePtAsym(0x0), fhRePtAsymPi0(0x0),fhRePtAsymEta(0x0),
71fhEvents(0x0), fhRealOpeningAngle(0x0),fhRealCosOpeningAngle(0x0), fhMixedOpeningAngle(0x0),fhMixedCosOpeningAngle(0x0),
72fhPrimPi0Pt(0x0), fhPrimPi0AccPt(0x0), fhPrimPi0Y(0x0), fhPrimPi0AccY(0x0), fhPrimPi0Phi(0x0), fhPrimPi0AccPhi(0x0),
73fhPrimPi0OpeningAngle(0x0),fhPrimPi0CosOpeningAngle(0x0),
74fhPrimEtaPt(0x0), fhPrimEtaAccPt(0x0), fhPrimEtaY(0x0), fhPrimEtaAccY(0x0), fhPrimEtaPhi(0x0), fhPrimEtaAccPhi(0x0),
75fhMCOrgMass(),fhMCOrgAsym(), fhMCOrgDeltaEta(),fhMCOrgDeltaPhi(),
76fhMCPi0MassPtRec(), fhMCPi0MassPtTrue(), fhMCPi0PtTruePtRec(), fhMCEtaMassPtRec(), fhMCEtaMassPtTrue(), fhMCEtaPtTruePtRec()
1c5acb87 77{
78//Default Ctor
79 InitParameters();
80
81}
1c5acb87 82
83//________________________________________________________________________________________________________________________________________________
84AliAnaPi0::~AliAnaPi0() {
477d6cee 85 // Remove event containers
7e7694bb 86
87 if(fDoOwnMix && fEventsList){
477d6cee 88 for(Int_t ic=0; ic<fNCentrBin; ic++){
5025c139 89 for(Int_t iz=0; iz<GetNZvertBin(); iz++){
90 for(Int_t irp=0; irp<GetNRPBin(); irp++){
91 fEventsList[ic*GetNZvertBin()*GetNRPBin()+iz*GetNRPBin()+irp]->Delete() ;
92 delete fEventsList[ic*GetNZvertBin()*GetNRPBin()+iz*GetNRPBin()+irp] ;
7e7694bb 93 }
477d6cee 94 }
95 }
96 delete[] fEventsList;
97 fEventsList=0 ;
98 }
591cc579 99
1c5acb87 100}
101
102//________________________________________________________________________________________________________________________________________________
103void AliAnaPi0::InitParameters()
104{
105//Init parameters when first called the analysis
106//Set default parameters
a3aebfff 107 SetInputAODName("PWG4Particle");
108
109 AddToHistogramsName("AnaPi0_");
6921fa00 110 fNModules = 12; // set maximum to maximum number of EMCAL modules
477d6cee 111 fNCentrBin = 1;
5025c139 112// fNZvertBin = 1;
113// fNrpBin = 1;
477d6cee 114 fNmaxMixEv = 10;
5025c139 115
477d6cee 116 fCalorimeter = "PHOS";
50f39b97 117 fUseAngleCut = kFALSE;
6175da48 118 fUseAngleEDepCut = kFALSE;
119 fAngleCut = 0.;
120 fAngleMaxCut = TMath::Pi();
121
5ae09196 122 fMultiCutAna = kFALSE;
123
124 fNPtCuts = 3;
af7b3903 125 fPtCuts[0] = 0.; fPtCuts[1] = 0.3; fPtCuts[2] = 0.5;
126 for(Int_t i = fNPtCuts; i < 10; i++)fPtCuts[i] = 0.;
5ae09196 127
9c59b5fe 128 fNAsymCuts = 4;
129 fAsymCuts[0] = 1.; fAsymCuts[1] = 0.8; fAsymCuts[2] = 0.6; fAsymCuts[3] = 0.1;
af7b3903 130 for(Int_t i = fNAsymCuts; i < 10; i++)fAsymCuts[i] = 0.;
131
5ae09196 132 fNCellNCuts = 3;
af7b3903 133 fCellNCuts[0] = 0; fCellNCuts[1] = 1; fCellNCuts[2] = 2;
021c573a 134 for(Int_t i = fNCellNCuts; i < 10; i++)fCellNCuts[i] = 0;
af7b3903 135
136 fNPIDBits = 2;
137 fPIDBits[0] = 0; fPIDBits[1] = 2; // fPIDBits[2] = 4; fPIDBits[3] = 6;// check, no cut, dispersion, neutral, dispersion&&neutral
021c573a 138 for(Int_t i = fNPIDBits; i < 10; i++)fPIDBits[i] = 0;
af7b3903 139
1c5acb87 140}
1c5acb87 141
0c1383b5 142
143//________________________________________________________________________________________________________________________________________________
144TObjString * AliAnaPi0::GetAnalysisCuts()
145{
af7b3903 146 //Save parameters used for analysis
147 TString parList ; //this will be list of parameters used for this analysis.
148 const Int_t buffersize = 255;
149 char onePar[buffersize] ;
150 snprintf(onePar,buffersize,"--- AliAnaPi0 ---\n") ;
151 parList+=onePar ;
152 snprintf(onePar,buffersize,"Number of bins in Centrality: %d \n",fNCentrBin) ;
153 parList+=onePar ;
154 snprintf(onePar,buffersize,"Number of bins in Z vert. pos: %d \n",GetNZvertBin()) ;
155 parList+=onePar ;
156 snprintf(onePar,buffersize,"Number of bins in Reac. Plain: %d \n",GetNRPBin()) ;
157 parList+=onePar ;
158 snprintf(onePar,buffersize,"Depth of event buffer: %d \n",fNmaxMixEv) ;
159 parList+=onePar ;
6175da48 160 snprintf(onePar,buffersize,"Pair in same Module: %d ; TrackMult as centrality: %d; PhotonMult as centrality: %d; cluster E as centrality: %d; cell as centrality: %d; Fill InvPt histos %d\n",
161 fSameSM, fUseTrackMultBins, fUsePhotonMultBins, fUseAverClusterEBins, fUseAverCellEBins, fMakeInvPtPlots) ;
162 parList+=onePar ;
163 snprintf(onePar,buffersize,"Select pairs with their angle: %d, edep %d, min angle %2.3f, max angle %2.3f,\n",fUseAngleCut, fUseAngleEDepCut,fAngleCut,fAngleMaxCut) ;
af7b3903 164 parList+=onePar ;
165 snprintf(onePar,buffersize," Asymmetry cuts: n = %d, asymmetry < ",fNAsymCuts) ;
166 for(Int_t i = 0; i < fNAsymCuts; i++) snprintf(onePar,buffersize,"%s %2.2f;",onePar,fAsymCuts[i]);
167 parList+=onePar ;
168 snprintf(onePar,buffersize," PID selection bits: n = %d, PID bit =\n",fNPIDBits) ;
169 for(Int_t i = 0; i < fNPIDBits; i++) snprintf(onePar,buffersize,"%s %d;",onePar,fPIDBits[i]);
170 parList+=onePar ;
171 snprintf(onePar,buffersize,"Cuts: \n") ;
172 parList+=onePar ;
173 snprintf(onePar,buffersize,"Z vertex position: -%f < z < %f \n",GetZvertexCut(),GetZvertexCut()) ;
174 parList+=onePar ;
175 snprintf(onePar,buffersize,"Calorimeter: %s \n",fCalorimeter.Data()) ;
176 parList+=onePar ;
177 snprintf(onePar,buffersize,"Number of modules: %d \n",fNModules) ;
178 parList+=onePar ;
db2bf6fd 179 if(fMultiCutAna){
180 snprintf(onePar, buffersize," pT cuts: n = %d, pt > ",fNPtCuts) ;
181 for(Int_t i = 0; i < fNPtCuts; i++) snprintf(onePar,buffersize,"%s %2.2f;",onePar,fPtCuts[i]);
182 parList+=onePar ;
db2bf6fd 183 snprintf(onePar,buffersize, " N cell in cluster cuts: n = %d, nCell > ",fNCellNCuts) ;
184 for(Int_t i = 0; i < fNCellNCuts; i++) snprintf(onePar,buffersize,"%s %d;",onePar,fCellNCuts[i]);
185 parList+=onePar ;
db2bf6fd 186 }
187
af7b3903 188 return new TObjString(parList) ;
0c1383b5 189}
190
2e557d1c 191//________________________________________________________________________________________________________________________________________________
192TList * AliAnaPi0::GetCreateOutputObjects()
193{
477d6cee 194 // Create histograms to be saved in output file and
195 // store them in fOutputContainer
196
197 //create event containers
5025c139 198 fEventsList = new TList*[fNCentrBin*GetNZvertBin()*GetNRPBin()] ;
1c5acb87 199
477d6cee 200 for(Int_t ic=0; ic<fNCentrBin; ic++){
5025c139 201 for(Int_t iz=0; iz<GetNZvertBin(); iz++){
202 for(Int_t irp=0; irp<GetNRPBin(); irp++){
203 fEventsList[ic*GetNZvertBin()*GetNRPBin()+iz*GetNRPBin()+irp] = new TList() ;
af7b3903 204 fEventsList[ic*GetNZvertBin()*GetNRPBin()+iz*GetNRPBin()+irp]->SetOwner(kFALSE);
477d6cee 205 }
206 }
207 }
7e7694bb 208
477d6cee 209 TList * outputContainer = new TList() ;
210 outputContainer->SetName(GetName());
6921fa00 211
af7b3903 212 fhReMod = new TH2D*[fNModules] ;
6175da48 213 fhReDiffMod = new TH2D*[fNModules+3] ;
214
215 fhMiMod = new TH2D*[fNModules] ;
216 fhMiDiffMod = new TH2D*[fNModules+3] ;
af7b3903 217
218 fhRe1 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
af7b3903 219 fhMi1 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
6175da48 220 if(fFillBadDistHisto){
221 fhRe2 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
222 fhRe3 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
223 fhMi2 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
224 fhMi3 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
225 }
226 if(fMakeInvPtPlots) {
398c93cc 227 fhReInvPt1 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
398c93cc 228 fhMiInvPt1 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
6175da48 229 if(fFillBadDistHisto){
230 fhReInvPt2 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
231 fhReInvPt3 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
232 fhMiInvPt2 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
233 fhMiInvPt3 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
234 }
398c93cc 235 }
6175da48 236
5ae09196 237 const Int_t buffersize = 255;
238 char key[buffersize] ;
239 char title[buffersize] ;
477d6cee 240
5a2dbc3c 241 Int_t nptbins = GetHistoPtBins();
242 Int_t nphibins = GetHistoPhiBins();
243 Int_t netabins = GetHistoEtaBins();
244 Float_t ptmax = GetHistoPtMax();
245 Float_t phimax = GetHistoPhiMax();
246 Float_t etamax = GetHistoEtaMax();
247 Float_t ptmin = GetHistoPtMin();
248 Float_t phimin = GetHistoPhiMin();
249 Float_t etamin = GetHistoEtaMin();
250
251 Int_t nmassbins = GetHistoMassBins();
252 Int_t nasymbins = GetHistoAsymmetryBins();
253 Float_t massmax = GetHistoMassMax();
254 Float_t asymmax = GetHistoAsymmetryMax();
255 Float_t massmin = GetHistoMassMin();
256 Float_t asymmin = GetHistoAsymmetryMin();
af7b3903 257 Int_t ntrmbins = GetHistoTrackMultiplicityBins();
258 Int_t ntrmmax = GetHistoTrackMultiplicityMax();
259 Int_t ntrmmin = GetHistoTrackMultiplicityMin();
6175da48 260
261 fhAverTotECluster = new TH1F("hAverTotECluster","hAverTotECluster",200,0,50) ;
262 fhAverTotECluster->SetXTitle("E_{cluster, aver. SM} (GeV)");
263 outputContainer->Add(fhAverTotECluster) ;
264
265 fhAverTotECell = new TH1F("hAverTotECell","hAverTotECell",200,0,50) ;
266 fhAverTotECell->SetXTitle("E_{cell, aver. SM} (GeV)");
267 outputContainer->Add(fhAverTotECell) ;
268
269 fhReConv = new TH2D("hReConv","Real Pair with one recombined conversion ",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
270 fhReConv->SetXTitle("p_{T} (GeV/c)");
271 fhReConv->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
272 outputContainer->Add(fhReConv) ;
273
274 fhMiConv = new TH2D("hMiConv","Mixed Pair with one recombined conversion ",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
275 fhMiConv->SetXTitle("p_{T} (GeV/c)");
276 fhMiConv->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
277 outputContainer->Add(fhMiConv) ;
278
279 fhReConv2 = new TH2D("hReConv2","Real Pair with 2 recombined conversion ",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
280 fhReConv2->SetXTitle("p_{T} (GeV/c)");
281 fhReConv2->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
282 outputContainer->Add(fhReConv2) ;
283
284 fhMiConv2 = new TH2D("hMiConv2","Mixed Pair with 2 recombined conversion ",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
285 fhMiConv2->SetXTitle("p_{T} (GeV/c)");
286 fhMiConv2->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
287 outputContainer->Add(fhMiConv2) ;
288
477d6cee 289 for(Int_t ic=0; ic<fNCentrBin; ic++){
6175da48 290 for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
291 for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
292 Int_t index = ((ic*fNPIDBits)+ipid)*fNAsymCuts + iasym;
293 //printf("cen %d, pid %d, asy %d, Index %d\n",ic,ipid,iasym,index);
398c93cc 294 //Distance to bad module 1
6175da48 295 snprintf(key, buffersize,"hRe_cen%d_pidbit%d_asy%d_dist1",ic,ipid,iasym) ;
398c93cc 296 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 1",
297 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
6175da48 298 fhRe1[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
299 fhRe1[index]->SetXTitle("p_{T} (GeV/c)");
300 fhRe1[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
301 //printf("name: %s\n ",fhRe1[index]->GetName());
302 outputContainer->Add(fhRe1[index]) ;
af7b3903 303
6175da48 304 if(fFillBadDistHisto){
398c93cc 305 //Distance to bad module 2
6175da48 306 snprintf(key, buffersize,"hRe_cen%d_pidbit%d_asy%d_dist2",ic,ipid,iasym) ;
307 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 2",
398c93cc 308 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
6175da48 309 fhRe2[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
310 fhRe2[index]->SetXTitle("p_{T} (GeV/c)");
311 fhRe2[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
312 outputContainer->Add(fhRe2[index]) ;
398c93cc 313
314 //Distance to bad module 3
6175da48 315 snprintf(key, buffersize,"hRe_cen%d_pidbit%d_asy%d_dist3",ic,ipid,iasym) ;
316 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 3",
398c93cc 317 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
6175da48 318 fhRe3[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
319 fhRe3[index]->SetXTitle("p_{T} (GeV/c)");
320 fhRe3[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
321 outputContainer->Add(fhRe3[index]) ;
398c93cc 322 }
6175da48 323
324 //Inverse pT
325 if(fMakeInvPtPlots){
326 //Distance to bad module 1
327 snprintf(key, buffersize,"hReInvPt_cen%d_pidbit%d_asy%d_dist1",ic,ipid,iasym) ;
328 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 1",
329 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
330 fhReInvPt1[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
331 fhReInvPt1[index]->SetXTitle("p_{T} (GeV/c)");
332 fhReInvPt1[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
333 outputContainer->Add(fhReInvPt1[index]) ;
334
335 if(fFillBadDistHisto){
336 //Distance to bad module 2
337 snprintf(key, buffersize,"hReInvPt_cen%d_pidbit%d_asy%d_dist2",ic,ipid,iasym) ;
338 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 2",
339 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
340 fhReInvPt2[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
341 fhReInvPt2[index]->SetXTitle("p_{T} (GeV/c)");
342 fhReInvPt2[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
343 outputContainer->Add(fhReInvPt2[index]) ;
344
345 //Distance to bad module 3
346 snprintf(key, buffersize,"hReInvPt_cen%d_pidbit%d_asy%d_dist3",ic,ipid,iasym) ;
347 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 3",
348 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
349 fhReInvPt3[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
350 fhReInvPt3[index]->SetXTitle("p_{T} (GeV/c)");
351 fhReInvPt3[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
352 outputContainer->Add(fhReInvPt3[index]) ;
353 }
354 }
355 if(fDoOwnMix){
356 //Distance to bad module 1
357 snprintf(key, buffersize,"hMi_cen%d_pidbit%d_asy%d_dist1",ic,ipid,iasym) ;
358 snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 1",
359 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
360 fhMi1[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
361 fhMi1[index]->SetXTitle("p_{T} (GeV/c)");
362 fhMi1[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
363 outputContainer->Add(fhMi1[index]) ;
364 if(fFillBadDistHisto){
365 //Distance to bad module 2
366 snprintf(key, buffersize,"hMi_cen%d_pidbit%d_asy%d_dist2",ic,ipid,iasym) ;
367 snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 2",
368 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
369 fhMi2[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
370 fhMi2[index]->SetXTitle("p_{T} (GeV/c)");
371 fhMi2[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
372 outputContainer->Add(fhMi2[index]) ;
373
374 //Distance to bad module 3
375 snprintf(key, buffersize,"hMi_cen%d_pidbit%d_asy%d_dist3",ic,ipid,iasym) ;
376 snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 3",
377 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
378 fhMi3[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
379 fhMi3[index]->SetXTitle("p_{T} (GeV/c)");
380 fhMi3[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
381 outputContainer->Add(fhMi3[index]) ;
382 }
383 //Inverse pT
384 if(fMakeInvPtPlots){
385 //Distance to bad module 1
386 snprintf(key, buffersize,"hMiInvPt_cen%d_pidbit%d_asy%d_dist1",ic,ipid,iasym) ;
387 snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 1",
388 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
389 fhMiInvPt1[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
390 fhMiInvPt1[index]->SetXTitle("p_{T} (GeV/c)");
391 fhMiInvPt1[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
392 outputContainer->Add(fhMiInvPt1[index]) ;
393 if(fFillBadDistHisto){
394 //Distance to bad module 2
395 snprintf(key, buffersize,"hMiInvPt_cen%d_pidbit%d_asy%d_dist2",ic,ipid,iasym) ;
396 snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 2",
397 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
398 fhMiInvPt2[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
399 fhMiInvPt2[index]->SetXTitle("p_{T} (GeV/c)");
400 fhMiInvPt2[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
401 outputContainer->Add(fhMiInvPt2[index]) ;
402
403 //Distance to bad module 3
404 snprintf(key, buffersize,"hMiInvPt_cen%d_pidbit%d_asy%d_dist3",ic,ipid,iasym) ;
405 snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f,dist bad 3",
406 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
407 fhMiInvPt3[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
408 fhMiInvPt3[index]->SetXTitle("p_{T} (GeV/c)");
409 fhMiInvPt3[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
410 outputContainer->Add(fhMiInvPt3[index]) ;
411 }
412 }
413 }
414 }
7e7694bb 415 }
1c5acb87 416 }
477d6cee 417
9c59b5fe 418 fhRePtAsym = new TH2D("hRePtAsym","Asymmetry vs pt, for pairs",nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax) ;
af7b3903 419 fhRePtAsym->SetXTitle("p_{T} (GeV/c)");
420 fhRePtAsym->SetYTitle("Asymmetry");
421 outputContainer->Add(fhRePtAsym);
422
9c59b5fe 423 fhRePtAsymPi0 = new TH2D("hRePtAsymPi0","Asymmetry vs pt, for pairs close to #pi^{0} mass",nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax) ;
af7b3903 424 fhRePtAsymPi0->SetXTitle("p_{T} (GeV/c)");
425 fhRePtAsymPi0->SetYTitle("Asymmetry");
426 outputContainer->Add(fhRePtAsymPi0);
427
9c59b5fe 428 fhRePtAsymEta = new TH2D("hRePtAsymEta","Asymmetry vs pt, for pairs close to #eta mass",nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax) ;
af7b3903 429 fhRePtAsymEta->SetXTitle("p_{T} (GeV/c)");
430 fhRePtAsymEta->SetYTitle("Asymmetry");
431 outputContainer->Add(fhRePtAsymEta);
432
5ae09196 433 if(fMultiCutAna){
434
435 fhRePIDBits = new TH2D*[fNPIDBits];
436 for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
437 snprintf(key, buffersize,"hRe_pidbit%d",ipid) ;
438 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for PIDBit=%d",fPIDBits[ipid]) ;
439 fhRePIDBits[ipid] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
af7b3903 440 fhRePIDBits[ipid]->SetXTitle("p_{T} (GeV/c)");
441 fhRePIDBits[ipid]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
5ae09196 442 outputContainer->Add(fhRePIDBits[ipid]) ;
443 }// pid bit loop
444
6175da48 445 fhRePtNCellAsymCuts = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
446 fhRePtNCellAsymCutsSM0 = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
447 fhRePtNCellAsymCutsSM1 = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
448 fhRePtNCellAsymCutsSM2 = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
449 fhRePtNCellAsymCutsSM3 = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
450 fhMiPtNCellAsymCuts = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
5ae09196 451 for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
452 for(Int_t icell=0; icell<fNCellNCuts; icell++){
453 for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
454 snprintf(key, buffersize,"hRe_pt%d_cell%d_asym%d",ipt,icell,iasym) ;
af7b3903 455 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for pt >%2.2f, ncell>%d and asym >%1.2f ",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]) ;
5ae09196 456 Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
457 //printf("ipt %d, icell %d, iassym %d, index %d\n",ipt, icell, iasym, index);
458 fhRePtNCellAsymCuts[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
af7b3903 459 fhRePtNCellAsymCuts[index]->SetXTitle("p_{T} (GeV/c)");
460 fhRePtNCellAsymCuts[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
5ae09196 461 outputContainer->Add(fhRePtNCellAsymCuts[index]) ;
6175da48 462
463 snprintf(key, buffersize,"hRe_pt%d_cell%d_asym%d_SM0",ipt,icell,iasym) ;
464 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for pt >%2.2f, ncell>%d and asym >%1.2f, SM 0 ",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]) ;
465 fhRePtNCellAsymCutsSM0[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
466 fhRePtNCellAsymCutsSM0[index]->SetXTitle("p_{T} (GeV/c)");
467 fhRePtNCellAsymCutsSM0[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
468 outputContainer->Add(fhRePtNCellAsymCutsSM0[index]) ;
469
470 snprintf(key, buffersize,"hRe_pt%d_cell%d_asym%d_SM1",ipt,icell,iasym) ;
471 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for pt >%2.2f, ncell>%d and asym >%1.2f, SM 1 ",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]) ;
472 fhRePtNCellAsymCutsSM1[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
473 fhRePtNCellAsymCutsSM1[index]->SetXTitle("p_{T} (GeV/c)");
474 fhRePtNCellAsymCutsSM1[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
475 outputContainer->Add(fhRePtNCellAsymCutsSM1[index]) ;
476
477 snprintf(key, buffersize,"hRe_pt%d_cell%d_asym%d_SM2",ipt,icell,iasym) ;
478 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for pt >%2.2f, ncell>%d and asym >%1.2f, SM 2 ",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]) ;
479 fhRePtNCellAsymCutsSM2[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
480 fhRePtNCellAsymCutsSM2[index]->SetXTitle("p_{T} (GeV/c)");
481 fhRePtNCellAsymCutsSM2[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
482 outputContainer->Add(fhRePtNCellAsymCutsSM2[index]) ;
483
484 snprintf(key, buffersize,"hRe_pt%d_cell%d_asym%d_SM3",ipt,icell,iasym) ;
485 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for pt >%2.2f, ncell>%d and asym >%1.2f, SM 3 ",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]) ;
486 fhRePtNCellAsymCutsSM3[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
487 fhRePtNCellAsymCutsSM3[index]->SetXTitle("p_{T} (GeV/c)");
488 fhRePtNCellAsymCutsSM3[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
489 outputContainer->Add(fhRePtNCellAsymCutsSM3[index]) ;
490
491 snprintf(key, buffersize,"hMi_pt%d_cell%d_asym%d",ipt,icell,iasym) ;
492 snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for pt >%2.2f, ncell>%d and asym >%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]) ;
493 fhMiPtNCellAsymCuts[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
494 fhMiPtNCellAsymCuts[index]->SetXTitle("p_{T} (GeV/c)");
495 fhMiPtNCellAsymCuts[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
496 outputContainer->Add(fhMiPtNCellAsymCuts[index]) ;
497
5ae09196 498 }
499 }
500 }
821c8090 501
af7b3903 502 fhRePtMult = new TH3D*[fNAsymCuts] ;
503 for(Int_t iasym = 0; iasym<fNAsymCuts; iasym++){
504 fhRePtMult[iasym] = new TH3D(Form("hRePtMult_asym%d",iasym),Form("(p_{T},C,M)_{#gamma#gamma}, A<%1.2f",fAsymCuts[iasym]),
505 nptbins,ptmin,ptmax,ntrmbins,ntrmmin,ntrmmax,nmassbins,massmin,massmax);
506 fhRePtMult[iasym]->SetXTitle("p_{T} (GeV/c)");
507 fhRePtMult[iasym]->SetYTitle("Track multiplicity");
508 fhRePtMult[iasym]->SetZTitle("m_{#gamma,#gamma} (GeV/c^{2})");
509 outputContainer->Add(fhRePtMult[iasym]) ;
510 }
511
5ae09196 512 }// multi cuts analysis
513
477d6cee 514 fhEvents=new TH3D("hEvents","Number of events",fNCentrBin,0.,1.*fNCentrBin,
5025c139 515 GetNZvertBin(),0.,1.*GetNZvertBin(),GetNRPBin(),0.,1.*GetNRPBin()) ;
477d6cee 516 outputContainer->Add(fhEvents) ;
50f39b97 517
518 fhRealOpeningAngle = new TH2D
6175da48 519 ("hRealOpeningAngle","Angle between all #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,300,0,TMath::Pi());
50f39b97 520 fhRealOpeningAngle->SetYTitle("#theta(rad)");
521 fhRealOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
522 outputContainer->Add(fhRealOpeningAngle) ;
7e7694bb 523
50f39b97 524 fhRealCosOpeningAngle = new TH2D
6175da48 525 ("hRealCosOpeningAngle","Cosinus of angle between all #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,100,0,1);
50f39b97 526 fhRealCosOpeningAngle->SetYTitle("cos (#theta) ");
527 fhRealCosOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
528 outputContainer->Add(fhRealCosOpeningAngle) ;
529
6175da48 530 if(fDoOwnMix){
531
532 fhMixedOpeningAngle = new TH2D
533 ("hMixedOpeningAngle","Angle between all #gamma pair vs E_{#pi^{0}}, Mixed pairs",nptbins,ptmin,ptmax,300,0,TMath::Pi());
534 fhMixedOpeningAngle->SetYTitle("#theta(rad)");
535 fhMixedOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
536 outputContainer->Add(fhMixedOpeningAngle) ;
537
538 fhMixedCosOpeningAngle = new TH2D
539 ("hMixedCosOpeningAngle","Cosinus of angle between all #gamma pair vs E_{#pi^{0}}, Mixed pairs",nptbins,ptmin,ptmax,100,0,1);
540 fhMixedCosOpeningAngle->SetYTitle("cos (#theta) ");
541 fhMixedCosOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
542 outputContainer->Add(fhMixedCosOpeningAngle) ;
543
544 }
545
477d6cee 546 //Histograms filled only if MC data is requested
0ae57829 547 if(IsDataMC()){
6175da48 548 //Pi0
549 fhPrimPi0Pt = new TH1D("hPrimPi0Pt","Primary pi0 pt",nptbins,ptmin,ptmax) ;
550 fhPrimPi0AccPt = new TH1D("hPrimPi0AccPt","Primary pi0 pt with both photons in acceptance",nptbins,ptmin,ptmax) ;
551 outputContainer->Add(fhPrimPi0Pt) ;
552 outputContainer->Add(fhPrimPi0AccPt) ;
553
554 fhPrimPi0Y = new TH1D("hPrimPi0Rapidity","Rapidity of primary pi0",netabins,etamin,etamax) ;
555 outputContainer->Add(fhPrimPi0Y) ;
556
557 fhPrimPi0AccY = new TH1D("hPrimPi0AccRapidity","Rapidity of primary pi0",netabins,etamin,etamax) ;
558 outputContainer->Add(fhPrimPi0AccY) ;
477d6cee 559
6175da48 560 fhPrimPi0Phi = new TH1D("hPrimPi0Phi","Azimithal of primary pi0",nphibins,phimin*TMath::RadToDeg(),phimax*TMath::RadToDeg()) ;
561 outputContainer->Add(fhPrimPi0Phi) ;
477d6cee 562
6175da48 563 fhPrimPi0AccPhi = new TH1D("hPrimPi0AccPhi","Azimithal of primary pi0 with accepted daughters",nphibins,phimin*TMath::RadToDeg(),phimax*TMath::RadToDeg()) ;
564 outputContainer->Add(fhPrimPi0AccPhi) ;
477d6cee 565
6175da48 566 //Eta
567 fhPrimEtaPt = new TH1D("hPrimEtaPt","Primary eta pt",nptbins,ptmin,ptmax) ;
568 fhPrimEtaAccPt = new TH1D("hPrimEtaAccPt","Primary eta pt with both photons in acceptance",nptbins,ptmin,ptmax) ;
569 outputContainer->Add(fhPrimEtaPt) ;
570 outputContainer->Add(fhPrimEtaAccPt) ;
477d6cee 571
6175da48 572 fhPrimEtaY = new TH1D("hPrimEtaRapidity","Rapidity of primary eta",netabins,etamin,etamax) ;
573 outputContainer->Add(fhPrimEtaY) ;
50f39b97 574
6175da48 575 fhPrimEtaAccY = new TH1D("hPrimEtaAccRapidity","Rapidity of primary eta",netabins,etamin,etamax) ;
576 outputContainer->Add(fhPrimEtaAccY) ;
50f39b97 577
6175da48 578 fhPrimEtaPhi = new TH1D("hPrimEtaPhi","Azimithal of primary eta",nphibins,phimin*TMath::RadToDeg(),phimax*TMath::RadToDeg()) ;
579 outputContainer->Add(fhPrimEtaPhi) ;
580
581 fhPrimEtaAccPhi = new TH1D("hPrimEtaAccPhi","Azimithal of primary eta with accepted daughters",nphibins,phimin*TMath::RadToDeg(),phimax*TMath::RadToDeg()) ;
582 outputContainer->Add(fhPrimEtaAccPhi) ;
583
50f39b97 584
6175da48 585 fhPrimPi0OpeningAngle = new TH2D
586 ("hPrimPi0OpeningAngle","Angle between all primary #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,100,0,0.5);
587 fhPrimPi0OpeningAngle->SetYTitle("#theta(rad)");
588 fhPrimPi0OpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
589 outputContainer->Add(fhPrimPi0OpeningAngle) ;
590
591 fhPrimPi0CosOpeningAngle = new TH2D
592 ("hPrimPi0CosOpeningAngle","Cosinus of angle between all primary #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,100,-1,1);
593 fhPrimPi0CosOpeningAngle->SetYTitle("cos (#theta) ");
594 fhPrimPi0CosOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
595 outputContainer->Add(fhPrimPi0CosOpeningAngle) ;
596
597 for(Int_t i = 0; i<13; i++){
598 fhMCOrgMass[i] = new TH2D(Form("hMCOrgMass_%d",i),Form("mass vs pt, origin %d",i),nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
599 fhMCOrgMass[i]->SetXTitle("p_{T} (GeV/c)");
600 fhMCOrgMass[i]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
601 outputContainer->Add(fhMCOrgMass[i]) ;
602
603 fhMCOrgAsym[i]= new TH2D(Form("hMCOrgAsym_%d",i),Form("asymmetry vs pt, origin %d",i),nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax) ;
604 fhMCOrgAsym[i]->SetXTitle("p_{T} (GeV/c)");
605 fhMCOrgAsym[i]->SetYTitle("A");
606 outputContainer->Add(fhMCOrgAsym[i]) ;
607
608 fhMCOrgDeltaEta[i] = new TH2D(Form("hMCOrgDeltaEta_%d",i),Form("mass vs pt, origin %d",i),nptbins,ptmin,ptmax,netabins,-1.4,1.4) ;
609 fhMCOrgDeltaEta[i]->SetXTitle("p_{T} (GeV/c)");
610 fhMCOrgDeltaEta[i]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
611 outputContainer->Add(fhMCOrgDeltaEta[i]) ;
612
613 fhMCOrgDeltaPhi[i]= new TH2D(Form("hMCOrgDeltaPhi_%d",i),Form("asymmetry vs pt, origin %d",i),nptbins,ptmin,ptmax,nphibins,-0.7,0.7) ;
614 fhMCOrgDeltaPhi[i]->SetXTitle("p_{T} (GeV/c)");
615 fhMCOrgDeltaPhi[i]->SetYTitle("A");
616 outputContainer->Add(fhMCOrgDeltaPhi[i]) ;
617
618 }
50f39b97 619
6175da48 620 if(fMultiCutAnaSim){
621 fhMCPi0MassPtTrue = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
622 fhMCPi0MassPtRec = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
623 fhMCPi0PtTruePtRec = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
624 fhMCEtaMassPtRec = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
625 fhMCEtaMassPtTrue = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
626 fhMCEtaPtTruePtRec = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
627 for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
628 for(Int_t icell=0; icell<fNCellNCuts; icell++){
629 for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
630 Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
631
632 fhMCPi0MassPtRec[index] = new TH2D(Form("hMCPi0MassPtRec_pt%d_cell%d_asym%d",ipt,icell,iasym),
633 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]),
634 nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
635 fhMCPi0MassPtRec[index]->SetXTitle("p_{T, reconstructed} (GeV/c)");
636 fhMCPi0MassPtRec[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
637 outputContainer->Add(fhMCPi0MassPtRec[index]) ;
638
639 fhMCPi0MassPtTrue[index] = new TH2D(Form("hMCPi0MassPtTrue_pt%d_cell%d_asym%d",ipt,icell,iasym),
640 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]),
641 nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
642 fhMCPi0MassPtTrue[index]->SetXTitle("p_{T, generated} (GeV/c)");
643 fhMCPi0MassPtTrue[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
644 outputContainer->Add(fhMCPi0MassPtTrue[index]) ;
645
646 fhMCPi0PtTruePtRec[index] = new TH2D(Form("hMCPi0PtTruePtRec_pt%d_cell%d_asym%d",ipt,icell,iasym),
647 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]),
648 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax) ;
649 fhMCPi0PtTruePtRec[index]->SetXTitle("p_{T, generated} (GeV/c)");
650 fhMCPi0PtTruePtRec[index]->SetYTitle("p_{T, reconstructed} (GeV/c)");
651 outputContainer->Add(fhMCPi0PtTruePtRec[index]) ;
652
653 fhMCEtaMassPtRec[index] = new TH2D(Form("hMCEtaMassPtRec_pt%d_cell%d_asym%d",ipt,icell,iasym),
654 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]),
655 nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
656 fhMCEtaMassPtRec[index]->SetXTitle("p_{T, generated} (GeV/c)");
657 fhMCEtaMassPtRec[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
658 outputContainer->Add(fhMCEtaMassPtRec[index]) ;
659
660 fhMCEtaMassPtTrue[index] = new TH2D(Form("hMCEtaMassPtTrue_pt%d_cell%d_asym%d",ipt,icell,iasym),
661 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]),
662 nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
663 fhMCEtaMassPtTrue[index]->SetXTitle("p_{T, generated} (GeV/c)");
664 fhMCEtaMassPtTrue[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
665 outputContainer->Add(fhMCEtaMassPtTrue[index]) ;
666
667 fhMCEtaPtTruePtRec[index] = new TH2D(Form("hMCEtaPtTruePtRec_pt%d_cell%d_asym%d",ipt,icell,iasym),
668 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]),
669 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax) ;
670 fhMCEtaPtTruePtRec[index]->SetXTitle("p_{T, generated} (GeV/c)");
671 fhMCEtaPtTruePtRec[index]->SetYTitle("p_{T, reconstructed} (GeV/c)");
672 outputContainer->Add(fhMCEtaPtTruePtRec[index]) ;
673 }
674 }
675 }
676 }//multi cut ana
677 else {
678 fhMCPi0MassPtTrue = new TH2D*[1];
679 fhMCPi0PtTruePtRec = new TH2D*[1];
680 fhMCEtaMassPtTrue = new TH2D*[1];
681 fhMCEtaPtTruePtRec = new TH2D*[1];
682
683 fhMCPi0MassPtTrue[0] = new TH2D("hMCPi0MassPtTrue","Reconstructed Mass vs generated p_T of true #pi^{0} cluster pairs",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
684 fhMCPi0MassPtTrue[0]->SetXTitle("p_{T, generated} (GeV/c)");
685 fhMCPi0MassPtTrue[0]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
686 outputContainer->Add(fhMCPi0MassPtTrue[0]) ;
687
688 fhMCPi0PtTruePtRec[0]= new TH2D("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) ;
689 fhMCPi0PtTruePtRec[0]->SetXTitle("p_{T, generated} (GeV/c)");
690 fhMCPi0PtTruePtRec[0]->SetYTitle("p_{T, reconstructed} (GeV/c)");
691 outputContainer->Add(fhMCPi0PtTruePtRec[0]) ;
692
693 fhMCEtaMassPtTrue[0] = new TH2D("hMCEtaMassPtTrue","Reconstructed Mass vs generated p_T of true #eta cluster pairs",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
694 fhMCEtaMassPtTrue[0]->SetXTitle("p_{T, generated} (GeV/c)");
695 fhMCEtaMassPtTrue[0]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
696 outputContainer->Add(fhMCEtaMassPtTrue[0]) ;
697
698 fhMCEtaPtTruePtRec[0]= new TH2D("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) ;
699 fhMCEtaPtTruePtRec[0]->SetXTitle("p_{T, generated} (GeV/c)");
700 fhMCEtaPtTruePtRec[0]->SetYTitle("p_{T, reconstructed} (GeV/c)");
701 outputContainer->Add(fhMCEtaPtTruePtRec[0]) ;
702 }
477d6cee 703 }
50f39b97 704
6175da48 705 TString * pairname = new TString[fNModules+3];
821c8090 706 if(fCalorimeter=="EMCAL"){
707 pairname[0]="A side (0-2)";
708 pairname[1]="C side (1-3)";
709 pairname[2]="Sector 0 (0-1)";
710 pairname[3]="Sector 1 (2-3)";
6175da48 711 pairname[4]="Cluster in different SM";
712 pairname[5]="SM 0 and SM3";
713 pairname[6]="SM 1 and SM2";
714 for(Int_t i = 7 ; i < fNModules ; i++) pairname[i]="";}
821c8090 715 if(fCalorimeter=="PHOS") {
716 pairname[0]="(0-1)";
717 pairname[1]="(0-2)";
af7b3903 718 pairname[2]="(1-2)";
821c8090 719 for(Int_t i = 3 ; i < fNModules ; i++) pairname[i]="";}
720
6921fa00 721 for(Int_t imod=0; imod<fNModules; imod++){
50f39b97 722 //Module dependent invariant mass
5ae09196 723 snprintf(key, buffersize,"hReMod_%d",imod) ;
724 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for Module %d",imod) ;
af7b3903 725 fhReMod[imod] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
726 fhReMod[imod]->SetXTitle("p_{T} (GeV/c)");
727 fhReMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
50f39b97 728 outputContainer->Add(fhReMod[imod]) ;
821c8090 729
730 snprintf(key, buffersize,"hReDiffMod_%d",imod) ;
731 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for Different Modules: %s",(pairname[imod]).Data()) ;
af7b3903 732 fhReDiffMod[imod] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
733 fhReDiffMod[imod]->SetXTitle("p_{T} (GeV/c)");
734 fhReDiffMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
821c8090 735 outputContainer->Add(fhReDiffMod[imod]) ;
6175da48 736
737 if(fDoOwnMix){
738 snprintf(key, buffersize,"hMiMod_%d",imod) ;
739 snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for Module %d",imod) ;
740 fhMiMod[imod] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
741 fhMiMod[imod]->SetXTitle("p_{T} (GeV/c)");
742 fhMiMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
743 outputContainer->Add(fhMiMod[imod]) ;
744
745 snprintf(key, buffersize,"hMiDiffMod_%d",imod) ;
746 snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for Different Modules: %s",(pairname[imod]).Data()) ;
747 fhMiDiffMod[imod] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
748 fhMiDiffMod[imod]->SetXTitle("p_{T} (GeV/c)");
749 fhMiDiffMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
750 outputContainer->Add(fhMiDiffMod[imod]) ;
751 }
752
6921fa00 753 }
50f39b97 754
6175da48 755 for (Int_t imod=4; imod<7; imod++) {
756
757 snprintf(key, buffersize,"hReDiffMod_%d",imod) ;
758 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for Different Modules: %s",(pairname[imod]).Data()) ;
759 fhReDiffMod[imod] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
760 fhReDiffMod[imod]->SetXTitle("p_{T} (GeV/c)");
761 fhReDiffMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
762 outputContainer->Add(fhReDiffMod[imod]) ;
821c8090 763
6175da48 764 if(fDoOwnMix){
765 snprintf(key, buffersize,"hMiDiffMod_%d",imod) ;
766 snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for Different Modules: %s",(pairname[imod]).Data()) ;
767 fhMiDiffMod[imod] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
768 fhMiDiffMod[imod]->SetXTitle("p_{T} (GeV/c)");
769 fhMiDiffMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
770 outputContainer->Add(fhMiDiffMod[imod]) ;
771 }
772 }
821c8090 773
6175da48 774 delete [] pairname;
821c8090 775
eee5fcf1 776// for(Int_t i = 0; i < outputContainer->GetEntries() ; i++){
777//
778// printf("Histogram %d, name: %s\n ",i, outputContainer->At(i)->GetName());
779//
780// }
781
477d6cee 782 return outputContainer;
1c5acb87 783}
784
785//_________________________________________________________________________________________________________________________________________________
786void AliAnaPi0::Print(const Option_t * /*opt*/) const
787{
477d6cee 788 //Print some relevant parameters set for the analysis
a3aebfff 789 printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
477d6cee 790 AliAnaPartCorrBaseClass::Print(" ");
a3aebfff 791
477d6cee 792 printf("Number of bins in Centrality: %d \n",fNCentrBin) ;
5025c139 793 printf("Number of bins in Z vert. pos: %d \n",GetNZvertBin()) ;
794 printf("Number of bins in Reac. Plain: %d \n",GetNRPBin()) ;
477d6cee 795 printf("Depth of event buffer: %d \n",fNmaxMixEv) ;
af7b3903 796 printf("Pair in same Module: %d \n",fSameSM) ;
477d6cee 797 printf("Cuts: \n") ;
5025c139 798 printf("Z vertex position: -%2.3f < z < %2.3f \n",GetZvertexCut(),GetZvertexCut()) ;
50f39b97 799 printf("Number of modules: %d \n",fNModules) ;
6175da48 800 printf("Select pairs with their angle: %d, edep %d, min angle %2.3f, max angle %2.3f \n",fUseAngleCut, fUseAngleEDepCut, fAngleCut, fAngleMaxCut) ;
af7b3903 801 printf("Asymmetry cuts: n = %d, \n",fNAsymCuts) ;
802 printf("\tasymmetry < ");
803 for(Int_t i = 0; i < fNAsymCuts; i++) printf("%2.2f ",fAsymCuts[i]);
804 printf("\n");
805
806 printf("PID selection bits: n = %d, \n",fNPIDBits) ;
807 printf("\tPID bit = ");
808 for(Int_t i = 0; i < fNPIDBits; i++) printf("%d ",fPIDBits[i]);
809 printf("\n");
810
db2bf6fd 811 if(fMultiCutAna){
812 printf("pT cuts: n = %d, \n",fNPtCuts) ;
813 printf("\tpT > ");
814 for(Int_t i = 0; i < fNPtCuts; i++) printf("%2.2f ",fPtCuts[i]);
815 printf("GeV/c\n");
816
817 printf("N cell in cluster cuts: n = %d, \n",fNCellNCuts) ;
818 printf("\tnCell > ");
819 for(Int_t i = 0; i < fNCellNCuts; i++) printf("%d ",fCellNCuts[i]);
820 printf("\n");
821
db2bf6fd 822 }
477d6cee 823 printf("------------------------------------------------------\n") ;
1c5acb87 824}
825
5ae09196 826//_____________________________________________________________
827void AliAnaPi0::FillAcceptanceHistograms(){
828 //Fill acceptance histograms if MC data is available
c8fe2783 829
6175da48 830 if(GetReader()->ReadStack()){
5ae09196 831 AliStack * stack = GetMCStack();
6175da48 832 if(stack){
5ae09196 833 for(Int_t i=0 ; i<stack->GetNprimary(); i++){
834 TParticle * prim = stack->Particle(i) ;
6175da48 835 Int_t pdg = prim->GetPdgCode();
836 if( pdg == 111 || pdg == 221){
5ae09196 837 Double_t pi0Pt = prim->Pt() ;
838 //printf("pi0, pt %2.2f\n",pi0Pt);
839 if(prim->Energy() == TMath::Abs(prim->Pz())) continue ; //Protection against floating point exception
840 Double_t pi0Y = 0.5*TMath::Log((prim->Energy()-prim->Pz())/(prim->Energy()+prim->Pz())) ;
841 Double_t phi = TMath::RadToDeg()*prim->Phi() ;
6175da48 842 if(pdg == 111){
843 if(TMath::Abs(pi0Y) < 0.5){
844 fhPrimPi0Pt->Fill(pi0Pt) ;
845 }
846 fhPrimPi0Y ->Fill(pi0Y) ;
847 fhPrimPi0Phi->Fill(phi) ;
848 }
849 else if(pdg == 221){
850 if(TMath::Abs(pi0Y) < 0.5){
851 fhPrimEtaPt->Fill(pi0Pt) ;
852 }
853 fhPrimEtaY ->Fill(pi0Y) ;
854 fhPrimEtaPhi->Fill(phi) ;
5ae09196 855 }
5ae09196 856 //Check if both photons hit Calorimeter
6175da48 857 if(prim->GetNDaughters()!=2) return; //Only interested in 2 gamma decay
5ae09196 858 Int_t iphot1=prim->GetFirstDaughter() ;
859 Int_t iphot2=prim->GetLastDaughter() ;
860 if(iphot1>-1 && iphot1<stack->GetNtrack() && iphot2>-1 && iphot2<stack->GetNtrack()){
861 TParticle * phot1 = stack->Particle(iphot1) ;
862 TParticle * phot2 = stack->Particle(iphot2) ;
863 if(phot1 && phot2 && phot1->GetPdgCode()==22 && phot2->GetPdgCode()==22){
864 //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",
865 // phot1->Pt(), phot1->Phi()*180./3.1415, phot1->Eta(), phot2->Pt(), phot2->Phi()*180./3.1415, phot2->Eta());
866
867 TLorentzVector lv1, lv2;
868 phot1->Momentum(lv1);
869 phot2->Momentum(lv2);
870
871 Bool_t inacceptance = kFALSE;
872 if(fCalorimeter == "PHOS"){
873 if(GetPHOSGeometry() && GetCaloUtils()->IsPHOSGeoMatrixSet()){
874 Int_t mod ;
875 Double_t x,z ;
876 if(GetPHOSGeometry()->ImpactOnEmc(phot1,mod,z,x) && GetPHOSGeometry()->ImpactOnEmc(phot2,mod,z,x))
877 inacceptance = kTRUE;
878 if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
879 }
880 else{
881
882 if(GetFiducialCut()->IsInFiducialCut(lv1,fCalorimeter) && GetFiducialCut()->IsInFiducialCut(lv2,fCalorimeter))
883 inacceptance = kTRUE ;
884 if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
885 }
886
887 }
888 else if(fCalorimeter == "EMCAL" && GetCaloUtils()->IsEMCALGeoMatrixSet()){
889 if(GetEMCALGeometry()){
890 if(GetEMCALGeometry()->Impact(phot1) && GetEMCALGeometry()->Impact(phot2))
891 inacceptance = kTRUE;
892 if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
893 }
894 else{
895 if(GetFiducialCut()->IsInFiducialCut(lv1,fCalorimeter) && GetFiducialCut()->IsInFiducialCut(lv2,fCalorimeter))
896 inacceptance = kTRUE ;
897 if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
898 }
899 }
900
901 if(inacceptance){
6175da48 902 if(pdg==111){
903 fhPrimPi0AccPt->Fill(pi0Pt) ;
904 fhPrimPi0AccPhi->Fill(phi) ;
905 fhPrimPi0AccY->Fill(pi0Y) ;
906 Double_t angle = lv1.Angle(lv2.Vect());
907 fhPrimPi0OpeningAngle ->Fill(pi0Pt,angle);
908 fhPrimPi0CosOpeningAngle->Fill(pi0Pt,TMath::Cos(angle));
909 }
910 else if(pdg==221){
911 fhPrimEtaAccPt->Fill(pi0Pt) ;
912 fhPrimEtaAccPhi->Fill(phi) ;
913 fhPrimEtaAccY->Fill(pi0Y) ;
914 }
5ae09196 915 }//Accepted
916 }// 2 photons
917 }//Check daughters exist
918 }// Primary pi0
919 }//loop on primaries
920 }//stack exists and data is MC
921 }//read stack
922 else if(GetReader()->ReadAODMCParticles()){
6175da48 923
924 TClonesArray * mcparticles = GetReader()->GetAODMCParticles(0);
925 if(mcparticles){
926 Int_t nprim = mcparticles->GetEntriesFast();
927 for(Int_t i=0 ; i < nprim; i++){
928 AliAODMCParticle * prim = (AliAODMCParticle *) mcparticles->At(i);
929 Int_t pdg = prim->GetPdgCode();
930 if( pdg == 111 || pdg == 221){
931 Double_t pi0Pt = prim->Pt() ;
932 //printf("pi0, pt %2.2f\n",pi0Pt);
933 if(prim->E() == TMath::Abs(prim->Pz())) continue ; //Protection against floating point exception
934 Double_t pi0Y = 0.5*TMath::Log((prim->E()-prim->Pz())/(prim->E()+prim->Pz())) ;
935 Double_t phi = TMath::RadToDeg()*prim->Phi() ;
936 if(pdg == 111){
937 if(TMath::Abs(pi0Y) < 0.5){
938 fhPrimPi0Pt->Fill(pi0Pt) ;
939 }
940 fhPrimPi0Y ->Fill(pi0Y) ;
941 fhPrimPi0Phi->Fill(phi) ;
942 }
943 else if(pdg == 221){
944 if(TMath::Abs(pi0Y) < 0.5){
945 fhPrimEtaPt->Fill(pi0Pt) ;
946 }
947 fhPrimEtaY ->Fill(pi0Y) ;
948 fhPrimEtaPhi->Fill(phi) ;
949 }
950 //Check if both photons hit Calorimeter
951 if(prim->GetNDaughters()!=2) return; //Only interested in 2 gamma decay
952 Int_t iphot1=prim->GetDaughter(0) ;
953 Int_t iphot2=prim->GetDaughter(1) ;
954 if(iphot1>-1 && iphot1<nprim && iphot2>-1 && iphot2<nprim){
955 AliAODMCParticle * phot1 = (AliAODMCParticle *) mcparticles->At(iphot1);
956 AliAODMCParticle * phot2 = (AliAODMCParticle *) mcparticles->At(iphot2);
957 if(phot1 && phot2 && phot1->GetPdgCode()==22 && phot2->GetPdgCode()==22){
958 //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",
959 // phot1->Pt(), phot1->Phi()*180./3.1415, phot1->Eta(), phot2->Pt(), phot2->Phi()*180./3.1415, phot2->Eta());
960
961 TLorentzVector lv1, lv2;
962 lv1.SetPxPyPzE(phot1->Px(),phot1->Py(),phot1->Pz(),phot1->E());
963 lv2.SetPxPyPzE(phot2->Px(),phot2->Py(),phot2->Pz(),phot2->E());
964
965 Bool_t inacceptance = kFALSE;
966 if(fCalorimeter == "PHOS"){
967 if(GetPHOSGeometry() && GetCaloUtils()->IsPHOSGeoMatrixSet()){
968 Int_t mod ;
969 Double_t x,z ;
970 Double_t vtx []={phot1->Xv(),phot1->Yv(),phot1->Zv()};
971 Double_t vtx2[]={phot2->Xv(),phot2->Yv(),phot2->Zv()};
972 if(GetPHOSGeometry()->ImpactOnEmc(vtx, phot1->Theta(),phot1->Phi(),mod,z,x) &&
973 GetPHOSGeometry()->ImpactOnEmc(vtx2,phot2->Theta(),phot2->Phi(),mod,z,x))
974 inacceptance = kTRUE;
975 if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
976 }
977 else{
978
979 if(GetFiducialCut()->IsInFiducialCut(lv1,fCalorimeter) && GetFiducialCut()->IsInFiducialCut(lv2,fCalorimeter))
980 inacceptance = kTRUE ;
981 if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
982 }
983
984 }
985 else if(fCalorimeter == "EMCAL" && GetCaloUtils()->IsEMCALGeoMatrixSet()){
986 if(GetEMCALGeometry()){
987 TVector3 vtx(phot1->Xv(),phot1->Yv(),phot1->Zv());
988 TVector3 vimpact(0,0,0);
989 Int_t absID1=0;
990
991 GetEMCALGeometry()->ImpactOnEmcal(vtx,phot1->Theta(),phot1->Phi(),absID1,vimpact);
992 TVector3 vtx2(phot2->Xv(),phot2->Yv(),phot2->Zv());
993 TVector3 vimpact2(0,0,0);
994 Int_t absID2=0;
995 GetEMCALGeometry()->ImpactOnEmcal(vtx2,phot2->Theta(),phot2->Phi(),absID2,vimpact2);
996// if(TMath::Abs(phot1->Eta()) < 0.7 && phot1->Phi() > 80*TMath::DegToRad() && phot1->Phi() < 120*TMath::DegToRad() )
997// printf("Phot1 accepted? %d\n",absID1);
998// if(TMath::Abs(phot2->Eta()) < 0.7 && phot2->Phi() > 80*TMath::DegToRad() && phot2->Phi() < 120*TMath::DegToRad() )
999// printf("Phot2 accepted? %d\n",absID2);
1000
1001 if( absID1 >= 0 && absID2 >= 0)
1002 inacceptance = kTRUE;
1003 if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
1004 }
1005 else{
1006 if(GetFiducialCut()->IsInFiducialCut(lv1,fCalorimeter) && GetFiducialCut()->IsInFiducialCut(lv2,fCalorimeter))
1007 inacceptance = kTRUE ;
1008 if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
1009 }
1010 }
1011
1012 if(inacceptance){
1013 if(pdg==111){
1014 fhPrimPi0AccPt->Fill(pi0Pt) ;
1015 fhPrimPi0AccPhi->Fill(phi) ;
1016 fhPrimPi0AccY->Fill(pi0Y) ;
1017 Double_t angle = lv1.Angle(lv2.Vect());
1018 fhPrimPi0OpeningAngle ->Fill(pi0Pt,angle);
1019 fhPrimPi0CosOpeningAngle->Fill(pi0Pt,TMath::Cos(angle));
1020 }
1021 else if(pdg==221){
1022 fhPrimEtaAccPt->Fill(pi0Pt) ;
1023 fhPrimEtaAccPhi->Fill(phi) ;
1024 fhPrimEtaAccY->Fill(pi0Y) ;
1025 }
1026 }//Accepted
1027 }// 2 photons
1028 }//Check daughters exist
1029 }// Primary pi0
1030 }//loop on primaries
1031 }//stack exists and data is MC
1032
1033
1034 } // read AOD MC
5ae09196 1035}
1036
6175da48 1037//_____________________________________________________________
1038void AliAnaPi0::FillMCVersusRecDataHistograms(const Int_t index1, const Int_t index2,
1039 const Float_t pt1, const Float_t pt2,
1040 const Int_t ncell1, const Int_t ncell2,
1041 const Double_t mass, const Double_t pt, const Double_t asym,
1042 const Double_t deta, const Double_t dphi){
1043 //Do some MC checks on the origin of the pair, is there any common ancestor and if there is one, who?
1044 //Adjusted for Pythia, need to see what to do for other generators.
1045 //Array of histograms ordered as follows: 0-Photon, 1-electron, 2-pi0, 3-eta, 4-a-proton, 5-a-neutron, 6-stable particles,
1046 // 7-other decays, 8-string, 9-final parton, 10-initial parton, intermediate, 11-colliding proton, 12-unrelated
1047
1048 Int_t ancPDG = 0;
1049 Int_t ancStatus = 0;
1050 TLorentzVector ancMomentum;
1051 Int_t ancLabel = GetMCAnalysisUtils()->CheckCommonAncestor(index1, index2,
1052 GetReader(), ancPDG, ancStatus,ancMomentum);
1053
1054 if(GetDebug() > 1) printf("AliAnaPi0::FillMCVersusRecDataHistograms() - Common ancestor label %d, pdg %d, name %s, status %d; \n",
1055 ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus);
1056
1057 if(ancLabel > -1){
1058 if(ancPDG==22){//gamma
1059 fhMCOrgMass[0]->Fill(pt,mass);
1060 fhMCOrgAsym[0]->Fill(pt,asym);
1061 fhMCOrgDeltaEta[0]->Fill(pt,deta);
1062 fhMCOrgDeltaPhi[0]->Fill(pt,dphi);
1063 }
1064 else if(TMath::Abs(ancPDG)==11){//e
1065 fhMCOrgMass[1]->Fill(pt,mass);
1066 fhMCOrgAsym[1]->Fill(pt,asym);
1067 fhMCOrgDeltaEta[1]->Fill(pt,deta);
1068 fhMCOrgDeltaPhi[1]->Fill(pt,dphi);
1069 }
1070 else if(ancPDG==111){//Pi0
1071 fhMCOrgMass[2]->Fill(pt,mass);
1072 fhMCOrgAsym[2]->Fill(pt,asym);
1073 fhMCOrgDeltaEta[2]->Fill(pt,deta);
1074 fhMCOrgDeltaPhi[2]->Fill(pt,dphi);
1075 if(fMultiCutAnaSim){
1076 for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
1077 for(Int_t icell=0; icell<fNCellNCuts; icell++){
1078 for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
1079 Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
1080 if(pt1 > fPtCuts[ipt] && pt2 > fPtCuts[ipt] &&
1081 asym < fAsymCuts[iasym] &&
1082 ncell1 >= fCellNCuts[icell] && ncell2 >= fCellNCuts[icell]){
1083 fhMCPi0MassPtRec [index]->Fill(pt,mass);
1084 fhMCPi0MassPtTrue[index]->Fill(ancMomentum.Pt(),mass);
1085 if(mass < 0.17 && mass > 0.1) fhMCPi0PtTruePtRec[index]->Fill(ancMomentum.Pt(),pt);
1086 }//pass the different cuts
1087 }// pid bit cut loop
1088 }// icell loop
1089 }// pt cut loop
1090 }//Multi cut ana sim
1091 else {
1092 fhMCPi0MassPtTrue[0]->Fill(ancMomentum.Pt(),mass);
1093 if(mass < 0.17 && mass > 0.1) fhMCPi0PtTruePtRec[0]->Fill(ancMomentum.Pt(),pt);
1094 }
1095 }
1096 else if(ancPDG==221){//Eta
1097 fhMCOrgMass[3]->Fill(pt,mass);
1098 fhMCOrgAsym[3]->Fill(pt,asym);
1099 fhMCOrgDeltaEta[3]->Fill(pt,deta);
1100 fhMCOrgDeltaPhi[3]->Fill(pt,dphi);
1101 if(fMultiCutAnaSim){
1102 for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
1103 for(Int_t icell=0; icell<fNCellNCuts; icell++){
1104 for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
1105 Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
1106 if(pt1 > fPtCuts[ipt] && pt2 > fPtCuts[ipt] &&
1107 asym < fAsymCuts[iasym] &&
1108 ncell1 >= fCellNCuts[icell] && ncell2 >= fCellNCuts[icell]){
1109 fhMCEtaMassPtRec [index]->Fill(pt,mass);
1110 fhMCEtaMassPtTrue[index]->Fill(ancMomentum.Pt(),mass);
1111 if(mass < 0.17 && mass > 0.1) fhMCEtaPtTruePtRec[index]->Fill(ancMomentum.Pt(),pt);
1112 }//pass the different cuts
1113 }// pid bit cut loop
1114 }// icell loop
1115 }// pt cut loop
1116 } //Multi cut ana sim
1117 else {
1118 fhMCEtaMassPtTrue[0]->Fill(ancMomentum.Pt(),mass);
1119 if(mass < 0.17 && mass > 0.1) fhMCEtaPtTruePtRec[0]->Fill(ancMomentum.Pt(),pt);
1120 }
1121 }
1122 else if(ancPDG==-2212){//AProton
1123 fhMCOrgMass[4]->Fill(pt,mass);
1124 fhMCOrgAsym[4]->Fill(pt,asym);
1125 fhMCOrgDeltaEta[4]->Fill(pt,deta);
1126 fhMCOrgDeltaPhi[4]->Fill(pt,dphi);
1127 }
1128 else if(ancPDG==-2112){//ANeutron
1129 fhMCOrgMass[5]->Fill(pt,mass);
1130 fhMCOrgAsym[5]->Fill(pt,asym);
1131 fhMCOrgDeltaEta[5]->Fill(pt,deta);
1132 fhMCOrgDeltaPhi[5]->Fill(pt,dphi);
1133 }
1134 else if(TMath::Abs(ancPDG)==13){//muons
1135 fhMCOrgMass[6]->Fill(pt,mass);
1136 fhMCOrgAsym[6]->Fill(pt,asym);
1137 fhMCOrgDeltaEta[6]->Fill(pt,deta);
1138 fhMCOrgDeltaPhi[6]->Fill(pt,dphi);
1139 }
1140 else if (TMath::Abs(ancPDG) > 100 && ancLabel > 7) {
1141 if(ancStatus==1){//Stable particles, converted? not decayed resonances
1142 fhMCOrgMass[6]->Fill(pt,mass);
1143 fhMCOrgAsym[6]->Fill(pt,asym);
1144 fhMCOrgDeltaEta[6]->Fill(pt,deta);
1145 fhMCOrgDeltaPhi[6]->Fill(pt,dphi);
1146 }
1147 else{//resonances and other decays, more hadron conversions?
1148 fhMCOrgMass[7]->Fill(pt,mass);
1149 fhMCOrgAsym[7]->Fill(pt,asym);
1150 fhMCOrgDeltaEta[7]->Fill(pt,deta);
1151 fhMCOrgDeltaPhi[7]->Fill(pt,dphi);
1152 }
1153 }
1154 else {//Partons, colliding protons, strings, intermediate corrections
1155 if(ancStatus==11 || ancStatus==12){//String fragmentation
1156 fhMCOrgMass[8]->Fill(pt,mass);
1157 fhMCOrgAsym[8]->Fill(pt,asym);
1158 fhMCOrgDeltaEta[8]->Fill(pt,deta);
1159 fhMCOrgDeltaPhi[8]->Fill(pt,dphi);
1160 }
1161 else if (ancStatus==21){
1162 if(ancLabel < 2) {//Colliding protons
1163 fhMCOrgMass[11]->Fill(pt,mass);
1164 fhMCOrgAsym[11]->Fill(pt,asym);
1165 fhMCOrgDeltaEta[11]->Fill(pt,deta);
1166 fhMCOrgDeltaPhi[11]->Fill(pt,dphi);
1167 }//colliding protons
1168 else if(ancLabel < 6){//partonic initial states interactions
1169 fhMCOrgMass[9]->Fill(pt,mass);
1170 fhMCOrgAsym[9]->Fill(pt,asym);
1171 fhMCOrgDeltaEta[9]->Fill(pt,deta);
1172 fhMCOrgDeltaPhi[9]->Fill(pt,dphi);
1173 }
1174 else if(ancLabel < 8){//Final state partons radiations?
1175 fhMCOrgMass[10]->Fill(pt,mass);
1176 fhMCOrgAsym[10]->Fill(pt,asym);
1177 fhMCOrgDeltaEta[10]->Fill(pt,deta);
1178 fhMCOrgDeltaPhi[10]->Fill(pt,dphi);
1179 }
1180 else {
1181 printf("AliAnaPi0::FillMCVersusRecDataHistograms() - Check ** Common ancestor label %d, pdg %d, name %s, status %d; \n",
1182 ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus);
1183 }
1184 }//status 21
1185 else {
1186 printf("AliAnaPi0::FillMCVersusRecDataHistograms() - Check *** Common ancestor label %d, pdg %d, name %s, status %d; \n",
1187 ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus);
1188 }
1189 }////Partons, colliding protons, strings, intermediate corrections
1190 }//ancLabel > -1
1191 else { //ancLabel <= -1
1192 //printf("Not related at all label = %d\n",ancLabel);
1193 fhMCOrgMass[12]->Fill(pt,mass);
1194 fhMCOrgAsym[12]->Fill(pt,asym);
1195 fhMCOrgDeltaEta[12]->Fill(pt,deta);
1196 fhMCOrgDeltaPhi[12]->Fill(pt,dphi);
1197 }
1198}
1199
1c5acb87 1200//____________________________________________________________________________________________________________________________________________________
6639984f 1201void AliAnaPi0::MakeAnalysisFillHistograms()
1c5acb87 1202{
477d6cee 1203 //Process one event and extract photons from AOD branch
1204 // filled with AliAnaPhoton and fill histos with invariant mass
1205
6175da48 1206 //In case of simulated data, fill acceptance histograms
1207 if(IsDataMC())FillAcceptanceHistograms();
1208
1209 //Init some variables
1210//Int_t iRun = (GetReader()->GetInputEvent())->GetRunNumber() ;
1211 Int_t nPhot = GetInputAODBranch()->GetEntriesFast() ;
1212 Int_t nClus = 0;
1213 Int_t nCell = 0;
1214 Float_t eClusTot = 0;
1215 Float_t eCellTot = 0;
477d6cee 1216
6175da48 1217 // Count the number of clusters and cells, in case multiplicity bins dependent on such numbers
1218 // are requested
1219 if(fCalorimeter=="EMCAL"){
1220 nClus = GetAODEMCAL() ->GetEntriesFast();
1221 nCell = GetEMCALCells()->GetNumberOfCells();
1222 for(Int_t icl=0; icl < nClus; icl++) eClusTot += ((AliVCluster*) GetAODEMCAL()->At(icl))->E();
1223 for(Int_t jce=0; jce < nCell; jce++) eCellTot += GetEMCALCells()->GetAmplitude(jce);
1224
1225 }
1226 else {
1227 nClus = GetAODPHOS() ->GetEntriesFast();
1228 nCell = GetPHOSCells()->GetNumberOfCells();
1229 for(Int_t icl=0; icl < nClus; icl++) eClusTot += ((AliVCluster*)GetAODPHOS()->At(icl))->E();
1230 for(Int_t jce=0; jce < nCell; jce++) eCellTot += GetPHOSCells()->GetAmplitude(jce);
1231 }
1232
1233 //Fill the average number of cells or clusters per SM
1234 eClusTot /=fNModules;
1235 eCellTot /=fNModules;
1236 fhAverTotECluster->Fill(eClusTot);
1237 fhAverTotECell ->Fill(eCellTot);
1238
c8fe2783 1239 if(GetDebug() > 1)
1240 printf("AliAnaPi0::MakeAnalysisFillHistograms() - Photon entries %d\n", nPhot);
6175da48 1241
1242 //If less than photon 2 entries in the list, skip this event
1243 if(nPhot < 2 ) return ;
7e7694bb 1244
6175da48 1245 if(GetDebug() > 1)
1246 printf("AliAnaPi0::MakeAnalysisFillHistograms() - # Clusters %d, sum cluster E per SM %f,# Cells %d, sum cell E per SM %f\n", nClus,eClusTot,nCell,eCellTot);
1247
1248 //Init variables
1249 Int_t module1 = -1;
1250 Int_t module2 = -1;
1251 Double_t vert[] = {0.0, 0.0, 0.0} ; //vertex
1252 Int_t evtIndex1 = 0 ;
1253 Int_t currentEvtIndex = -1;
1254 Int_t curCentrBin = 0 ;
1255 Int_t curRPBin = 0 ;
1256 Int_t curZvertBin = 0 ;
1257
1258 //---------------------------------
1259 //First loop on photons/clusters
1260 //---------------------------------
477d6cee 1261 for(Int_t i1=0; i1<nPhot-1; i1++){
1262 AliAODPWG4Particle * p1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i1)) ;
6175da48 1263 //printf("AliAnaPi0::MakeAnalysisFillHistograms() : cluster1 id %d\n",p1->GetCaloLabel(0));
1264
7e7694bb 1265 // get the event index in the mixed buffer where the photon comes from
1266 // in case of mixing with analysis frame, not own mixing
c8fe2783 1267 evtIndex1 = GetEventIndex(p1, vert) ;
5025c139 1268 //printf("charge = %d\n", track->Charge());
c8fe2783 1269 if ( evtIndex1 == -1 )
1270 return ;
1271 if ( evtIndex1 == -2 )
1272 continue ;
2244659d 1273 if(TMath::Abs(vert[2]) > GetZvertexCut()) continue ; //vertex cut
6175da48 1274
1275 //----------------------------------------------------------------------------
1276 // Get the multiplicity bin. Different cases: centrality (PbPb),
1277 // average cluster multiplicity, average cell multiplicity, track multiplicity
1278 // default is centrality bins
1279 //----------------------------------------------------------------------------
c8fe2783 1280 if (evtIndex1 != currentEvtIndex) {
6175da48 1281 if(fUseTrackMultBins){ // Track multiplicity bins
1282 //printf("track mult %d\n",GetTrackMultiplicity());
1283 curCentrBin = (GetTrackMultiplicity()-1)/5;
1284 if(curCentrBin > fNCentrBin-1) curCentrBin=fNCentrBin-1;
1285 //printf("track mult bin %d\n",curCentrBin);
1286 }
1287 else if(fUsePhotonMultBins){ // Photon multiplicity bins
1288 //printf("photon mult %d cluster mult %d\n",nPhot, nClus);
1289 curCentrBin = (nClus-1)/3;
1290 if(curCentrBin > fNCentrBin-1) curCentrBin=fNCentrBin-1;
1291 //printf("photon mult bin %d\n",curCentrBin);
1292 }
1293 else if(fUseAverClusterEBins){ // Cluster multiplicity bins
1294 //Bins for pp, if needed can be done in a more general way
1295 if (eClusTot < 0.5 )curCentrBin = 0;
1296 else if(eClusTot < 1.0) curCentrBin = 1;
1297 else if(eClusTot < 1.5) curCentrBin = 2;
1298 else if(eClusTot < 2.0) curCentrBin = 3;
1299 else if(eClusTot < 3.0) curCentrBin = 4;
1300 else if(eClusTot < 4.0) curCentrBin = 5;
1301 else if(eClusTot < 5.0) curCentrBin = 6;
1302 else if(eClusTot < 7.5) curCentrBin = 7;
1303 else if(eClusTot < 10.) curCentrBin = 8;
1304 else if(eClusTot < 15.) curCentrBin = 9;
1305 else if(eClusTot < 20.) curCentrBin = 10;
1306 if(curCentrBin > fNCentrBin-1) curCentrBin=fNCentrBin-1;
1307 //printf("cluster E average %f, bin %d \n",eClusTot,curCentrBin);
1308 }
1309 else if(fUseAverCellEBins){ // Cell multiplicity bins
1310 //Bins for pp, if needed can be done in a more general way
1311 if (eCellTot < 0.5) curCentrBin = 0;
1312 else if(eCellTot < 1.0) curCentrBin = 1;
1313 else if(eCellTot < 1.5) curCentrBin = 2;
1314 else if(eCellTot < 2.0) curCentrBin = 3;
1315 else if(eCellTot < 3.0) curCentrBin = 4;
1316 else if(eCellTot < 4.0) curCentrBin = 5;
1317 else if(eCellTot < 5.0) curCentrBin = 6;
1318 else if(eCellTot < 7.5) curCentrBin = 7;
1319 else if(eCellTot < 10.) curCentrBin = 8;
1320 else if(eCellTot < 15.) curCentrBin = 9;
1321 else if(eCellTot < 20.) curCentrBin = 10;
1322 if(curCentrBin > fNCentrBin-1) curCentrBin=fNCentrBin-1;
1323 //printf("cell E average %f, bin %d \n",eCellTot,curCentrBin);
1324 }
1325 else { //Event centrality
1326 curCentrBin = GetEventCentrality();
1327 }
1328
1329 //Get vertex z bin
c8fe2783 1330 curRPBin = 0 ;
5025c139 1331 curZvertBin = (Int_t)(0.5*GetNZvertBin()*(vert[2]+GetZvertexCut())/GetZvertexCut()) ;
6175da48 1332
1333 //Fill event bin info
c8fe2783 1334 fhEvents->Fill(curCentrBin+0.5,curZvertBin+0.5,curRPBin+0.5) ;
1335 currentEvtIndex = evtIndex1 ;
ca468d44 1336 if(GetDebug() > 1)
6175da48 1337 printf("AliAnaPi0::MakeAnalysisFillHistograms() - Centrality %d, Vertex Bin %d, RP bin %d \n",curCentrBin,curRPBin,curZvertBin);
c8fe2783 1338 }
7e7694bb 1339
f8006433 1340 //printf("AliAnaPi0::MakeAnalysisFillHistograms(): Photon 1 Evt %d Vertex : %f,%f,%f\n",evtIndex1, GetVertex(evtIndex1)[0] ,GetVertex(evtIndex1)[1],GetVertex(evtIndex1)[2]);
af7b3903 1341
6175da48 1342 //Get the momentum of this cluster
477d6cee 1343 TLorentzVector photon1(p1->Px(),p1->Py(),p1->Pz(),p1->E());
6175da48 1344
1345 //Get (Super)Module number of this cluster
59b6bd99 1346 module1 = GetModuleNumber(p1);
6175da48 1347
1348 //---------------------------------
1349 //Second loop on photons/clusters
1350 //---------------------------------
477d6cee 1351 for(Int_t i2=i1+1; i2<nPhot; i2++){
1352 AliAODPWG4Particle * p2 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i2)) ;
6175da48 1353
1354 //In case of mixing frame, check we are not in the same event as the first cluster
c8fe2783 1355 Int_t evtIndex2 = GetEventIndex(p2, vert) ;
1356 if ( evtIndex2 == -1 )
1357 return ;
1358 if ( evtIndex2 == -2 )
1359 continue ;
1360 if (GetMixedEvent() && (evtIndex1 == evtIndex2))
7e7694bb 1361 continue ;
6175da48 1362
f8006433 1363 //printf("AliAnaPi0::MakeAnalysisFillHistograms(): Photon 2 Evt %d Vertex : %f,%f,%f\n",evtIndex2, GetVertex(evtIndex2)[0] ,GetVertex(evtIndex2)[1],GetVertex(evtIndex2)[2]);
6175da48 1364
1365 //Get the momentum of this cluster
477d6cee 1366 TLorentzVector photon2(p2->Px(),p2->Py(),p2->Pz(),p2->E());
59b6bd99 1367 //Get module number
6175da48 1368 module2 = GetModuleNumber(p2);
1369
1370 //---------------------------------
1371 // Get pair kinematics
1372 //---------------------------------
1373 Double_t m = (photon1 + photon2).M() ;
1374 Double_t pt = (photon1 + photon2).Pt();
1375 Double_t deta = photon1.Eta() - photon2.Eta();
1376 Double_t dphi = photon1.Phi() - photon2.Phi();
1377 Double_t a = TMath::Abs(p1->E()-p2->E())/(p1->E()+p2->E()) ;
1378
477d6cee 1379 if(GetDebug() > 2)
6175da48 1380 printf(" E: photon1 %f, photon2 %f; Pair: pT %f, mass %f, a %f\n", p1->E(), p2->E(), (photon1 + photon2).E(),m,a);
1381
1382 //--------------------------------
1383 // Opening angle selection
1384 //--------------------------------
50f39b97 1385 //Check if opening angle is too large or too small compared to what is expected
1386 Double_t angle = photon1.Angle(photon2.Vect());
6175da48 1387 if(fUseAngleEDepCut && !GetNeutralMesonSelection()->IsAngleInWindow((photon1+photon2).E(),angle+0.05)) {
1388 if(GetDebug() > 2)
1389 printf("AliAnaPi0::MakeAnalysisFillHistograms() -Real pair angle %f not in E %f window\n",angle, (photon1+photon2).E());
c8fe2783 1390 continue;
6175da48 1391 }
af7b3903 1392
6175da48 1393 if(fUseAngleCut && (angle < fAngleCut || angle > fAngleMaxCut)) {
1394 if(GetDebug() > 2)
1395 printf("AliAnaPi0::MakeAnalysisFillHistograms() - Real pair cut %f < angle %f < cut %f\n",fAngleCut, angle, fAngleMaxCut);
1396 continue;
1397 }
1398
1399 //-------------------------------------------------------------------------------------------------
af7b3903 1400 //Fill module dependent histograms, put a cut on assymmetry on the first available cut in the array
6175da48 1401 //-------------------------------------------------------------------------------------------------
af7b3903 1402 if(a < fAsymCuts[0]){
1403 if(module1==module2 && module1 >=0 && module1<fNModules)
1404 fhReMod[module1]->Fill(pt,m) ;
1405 else
6175da48 1406 fhReDiffMod[fNModules+2]->Fill(pt,m) ;
af7b3903 1407
1408 if(fCalorimeter=="EMCAL"){
1409 if((module1==0 && module2==2) || (module1==2 && module2==0)) fhReDiffMod[0]->Fill(pt,m) ;
1410 if((module1==1 && module2==3) || (module1==3 && module2==1)) fhReDiffMod[1]->Fill(pt,m) ;
1411 if((module1==0 && module2==1) || (module1==1 && module2==0)) fhReDiffMod[2]->Fill(pt,m) ;
6175da48 1412 if((module1==2 && module2==3) || (module1==3 && module2==2)) fhReDiffMod[3]->Fill(pt,m) ;
1413 if((module1==0 && module2==3) || (module1==3 && module2==0)) fhReDiffMod[4]->Fill(pt,m) ;
1414 if((module1==2 && module2==1) || (module1==1 && module2==2)) fhReDiffMod[5]->Fill(pt,m) ;
af7b3903 1415 }
1416 else {
1417 if((module1==0 && module2==1) || (module1==1 && module2==0)) fhReDiffMod[0]->Fill(pt,m) ;
1418 if((module1==0 && module2==2) || (module1==2 && module2==0)) fhReDiffMod[1]->Fill(pt,m) ;
1419 if((module1==1 && module2==2) || (module1==2 && module2==1)) fhReDiffMod[2]->Fill(pt,m) ;
1420 }
821c8090 1421 }
7e7694bb 1422
af7b3903 1423 //In case we want only pairs in same (super) module, check their origin.
1424 Bool_t ok = kTRUE;
1425 if(fSameSM && module1!=module2) ok=kFALSE;
1426 if(ok){
6175da48 1427
1428 //Check if one of the clusters comes from a conversion
1429 if (p1->IsTagged() && p2->IsTagged()) fhReConv2->Fill(pt,m);
1430 else if(p1->IsTagged() || p2->IsTagged()) fhReConv ->Fill(pt,m);
1431
af7b3903 1432 //Fill histograms for different bad channel distance, centrality, assymmetry cut and pid bit
5ae09196 1433 for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
af7b3903 1434 if((p1->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton)) && (p2->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton))){
1435 for(Int_t iasym=0; iasym < fNAsymCuts; iasym++){
1436 if(a < fAsymCuts[iasym]){
1437 Int_t index = ((curCentrBin*fNPIDBits)+ipid)*fNAsymCuts + iasym;
6175da48 1438 //printf("index %d :(cen %d * nPID %d + ipid %d)*nasym %d + iasym %d\n",index,curCentrBin,fNPIDBits,ipid,fNAsymCuts,iasym);
af7b3903 1439 fhRe1 [index]->Fill(pt,m);
398c93cc 1440 if(fMakeInvPtPlots)fhReInvPt1[index]->Fill(pt,m,1./pt) ;
6175da48 1441 if(fFillBadDistHisto){
1442 if(p1->DistToBad()>0 && p2->DistToBad()>0){
1443 fhRe2 [index]->Fill(pt,m) ;
1444 if(fMakeInvPtPlots)fhReInvPt2[index]->Fill(pt,m,1./pt) ;
1445 if(p1->DistToBad()>1 && p2->DistToBad()>1){
1446 fhRe3 [index]->Fill(pt,m) ;
1447 if(fMakeInvPtPlots)fhReInvPt3[index]->Fill(pt,m,1./pt) ;
1448 }// bad 3
1449 }// bad2
1450 }// Fill bad dist histos
1451 }//assymetry cut
1452 }// asymmetry cut loop
af7b3903 1453 }// bad 1
1454 }// pid bit loop
5ae09196 1455
af7b3903 1456 //Fill histograms with opening angle
1457 fhRealOpeningAngle ->Fill(pt,angle);
1458 fhRealCosOpeningAngle->Fill(pt,TMath::Cos(angle));
1459
1460 //Fill histograms with pair assymmetry
1461 fhRePtAsym->Fill(pt,a);
6175da48 1462 if(m > 0.10 && m < 0.17) fhRePtAsymPi0->Fill(pt,a);
af7b3903 1463 if(m > 0.45 && m < 0.65) fhRePtAsymEta->Fill(pt,a);
1464
6175da48 1465 //-------------------------------------------------------
1466 //Get the number of cells needed for multi cut analysis.
1467 //-------------------------------------------------------
1468 Int_t ncell1 = 0;
1469 Int_t ncell2 = 0;
1470 if(fMultiCutAna || (IsDataMC() && fMultiCutAnaSim)){
1471
af7b3903 1472 AliVEvent * event = GetReader()->GetInputEvent();
1473 if(event){
1474 for(Int_t iclus = 0; iclus < event->GetNumberOfCaloClusters(); iclus++){
1475 AliVCluster *cluster = event->GetCaloCluster(iclus);
5ae09196 1476
af7b3903 1477 Bool_t is = kFALSE;
1478 if (fCalorimeter == "EMCAL" && GetReader()->IsEMCALCluster(cluster)) is = kTRUE;
1479 else if(fCalorimeter == "PHOS" && GetReader()->IsPHOSCluster (cluster)) is = kTRUE;
5ae09196 1480
af7b3903 1481 if(is){
1482 if (p1->GetCaloLabel(0) == cluster->GetID()) ncell1 = cluster->GetNCells();
1483 else if (p2->GetCaloLabel(0) == cluster->GetID()) ncell2 = cluster->GetNCells();
1484 } // PHOS or EMCAL cluster as requested in analysis
1485
1486 if(ncell2 > 0 && ncell1 > 0) break; // No need to continue the iteration
1487
1488 }
1489 //printf("e 1: %2.2f, e 2: %2.2f, ncells: n1 %d, n2 %d\n", p1->E(), p2->E(),ncell1,ncell2);
1490 }
6175da48 1491 }
1492
1493 //---------
1494 // MC data
1495 //---------
1496 //Do some MC checks on the origin of the pair, is there any common ancestor and if there is one, who?
1497 if(IsDataMC()) FillMCVersusRecDataHistograms(p1->GetLabel(), p2->GetLabel(),p1->Pt(), p2->Pt(),ncell1, ncell2, m, pt, a,deta, dphi);
1498
1499 //-----------------------
1500 //Multi cuts analysis
1501 //-----------------------
1502 if(fMultiCutAna){
1503 //Histograms for different PID bits selection
1504 for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
1505
1506 if(p1->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton) &&
1507 p2->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton)) fhRePIDBits[ipid]->Fill(pt,m) ;
1508
1509 //printf("ipt %d, ipid%d, name %s\n",ipt, ipid, fhRePtPIDCuts[ipt*fNPIDBitsBits+ipid]->GetName());
1510 } // pid bit cut loop
1511
1512 //Several pt,ncell and asymmetry cuts
af7b3903 1513 for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
1514 for(Int_t icell=0; icell<fNCellNCuts; icell++){
1515 for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
1516 Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
1517 if(p1->Pt() > fPtCuts[ipt] && p2->Pt() > fPtCuts[ipt] &&
1518 a < fAsymCuts[iasym] &&
6175da48 1519 ncell1 >= fCellNCuts[icell] && ncell2 >= fCellNCuts[icell]){
1520 fhRePtNCellAsymCuts[index]->Fill(pt,m) ;
1521 //printf("ipt %d, icell%d, iasym %d, name %s\n",ipt, icell, iasym, fhRePtNCellAsymCuts[((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym]->GetName());
1522 if(module1==module2){
1523 if (module1==0) fhRePtNCellAsymCutsSM0[index]->Fill(pt,m) ;
1524 else if(module1==1) fhRePtNCellAsymCutsSM1[index]->Fill(pt,m) ;
1525 else if(module1==2) fhRePtNCellAsymCutsSM2[index]->Fill(pt,m) ;
1526 else if(module1==3) fhRePtNCellAsymCutsSM3[index]->Fill(pt,m) ;
1527 else printf("AliAnaPi0::FillHistograms() - WRONG SM NUMBER\n");
1528 }
1529 }
af7b3903 1530 }// pid bit cut loop
1531 }// icell loop
1532 }// pt cut loop
1533 for(Int_t iasym = 0; iasym < fNAsymCuts; iasym++){
1534 if(a < fAsymCuts[iasym])fhRePtMult[iasym]->Fill(pt,GetTrackMultiplicity(),m) ;
1535 }
af7b3903 1536 }// multiple cuts analysis
1537 }// ok if same sm
7e7694bb 1538 }// second same event particle
1539 }// first cluster
6175da48 1540
1541 //-------------------------------------------------------------
1542 // Mixing
1543 //-------------------------------------------------------------
7e7694bb 1544 if(fDoOwnMix){
6175da48 1545 //Recover events in with same characteristics as the current event
5025c139 1546 TList * evMixList=fEventsList[curCentrBin*GetNZvertBin()*GetNRPBin()+curZvertBin*GetNRPBin()+curRPBin] ;
7e7694bb 1547 Int_t nMixed = evMixList->GetSize() ;
1548 for(Int_t ii=0; ii<nMixed; ii++){
1549 TClonesArray* ev2= (TClonesArray*) (evMixList->At(ii));
1550 Int_t nPhot2=ev2->GetEntriesFast() ;
1551 Double_t m = -999;
6175da48 1552 if(GetDebug() > 1)
1553 printf("AliAnaPi0::MakeAnalysisFillHistograms() - Mixed event %d photon entries %d, centrality bin %d\n", ii, nPhot2, curCentrBin);
7e7694bb 1554
6175da48 1555 //---------------------------------
1556 //First loop on photons/clusters
1557 //---------------------------------
7e7694bb 1558 for(Int_t i1=0; i1<nPhot; i1++){
1559 AliAODPWG4Particle * p1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i1)) ;
6175da48 1560 if(fSameSM && GetModuleNumber(p1)!=module1) continue;
1561
1562 //Get kinematics of cluster and (super) module of this cluster
7e7694bb 1563 TLorentzVector photon1(p1->Px(),p1->Py(),p1->Pz(),p1->E());
af7b3903 1564 module1 = GetModuleNumber(p1);
6175da48 1565
1566 //---------------------------------
1567 //First loop on photons/clusters
1568 //---------------------------------
7e7694bb 1569 for(Int_t i2=0; i2<nPhot2; i2++){
1570 AliAODPWG4Particle * p2 = (AliAODPWG4Particle*) (ev2->At(i2)) ;
1571
6175da48 1572 //Get kinematics of second cluster and calculate those of the pair
7e7694bb 1573 TLorentzVector photon2(p2->Px(),p2->Py(),p2->Pz(),p2->E());
6175da48 1574 m = (photon1+photon2).M() ;
7e7694bb 1575 Double_t pt = (photon1 + photon2).Pt();
1576 Double_t a = TMath::Abs(p1->E()-p2->E())/(p1->E()+p2->E()) ;
1577
1578 //Check if opening angle is too large or too small compared to what is expected
1579 Double_t angle = photon1.Angle(photon2.Vect());
6175da48 1580 if(fUseAngleEDepCut && !GetNeutralMesonSelection()->IsAngleInWindow((photon1+photon2).E(),angle+0.05)){
1581 if(GetDebug() > 2)
1582 printf("AliAnaPi0::MakeAnalysisFillHistograms() -Mix pair angle %f not in E %f window\n",angle, (photon1+photon2).E());
1583 continue;
1584 }
1585 if(fUseAngleCut && (angle < fAngleCut || angle > fAngleMaxCut)) {
1586 if(GetDebug() > 2)
1587 printf("AliAnaPi0::MakeAnalysisFillHistograms() -Mix pair angle %f < cut %f\n",angle,fAngleCut);
1588 continue;
1589
1590 }
7e7694bb 1591
1592 if(GetDebug() > 2)
1593 printf("AliAnaPi0::MakeAnalysisFillHistograms() - Mixed Event: pT: photon1 %2.2f, photon2 %2.2f; Pair: pT %2.2f, mass %2.3f, a %f2.3\n",
af7b3903 1594 p1->Pt(), p2->Pt(), pt,m,a);
6175da48 1595
af7b3903 1596 //In case we want only pairs in same (super) module, check their origin.
1597 module2 = GetModuleNumber(p2);
6175da48 1598
1599 //-------------------------------------------------------------------------------------------------
1600 //Fill module dependent histograms, put a cut on assymmetry on the first available cut in the array
1601 //-------------------------------------------------------------------------------------------------
1602 if(a < fAsymCuts[0]){
1603 if(module1==module2 && module1 >=0 && module1<fNModules)
1604 fhMiMod[module1]->Fill(pt,m) ;
1605 else
1606 fhMiDiffMod[fNModules+2]->Fill(pt,m) ;
1607
1608 if(fCalorimeter=="EMCAL"){
1609 if((module1==0 && module2==2) || (module1==2 && module2==0)) fhMiDiffMod[0]->Fill(pt,m) ;
1610 if((module1==1 && module2==3) || (module1==3 && module2==1)) fhMiDiffMod[1]->Fill(pt,m) ;
1611 if((module1==0 && module2==1) || (module1==1 && module2==0)) fhMiDiffMod[2]->Fill(pt,m) ;
1612 if((module1==2 && module2==3) || (module1==3 && module2==2)) fhMiDiffMod[3]->Fill(pt,m) ;
1613 if((module1==0 && module2==3) || (module1==3 && module2==0)) fhMiDiffMod[4]->Fill(pt,m) ;
1614 if((module1==2 && module2==1) || (module1==1 && module2==2)) fhMiDiffMod[5]->Fill(pt,m) ;
1615
1616 }
1617 else {
1618 if((module1==0 && module2==1) || (module1==1 && module2==0)) fhMiDiffMod[0]->Fill(pt,m) ;
1619 if((module1==0 && module2==2) || (module1==2 && module2==0)) fhMiDiffMod[1]->Fill(pt,m) ;
1620 if((module1==1 && module2==2) || (module1==2 && module2==1)) fhMiDiffMod[2]->Fill(pt,m) ;
1621 }
1622 }
1623
af7b3903 1624 Bool_t ok = kTRUE;
1625 if(fSameSM && module1!=module2) ok=kFALSE;
1626 if(ok){
6175da48 1627
1628 //Check if one of the clusters comes from a conversion
1629 if (p1->IsTagged() && p2->IsTagged()) fhMiConv2->Fill(pt,m);
1630 else if(p1->IsTagged() || p2->IsTagged()) fhMiConv ->Fill(pt,m);
1631
1632 //Fill histograms for different bad channel distance, centrality, assymmetry cut and pid bit
af7b3903 1633 for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
1634 if((p1->IsPIDOK(ipid,AliCaloPID::kPhoton)) && (p2->IsPIDOK(ipid,AliCaloPID::kPhoton))){
1635 for(Int_t iasym=0; iasym < fNAsymCuts; iasym++){
1636 if(a < fAsymCuts[iasym]){
1637 Int_t index = ((curCentrBin*fNPIDBits)+ipid)*fNAsymCuts + iasym;
1638 fhMi1 [index]->Fill(pt,m) ;
398c93cc 1639 if(fMakeInvPtPlots)fhMiInvPt1[index]->Fill(pt,m,1./pt) ;
6175da48 1640 if(fFillBadDistHisto){
1641 if(p1->DistToBad()>0 && p2->DistToBad()>0){
1642 fhMi2 [index]->Fill(pt,m) ;
1643 if(fMakeInvPtPlots)fhMiInvPt2[index]->Fill(pt,m,1./pt) ;
1644 if(p1->DistToBad()>1 && p2->DistToBad()>1){
1645 fhMi3 [index]->Fill(pt,m) ;
1646 if(fMakeInvPtPlots)fhMiInvPt3[index]->Fill(pt,m,1./pt) ;
1647 }
af7b3903 1648 }
6175da48 1649 }// Fill bad dist histo
af7b3903 1650 }//Asymmetry cut
1651 }// Asymmetry loop
1652 }//PID cut
1653 }//loop for histograms
6175da48 1654
1655 //-----------------------
1656 //Multi cuts analysis
1657 //-----------------------
1658 if(fMultiCutAna){
1659 //Several pt,ncell and asymmetry cuts
1660 for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
1661 for(Int_t icell=0; icell<fNCellNCuts; icell++){
1662 for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
1663 Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
1664 if(p1->Pt() > fPtCuts[ipt] && p2->Pt() > fPtCuts[ipt] &&
1665 a < fAsymCuts[iasym] &&
1666 p1->GetBtag() >= fCellNCuts[icell] && p2->GetBtag() >= fCellNCuts[icell]){
1667 fhMiPtNCellAsymCuts[index]->Fill(pt,m) ;
1668 //printf("ipt %d, icell%d, iasym %d, name %s\n",ipt, icell, iasym, fhRePtNCellAsymCuts[((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym]->GetName());
1669 }
1670 }// pid bit cut loop
1671 }// icell loop
1672 }// pt cut loop
1673 } // Multi cut ana
1674
1675 //Fill histograms with opening angle
1676 fhMixedOpeningAngle ->Fill(pt,angle);
1677 fhMixedCosOpeningAngle->Fill(pt,TMath::Cos(angle));
af7b3903 1678 }//ok
7e7694bb 1679 }// second cluster loop
1680 }//first cluster loop
1681 }//loop on mixed events
1682
6175da48 1683 //--------------------------------------------------------
1684 //Add the current event to the list of events for mixing
1685 //--------------------------------------------------------
7e7694bb 1686 TClonesArray *currentEvent = new TClonesArray(*GetInputAODBranch());
af7b3903 1687 //Add current event to buffer and Remove redundant events
7e7694bb 1688 if(currentEvent->GetEntriesFast()>0){
1689 evMixList->AddFirst(currentEvent) ;
1690 currentEvent=0 ; //Now list of particles belongs to buffer and it will be deleted with buffer
1691 if(evMixList->GetSize()>=fNmaxMixEv)
1692 {
1693 TClonesArray * tmp = (TClonesArray*) (evMixList->Last()) ;
1694 evMixList->RemoveLast() ;
1695 delete tmp ;
1696 }
1697 }
1698 else{ //empty event
1699 delete currentEvent ;
1700 currentEvent=0 ;
477d6cee 1701 }
7e7694bb 1702 }// DoOwnMix
c8fe2783 1703
1c5acb87 1704}
1705
a5cc4f03 1706//________________________________________________________________________
1707void AliAnaPi0::ReadHistograms(TList* outputList)
1708{
50f39b97 1709 // Needed when Terminate is executed in distributed environment
1710 // Refill analysis histograms of this class with corresponding histograms in output list.
1711
1712 // Histograms of this analsys are kept in the same list as other analysis, recover the position of
1713 // the first one and then add the next.
1714 Int_t index = outputList->IndexOf(outputList->FindObject(GetAddedHistogramsStringToName()+"hRe_cen0_pid0_dist1"));
1715
af7b3903 1716 if(!fhRe1) fhRe1 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
1717 if(!fhRe2) fhRe2 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
1718 if(!fhRe3) fhRe3 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
1719 if(!fhMi1) fhMi1 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
1720 if(!fhMi2) fhMi2 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
1721 if(!fhMi3) fhMi3 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
6175da48 1722 if(!fhReInvPt1) fhReInvPt1 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
1723 if(!fhReInvPt2) fhReInvPt2 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
1724 if(!fhReInvPt3) fhReInvPt3 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
1725 if(!fhMiInvPt1) fhMiInvPt1 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
1726 if(!fhMiInvPt2) fhMiInvPt2 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
1727 if(!fhMiInvPt3) fhMiInvPt3 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
af7b3903 1728 if(!fhReMod) fhReMod = new TH2D*[fNModules] ;
1729 if(!fhReDiffMod)fhReDiffMod = new TH2D*[fNModules+1] ;
6175da48 1730 if(!fhMiMod) fhReMod = new TH2D*[fNModules] ;
1731 if(!fhMiDiffMod)fhReDiffMod = new TH2D*[fNModules+1] ;
1732
1733 fhReConv = (TH2D*) outputList->At(index++);
1734 fhMiConv = (TH2D*) outputList->At(index++);
1735 fhReConv2 = (TH2D*) outputList->At(index++);
1736 fhMiConv2 = (TH2D*) outputList->At(index++);
821c8090 1737
50f39b97 1738 for(Int_t ic=0; ic<fNCentrBin; ic++){
af7b3903 1739 for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
1740 for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
1741 Int_t ihisto = ((ic*fNPIDBits)+ipid)*fNAsymCuts + iasym;
1742
1743 fhRe1[ihisto] = (TH2D*) outputList->At(index++);
1744 fhRe2[ihisto] = (TH2D*) outputList->At(index++);
1745 fhRe3[ihisto] = (TH2D*) outputList->At(index++);
6175da48 1746
1747 fhReInvPt1[ihisto] = (TH2D*) outputList->At(index++);
1748 fhReInvPt2[ihisto] = (TH2D*) outputList->At(index++);
1749 fhReInvPt3[ihisto] = (TH2D*) outputList->At(index++);
5ae09196 1750
af7b3903 1751 if(fDoOwnMix){
1752 fhMi1[ihisto] = (TH2D*) outputList->At(index++);
1753 fhMi2[ihisto] = (TH2D*) outputList->At(index++);
1754 fhMi3[ihisto] = (TH2D*) outputList->At(index++);
6175da48 1755
1756 fhMiInvPt1[ihisto] = (TH2D*) outputList->At(index++);
1757 fhMiInvPt2[ihisto] = (TH2D*) outputList->At(index++);
1758 fhMiInvPt3[ihisto] = (TH2D*) outputList->At(index++);
af7b3903 1759 }//Own mix
1760 }//asymmetry loop
1761 }// pid loop
1762 }// centrality loop
1763
1764 fhRePtAsym = (TH2D*)outputList->At(index++);
1765 fhRePtAsymPi0 = (TH2D*)outputList->At(index++);
1766 fhRePtAsymEta = (TH2D*)outputList->At(index++);
eee5fcf1 1767
5ae09196 1768 if(fMultiCutAna){
1769
eee5fcf1 1770 if(!fhRePtNCellAsymCuts) fhRePtNCellAsymCuts = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
1771 if(!fhRePIDBits) fhRePIDBits = new TH2D*[fNPIDBits];
1772
5ae09196 1773 for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
1774 fhRePIDBits[ipid] = (TH2D*) outputList->At(index++);
1775 }// ipid loop
1776
1777 for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
1778 for(Int_t icell=0; icell<fNCellNCuts; icell++){
1779 for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
1780 fhRePtNCellAsymCuts[((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym] = (TH2D*) outputList->At(index++);
1781 }// iasym
1782 }// icell loop
1783 }// ipt loop
af7b3903 1784
1785 if(!fhRePtMult) fhRePtMult = new TH3D*[fNAsymCuts] ;
1786 for(Int_t iasym = 0; iasym < fNAsymCuts; iasym++)
1787 fhRePtMult[iasym] = (TH3D*) outputList->At(index++);
5ae09196 1788 }// multi cut analysis
50f39b97 1789
1790 fhEvents = (TH3D *) outputList->At(index++);
1791
af7b3903 1792 fhRealOpeningAngle = (TH2D*) outputList->At(index++);
1793 fhRealCosOpeningAngle = (TH2D*) outputList->At(index++);
6175da48 1794 if(fDoOwnMix){
1795 fhMixedOpeningAngle = (TH2D*) outputList->At(index++);
1796 fhMixedCosOpeningAngle = (TH2D*) outputList->At(index++);
1797 }
af7b3903 1798
50f39b97 1799 //Histograms filled only if MC data is requested
1800 if(IsDataMC() || (GetReader()->GetDataType() == AliCaloTrackReader::kMC) ){
6175da48 1801 fhPrimPi0Pt = (TH1D*) outputList->At(index++);
1802 fhPrimPi0AccPt = (TH1D*) outputList->At(index++);
1803 fhPrimPi0Y = (TH1D*) outputList->At(index++);
1804 fhPrimPi0AccY = (TH1D*) outputList->At(index++);
1805 fhPrimPi0Phi = (TH1D*) outputList->At(index++);
1806 fhPrimPi0AccPhi = (TH1D*) outputList->At(index++);
1807 for(Int_t i = 0; i<13; i++){
1808 fhMCOrgMass[i] = (TH2D*) outputList->At(index++);
1809 fhMCOrgAsym[i] = (TH2D*) outputList->At(index++);
1810 fhMCOrgDeltaEta[i] = (TH2D*) outputList->At(index++);
1811 fhMCOrgDeltaPhi[i] = (TH2D*) outputList->At(index++);
1812 }
1813
1814 if(fMultiCutAnaSim){
1815 fhMCPi0MassPtTrue = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
1816 fhMCPi0MassPtRec = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
1817 fhMCPi0PtTruePtRec = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
1818 fhMCEtaMassPtTrue = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
1819 fhMCEtaMassPtRec = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
1820 fhMCEtaPtTruePtRec = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
1821 for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
1822 for(Int_t icell=0; icell<fNCellNCuts; icell++){
1823 for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
1824 Int_t in = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
1825 fhMCPi0MassPtTrue[in] = (TH2D*) outputList->At(index++);
1826 fhMCPi0PtTruePtRec[in] = (TH2D*) outputList->At(index++);
1827 fhMCEtaMassPtTrue[in] = (TH2D*) outputList->At(index++);
1828 fhMCEtaPtTruePtRec[in] = (TH2D*) outputList->At(index++);
1829 }
1830 }
1831 }
1832 }
1833 else{
1834 fhMCPi0MassPtTrue = new TH2D*[1];
1835 fhMCPi0PtTruePtRec = new TH2D*[1];
1836 fhMCEtaMassPtTrue = new TH2D*[1];
1837 fhMCEtaPtTruePtRec = new TH2D*[1];
1838
1839 fhMCPi0MassPtTrue[0] = (TH2D*) outputList->At(index++);
1840 fhMCPi0PtTruePtRec[0] = (TH2D*) outputList->At(index++);
1841 fhMCEtaMassPtTrue[0] = (TH2D*) outputList->At(index++);
1842 fhMCEtaPtTruePtRec[0] = (TH2D*) outputList->At(index++);
1843 }
50f39b97 1844 }
1845
6175da48 1846 for(Int_t imod=0; imod < fNModules; imod++){
1847 fhReMod[imod] = (TH2D*) outputList->At(index++);
1848 fhReDiffMod[imod] = (TH2D*) outputList->At(index++);
1849 if(fDoOwnMix){
1850 fhMiMod[imod] = (TH2D*) outputList->At(index++);
1851 fhMiDiffMod[imod] = (TH2D*) outputList->At(index++);
1852 }
1853 }
eee5fcf1 1854
a5cc4f03 1855}
1856
1857
6639984f 1858//____________________________________________________________________________________________________________________________________________________
a5cc4f03 1859void AliAnaPi0::Terminate(TList* outputList)
6639984f 1860{
1861 //Do some calculations and plots from the final histograms.
477d6cee 1862
fbeaf916 1863 printf(" *** %s Terminate:\n", GetName()) ;
50f39b97 1864
a5cc4f03 1865 //Recover histograms from output histograms list, needed for distributed analysis.
1866 ReadHistograms(outputList);
50f39b97 1867
2e557d1c 1868 if (!fhRe1) {
50f39b97 1869 printf("AliAnaPi0::Terminate() - Error: Remote output histograms not imported in AliAnaPi0 object");
1870 return;
2e557d1c 1871 }
50f39b97 1872
a3aebfff 1873 printf("AliAnaPi0::Terminate() Mgg Real : %5.3f , RMS : %5.3f \n", fhRe1[0]->GetMean(), fhRe1[0]->GetRMS() ) ;
5ae09196 1874
1875 const Int_t buffersize = 255;
1876
1877 char nameIM[buffersize];
1878 snprintf(nameIM, buffersize,"AliAnaPi0_%s_cPt",fCalorimeter.Data());
71dd883b 1879 TCanvas * cIM = new TCanvas(nameIM, "", 400, 10, 600, 700) ;
6639984f 1880 cIM->Divide(2, 2);
50f39b97 1881
6639984f 1882 cIM->cd(1) ;
1883 //gPad->SetLogy();
af7b3903 1884 TH1D * hIMAllPt = (TH1D*) fhRe1[0]->ProjectionY(Form("IMPtAll_%s",fCalorimeter.Data()));
6639984f 1885 hIMAllPt->SetLineColor(2);
1886 hIMAllPt->SetTitle("No cut on p_{T, #gamma#gamma} ");
1887 hIMAllPt->Draw();
1888
1889 cIM->cd(2) ;
af7b3903 1890 TH1D * hIMPt5 = (TH1D*) fhRe1[0]->ProjectionY(Form("IMPt0-5_%s",fCalorimeter.Data()),0, fhRe1[0]->GetXaxis()->FindBin(5.));
2244659d 1891// hRe1Pt5->GetXaxis()->SetRangeUser(0,5);
1892// TH1D * hIMPt5 = (TH1D*) hRe1Pt5->Project3D(Form("IMPt5_%s_pz",fCalorimeter.Data()));
6639984f 1893 hIMPt5->SetLineColor(2);
1894 hIMPt5->SetTitle("0 < p_{T, #gamma#gamma} < 5 GeV/c");
1895 hIMPt5->Draw();
1896
1897 cIM->cd(3) ;
af7b3903 1898 TH1D * hIMPt10 = (TH1D*) fhRe1[0]->ProjectionY(Form("IMPt5-10_%s",fCalorimeter.Data()), fhRe1[0]->GetXaxis()->FindBin(5.),fhRe1[0]->GetXaxis()->FindBin(10.));
2244659d 1899// hRe1Pt10->GetXaxis()->SetRangeUser(5,10);
1900// TH1D * hIMPt10 = (TH1D*) hRe1Pt10->Project3D(Form("IMPt10_%s_pz",fCalorimeter.Data()));
6639984f 1901 hIMPt10->SetLineColor(2);
1902 hIMPt10->SetTitle("5 < p_{T, #gamma#gamma} < 10 GeV/c");
1903 hIMPt10->Draw();
1904
1905 cIM->cd(4) ;
af7b3903 1906 TH1D * hIMPt20 = (TH1D*) fhRe1[0]->ProjectionY(Form("IMPt10-20_%s",fCalorimeter.Data()), fhRe1[0]->GetXaxis()->FindBin(10.),fhRe1[0]->GetXaxis()->FindBin(20.));
2244659d 1907 // TH3F * hRe1Pt20 = (TH3F*)fhRe1[0]->Clone(Form("IMPt20_%s",fCalorimeter.Data()));
1908// hRe1Pt20->GetXaxis()->SetRangeUser(10,20);
1909// TH1D * hIMPt20 = (TH1D*) hRe1Pt20->Project3D(Form("IMPt20_%s_pz",fCalorimeter.Data()));
6639984f 1910 hIMPt20->SetLineColor(2);
1911 hIMPt20->SetTitle("10 < p_{T, #gamma#gamma} < 20 GeV/c");
1912 hIMPt20->Draw();
1913
5ae09196 1914 char nameIMF[buffersize];
1915 snprintf(nameIMF,buffersize,"AliAnaPi0_%s_Mgg.eps",fCalorimeter.Data());
71dd883b 1916 cIM->Print(nameIMF);
6639984f 1917
5ae09196 1918 char namePt[buffersize];
1919 snprintf(namePt,buffersize,"AliAnaPi0_%s_cPt",fCalorimeter.Data());
71dd883b 1920 TCanvas * cPt = new TCanvas(namePt, "", 400, 10, 600, 700) ;
6639984f 1921 cPt->Divide(2, 2);
1922
1923 cPt->cd(1) ;
1924 //gPad->SetLogy();
af7b3903 1925 TH1D * hPt = (TH1D*) fhRe1[0]->ProjectionX(Form("Pt0_%s",fCalorimeter.Data()),-1,-1);
6639984f 1926 hPt->SetLineColor(2);
1927 hPt->SetTitle("No cut on M_{#gamma#gamma} ");
1928 hPt->Draw();
1929
1930 cPt->cd(2) ;
af7b3903 1931 TH1D * hPtIM1 = (TH1D*)fhRe1[0]->ProjectionX(Form("Pt1_%s",fCalorimeter.Data()), fhRe1[0]->GetZaxis()->FindBin(0.05),fhRe1[0]->GetZaxis()->FindBin(0.21));
2244659d 1932// TH3F * hRe1IM1 = (TH3F*)fhRe1[0]->Clone(Form("Pt1_%s",fCalorimeter.Data()));
1933// hRe1IM1->GetZaxis()->SetRangeUser(0.05,0.21);
1934// TH1D * hPtIM1 = (TH1D*) hRe1IM1->Project3D("x");
6639984f 1935 hPtIM1->SetLineColor(2);
1936 hPtIM1->SetTitle("0.05 < M_{#gamma#gamma} < 0.21 GeV/c^{2}");
1937 hPtIM1->Draw();
1938
1939 cPt->cd(3) ;
af7b3903 1940 TH1D * hPtIM2 = (TH1D*)fhRe1[0]->ProjectionX(Form("Pt2_%s",fCalorimeter.Data()), fhRe1[0]->GetZaxis()->FindBin(0.09),fhRe1[0]->GetZaxis()->FindBin(0.17));
2244659d 1941// TH3F * hRe1IM2 = (TH3F*)fhRe1[0]->Clone(Form("Pt2_%s",fCalorimeter.Data()));
1942// hRe1IM2->GetZaxis()->SetRangeUser(0.09,0.17);
1943// TH1D * hPtIM2 = (TH1D*) hRe1IM2->Project3D("x");
6639984f 1944 hPtIM2->SetLineColor(2);
1945 hPtIM2->SetTitle("0.09 < M_{#gamma#gamma} < 0.17 GeV/c^{2}");
1946 hPtIM2->Draw();
1947
1948 cPt->cd(4) ;
af7b3903 1949 TH1D * hPtIM3 = (TH1D*)fhRe1[0]->ProjectionX(Form("Pt3_%s",fCalorimeter.Data()), fhRe1[0]->GetZaxis()->FindBin(0.11),fhRe1[0]->GetZaxis()->FindBin(0.15));
2244659d 1950// TH3F * hRe1IM3 = (TH3F*)fhRe1[0]->Clone(Form("Pt3_%s",fCalorimeter.Data()));
1951// hRe1IM3->GetZaxis()->SetRangeUser(0.11,0.15);
1952// TH1D * hPtIM3 = (TH1D*) hRe1IM1->Project3D("x");
6639984f 1953 hPtIM3->SetLineColor(2);
1954 hPtIM3->SetTitle("0.11 < M_{#gamma#gamma} < 0.15 GeV/c^{2}");
1955 hPtIM3->Draw();
1956
164a1d84 1957 char namePtF[buffersize];
5ae09196 1958 snprintf(namePtF,buffersize,"AliAnaPi0_%s_Pt.eps",fCalorimeter.Data());
71dd883b 1959 cPt->Print(namePtF);
1c5acb87 1960
5ae09196 1961 char line[buffersize] ;
1962 snprintf(line,buffersize,".!tar -zcf %s_%s.tar.gz *.eps", GetName(),fCalorimeter.Data()) ;
6639984f 1963 gROOT->ProcessLine(line);
5ae09196 1964 snprintf(line, buffersize,".!rm -fR AliAnaPi0_%s*.eps",fCalorimeter.Data());
6639984f 1965 gROOT->ProcessLine(line);
1966
71dd883b 1967 printf(" AliAnaPi0::Terminate() - !! All the eps files are in %s_%s.tar.gz !!!\n", GetName(), fCalorimeter.Data());
1c5acb87 1968
6639984f 1969}
c8fe2783 1970 //____________________________________________________________________________________________________________________________________________________
1971Int_t AliAnaPi0::GetEventIndex(AliAODPWG4Particle * part, Double_t * vert)
1972{
f8006433 1973 // retieves the event index and checks the vertex
1974 // in the mixed buffer returns -2 if vertex NOK
1975 // for normal events returns 0 if vertex OK and -1 if vertex NOK
1976
1977 Int_t evtIndex = -1 ;
1978 if(GetReader()->GetDataType()!=AliCaloTrackReader::kMC){
1979
1980 if (GetMixedEvent()){
1981
1982 evtIndex = GetMixedEvent()->EventIndexForCaloCluster(part->GetCaloLabel(0)) ;
1983 GetVertex(vert,evtIndex);
1984
5025c139 1985 if(TMath::Abs(vert[2])> GetZvertexCut())
f8006433 1986 evtIndex = -2 ; //Event can not be used (vertex, centrality,... cuts not fulfilled)
1987 } else {// Single event
1988
1989 GetVertex(vert);
1990
5025c139 1991 if(TMath::Abs(vert[2])> GetZvertexCut())
f8006433 1992 evtIndex = -1 ; //Event can not be used (vertex, centrality,... cuts not fulfilled)
1993 else
1994 evtIndex = 0 ;
c8fe2783 1995 }
0ae57829 1996 }//No MC reader
f8006433 1997 else {
1998 evtIndex = 0;
1999 vert[0] = 0. ;
2000 vert[1] = 0. ;
2001 vert[2] = 0. ;
2002 }
0ae57829 2003
f8006433 2004 return evtIndex ;
c8fe2783 2005}
f8006433 2006