]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG4/PartCorrDep/AliAnaPi0.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Set the analysis parameters/cuts used in a folder in the output file, in TObjString...
[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
19// extractin raw pi0 yield.
20//
21//-- Author: Dmitri Peressounko (RRC "KI")
22//-- Adapted to PartCorr frame by Lamia Benhabib (SUBATECH)
23//-- and Gustavo Conesa (INFN-Frascati)
24//_________________________________________________________________________
25
26
27// --- ROOT system ---
28#include "TH3.h"
50f39b97 29#include "TH2D.h"
1c5acb87 30//#include "Riostream.h"
6639984f 31#include "TCanvas.h"
32#include "TPad.h"
33#include "TROOT.h"
477d6cee 34#include "TClonesArray.h"
0c1383b5 35#include "TObjString.h"
1c5acb87 36
37//---- AliRoot system ----
38#include "AliAnaPi0.h"
39#include "AliCaloTrackReader.h"
40#include "AliCaloPID.h"
6639984f 41#include "AliStack.h"
ff45398a 42#include "AliFiducialCut.h"
477d6cee 43#include "TParticle.h"
44#include "AliAODCaloCluster.h"
45#include "AliVEvent.h"
6921fa00 46#include "AliESDCaloCluster.h"
47#include "AliESDEvent.h"
48#include "AliAODEvent.h"
50f39b97 49#include "AliNeutralMesonSelection.h"
591cc579 50
1c5acb87 51ClassImp(AliAnaPi0)
52
53//________________________________________________________________________________________________________________________________________________
54AliAnaPi0::AliAnaPi0() : AliAnaPartCorrBaseClass(),
55fNCentrBin(0),fNZvertBin(0),fNrpBin(0),
56fNPID(0),fNmaxMixEv(0), fZvtxCut(0.),fCalorimeter(""),
50f39b97 57fNModules(12), fUseAngleCut(kFALSE), fEventsList(0x0), //fhEtalon(0x0),
5a2dbc3c 58fhReMod(0x0), fhRe1(0x0),fhMi1(0x0),fhRe2(0x0),fhMi2(0x0),fhRe3(0x0),fhMi3(0x0),fhEvents(0x0),
50f39b97 59fhRealOpeningAngle(0x0),fhRealCosOpeningAngle(0x0),
60fhPrimPt(0x0), fhPrimAccPt(0x0), fhPrimY(0x0), fhPrimAccY(0x0), fhPrimPhi(0x0), fhPrimAccPhi(0x0),
61fhPrimOpeningAngle(0x0),fhPrimCosOpeningAngle(0x0)
1c5acb87 62{
63//Default Ctor
64 InitParameters();
65
66}
78219bac 67/*
1c5acb87 68//________________________________________________________________________________________________________________________________________________
69AliAnaPi0::AliAnaPi0(const AliAnaPi0 & ex) : AliAnaPartCorrBaseClass(ex),
70fNCentrBin(ex.fNCentrBin),fNZvertBin(ex.fNZvertBin),fNrpBin(ex.fNrpBin),
c0649ccf 71fNPID(ex.fNPID),fNmaxMixEv(ex.fNmaxMixEv),fZvtxCut(ex.fZvtxCut), fCalorimeter(ex.fCalorimeter),
614701c6 72fNModules(ex.fNModules), fUseAngleCut(ex.fUseAngleCut), fEventsList(0x0), //fhEtalon(ex.fhEtalon),
5a2dbc3c 73fhReMod(ex.fhReMod), fhRe1(ex.fhRe1),fhMi1(ex.fhMi1),fhRe2(ex.fhRe2),fhMi2(ex.fhMi2),
74fhRe3(ex.fhRe3),fhMi3(ex.fhMi3),fhEvents(ex.fhEvents),
50f39b97 75fhRealOpeningAngle(ex.fhRealOpeningAngle),fhRealCosOpeningAngle(ex.fhRealCosOpeningAngle),
1c5acb87 76fhPrimPt(ex.fhPrimPt), fhPrimAccPt(ex.fhPrimAccPt), fhPrimY(ex.fhPrimY),
50f39b97 77fhPrimAccY(ex.fhPrimAccY), fhPrimPhi(ex.fhPrimPhi), fhPrimAccPhi(ex.fhPrimAccPhi),
78fhPrimOpeningAngle(ex.fhPrimOpeningAngle),fhPrimCosOpeningAngle(ex.fhPrimCosOpeningAngle)
1c5acb87 79{
477d6cee 80 // cpy ctor
81 //Do not need it
1c5acb87 82}
78219bac 83*/
614701c6 84//
85////________________________________________________________________________________________________________________________________________________
86//AliAnaPi0 & AliAnaPi0::operator = (const AliAnaPi0 & ex)
87//{
88// // assignment operator
89//
90// if(this == &ex)return *this;
91// ((AliAnaPartCorrBaseClass *)this)->operator=(ex);
92//
93// fNCentrBin = ex.fNCentrBin ; fNZvertBin = ex.fNZvertBin ; fNrpBin = ex.fNrpBin ;
94// fNPID = ex.fNPID ; fNmaxMixEv = ex.fNmaxMixEv ; fZvtxCut = ex.fZvtxCut ; fCalorimeter = ex.fCalorimeter ;
95// fNModules = ex.fNModules; fEventsList = ex.fEventsList ; //fhEtalon = ex.fhEtalon ;
96// fhRe1 = ex.fhRe1 ; fhMi1 = ex.fhMi1 ; fhRe2 = ex.fhRe2 ; fhMi2 = ex.fhMi2 ; fhReMod = ex.fhReMod;
97// fhRe3 = ex.fhRe3 ; fhMi3 = ex.fhMi3 ; fhEvents = ex.fhEvents ; fUseAngleCut = ex.fUseAngleCut;
98// fhPrimPt = ex.fhPrimPt ; fhPrimAccPt = ex.fhPrimAccPt ; fhPrimY = ex.fhPrimY ;
99// fhPrimAccY = ex.fhPrimAccY ; fhPrimPhi = ex.fhPrimPhi ; fhPrimAccPhi = ex.fhPrimAccPhi ;
100// fhRealOpeningAngle = ex.fhRealOpeningAngle; fhRealCosOpeningAngle = ex.fhRealCosOpeningAngle;
101// fhPrimOpeningAngle = ex.fhPrimOpeningAngle; fhPrimCosOpeningAngle = ex.fhPrimCosOpeningAngle;
102//
103// return *this;
104//
105//}
1c5acb87 106
107//________________________________________________________________________________________________________________________________________________
108AliAnaPi0::~AliAnaPi0() {
477d6cee 109 // Remove event containers
110 if(fEventsList){
111 for(Int_t ic=0; ic<fNCentrBin; ic++){
112 for(Int_t iz=0; iz<fNZvertBin; iz++){
113 for(Int_t irp=0; irp<fNrpBin; irp++){
114 fEventsList[ic*fNZvertBin*fNrpBin+iz*fNrpBin+irp]->Delete() ;
115 delete fEventsList[ic*fNZvertBin*fNrpBin+iz*fNrpBin+irp] ;
1c5acb87 116 }
477d6cee 117 }
118 }
119 delete[] fEventsList;
120 fEventsList=0 ;
121 }
591cc579 122
1c5acb87 123}
124
125//________________________________________________________________________________________________________________________________________________
126void AliAnaPi0::InitParameters()
127{
128//Init parameters when first called the analysis
129//Set default parameters
a3aebfff 130 SetInputAODName("PWG4Particle");
131
132 AddToHistogramsName("AnaPi0_");
6921fa00 133 fNModules = 12; // set maximum to maximum number of EMCAL modules
477d6cee 134 fNCentrBin = 1;
135 fNZvertBin = 1;
136 fNrpBin = 1;
137 fNPID = 9;
138 fNmaxMixEv = 10;
139 fZvtxCut = 40;
140 fCalorimeter = "PHOS";
50f39b97 141 fUseAngleCut = kFALSE;
142
1c5acb87 143}
144//________________________________________________________________________________________________________________________________________________
5a2dbc3c 145//void AliAnaPi0::Init()
146//{
477d6cee 147 //Init some data members needed in analysis
148
149 //Histograms binning and range
5a2dbc3c 150 // if(!fhEtalon){ // p_T alpha d m_gg
151 // fhEtalon = new TH3D("hEtalon","Histo with binning parameters",50,0.,25.,10,0.,1.,200,0.,1.) ;
152 // fhEtalon->SetXTitle("P_{T} (GeV)") ;
153 // fhEtalon->SetYTitle("#alpha") ;
154 // fhEtalon->SetZTitle("m_{#gamma#gamma} (GeV)") ;
155 // }
477d6cee 156
5a2dbc3c 157//}
1c5acb87 158
0c1383b5 159
160//________________________________________________________________________________________________________________________________________________
161TObjString * AliAnaPi0::GetAnalysisCuts()
162{
163 //Save parameters used for analysis
164 TString parList ; //this will be list of parameters used for this analysis.
165 char onePar[255] ;
166 sprintf(onePar,"--- AliAnaPi0 ---\n") ;
167 parList+=onePar ;
168 sprintf(onePar,"Number of bins in Centrality: %d \n",fNCentrBin) ;
169 parList+=onePar ;
170 sprintf(onePar,"Number of bins in Z vert. pos: %d \n",fNZvertBin) ;
171 parList+=onePar ;
172 sprintf(onePar,"Number of bins in Reac. Plain: %d \n",fNrpBin) ;
173 parList+=onePar ;
174 sprintf(onePar,"Depth of event buffer: %d \n",fNmaxMixEv) ;
175 parList+=onePar ;
176 sprintf(onePar,"Number of different PID used: %d \n",fNPID) ;
177 parList+=onePar ;
178 sprintf(onePar,"Cuts: \n") ;
179 parList+=onePar ;
180 sprintf(onePar,"Z vertex position: -%f < z < %f \n",fZvtxCut,fZvtxCut) ;
181 parList+=onePar ;
182 sprintf(onePar,"Calorimeter: %s \n",fCalorimeter.Data()) ;
183 parList+=onePar ;
184 sprintf(onePar,"Number of modules: %d \n",fNModules) ;
185 parList+=onePar ;
186
187 return new TObjString(parList) ;
188}
189
2e557d1c 190//________________________________________________________________________________________________________________________________________________
191TList * AliAnaPi0::GetCreateOutputObjects()
192{
477d6cee 193 // Create histograms to be saved in output file and
194 // store them in fOutputContainer
195
196 //create event containers
197 fEventsList = new TList*[fNCentrBin*fNZvertBin*fNrpBin] ;
1c5acb87 198
477d6cee 199 for(Int_t ic=0; ic<fNCentrBin; ic++){
200 for(Int_t iz=0; iz<fNZvertBin; iz++){
201 for(Int_t irp=0; irp<fNrpBin; irp++){
202 fEventsList[ic*fNZvertBin*fNrpBin+iz*fNrpBin+irp] = new TList() ;
203 }
204 }
205 }
765d44e7 206
477d6cee 207 TList * outputContainer = new TList() ;
208 outputContainer->SetName(GetName());
6921fa00 209
210 fhReMod = new TH3D*[fNModules] ;
50305ae9 211 fhRe1 = new TH3D*[fNCentrBin*fNPID] ;
212 fhRe2 = new TH3D*[fNCentrBin*fNPID] ;
213 fhRe3 = new TH3D*[fNCentrBin*fNPID] ;
214 fhMi1 = new TH3D*[fNCentrBin*fNPID] ;
215 fhMi2 = new TH3D*[fNCentrBin*fNPID] ;
216 fhMi3 = new TH3D*[fNCentrBin*fNPID] ;
477d6cee 217
218 char key[255] ;
219 char title[255] ;
5a2dbc3c 220 Int_t nptbins = GetHistoPtBins();
221 Int_t nphibins = GetHistoPhiBins();
222 Int_t netabins = GetHistoEtaBins();
223 Float_t ptmax = GetHistoPtMax();
224 Float_t phimax = GetHistoPhiMax();
225 Float_t etamax = GetHistoEtaMax();
226 Float_t ptmin = GetHistoPtMin();
227 Float_t phimin = GetHistoPhiMin();
228 Float_t etamin = GetHistoEtaMin();
229
230 Int_t nmassbins = GetHistoMassBins();
231 Int_t nasymbins = GetHistoAsymmetryBins();
232 Float_t massmax = GetHistoMassMax();
233 Float_t asymmax = GetHistoAsymmetryMax();
234 Float_t massmin = GetHistoMassMin();
235 Float_t asymmin = GetHistoAsymmetryMin();
236
477d6cee 237 for(Int_t ic=0; ic<fNCentrBin; ic++){
238 for(Int_t ipid=0; ipid<fNPID; ipid++){
6921fa00 239
477d6cee 240 //Distance to bad module 1
241 sprintf(key,"hRe_cen%d_pid%d_dist1",ic,ipid) ;
242 sprintf(title,"Real m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
243
5a2dbc3c 244 //fhEtalon->Clone(key);
245 //fhRe1[ic*fNPID+ipid]=(TH3D*)fhEtalon->Clone(key) ;
246 //fhRe1[ic*fNPID+ipid]->SetName(key) ;
247 //fhRe1[ic*fNPID+ipid]->SetTitle(title) ;
248 fhRe1[ic*fNPID+ipid] = new TH3D(key,title,nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax,nmassbins,massmin,massmax) ;
477d6cee 249 outputContainer->Add(fhRe1[ic*fNPID+ipid]) ;
250
251 sprintf(key,"hMi_cen%d_pid%d_dist1",ic,ipid) ;
252 sprintf(title,"Mixed m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
5a2dbc3c 253 //fhMi1[ic*fNPID+ipid]=(TH3D*)fhEtalon->Clone(key) ;
254 //fhMi1[ic*fNPID+ipid]->SetName(key) ;
255 //fhMi1[ic*fNPID+ipid]->SetTitle(title) ;
256 fhMi1[ic*fNPID+ipid] = new TH3D(key,title,nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax,nmassbins,massmin,massmax) ;
477d6cee 257 outputContainer->Add(fhMi1[ic*fNPID+ipid]) ;
258
259 //Distance to bad module 2
260 sprintf(key,"hRe_cen%d_pid%d_dist2",ic,ipid) ;
261 sprintf(title,"Real m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
5a2dbc3c 262 //fhRe2[ic*fNPID+ipid]=(TH3D*)fhEtalon->Clone(key) ;
263 //fhRe2[ic*fNPID+ipid]->SetName(key) ;
264 //fhRe2[ic*fNPID+ipid]->SetTitle(title) ;
265 fhRe2[ic*fNPID+ipid] = new TH3D(key,title,nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax,nmassbins,massmin,massmax) ;
477d6cee 266 outputContainer->Add(fhRe2[ic*fNPID+ipid]) ;
267
268 sprintf(key,"hMi_cen%d_pid%d_dist2",ic,ipid) ;
269 sprintf(title,"Mixed m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
5a2dbc3c 270 //fhMi2[ic*fNPID+ipid]=(TH3D*)fhEtalon->Clone(key) ;
271 //fhMi2[ic*fNPID+ipid]->SetName(key) ;
272 //fhMi2[ic*fNPID+ipid]->SetTitle(title) ;
273 fhMi2[ic*fNPID+ipid] = new TH3D(key,title,nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax,nmassbins,massmin,massmax) ;
477d6cee 274 outputContainer->Add(fhMi2[ic*fNPID+ipid]) ;
275
276 //Distance to bad module 3
277 sprintf(key,"hRe_cen%d_pid%d_dist3",ic,ipid) ;
278 sprintf(title,"Real m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
5a2dbc3c 279 //fhRe3[ic*fNPID+ipid]=(TH3D*)fhEtalon->Clone(key) ;
280 //fhRe3[ic*fNPID+ipid]->SetName(key) ;
281 //fhRe3[ic*fNPID+ipid]->SetTitle(title) ;
282 fhRe3[ic*fNPID+ipid] = new TH3D(key,title,nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax,nmassbins,massmin,massmax) ;
477d6cee 283 outputContainer->Add(fhRe3[ic*fNPID+ipid]) ;
284
285 sprintf(key,"hMi_cen%d_pid%d_dist3",ic,ipid) ;
286 sprintf(title,"Mixed m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
5a2dbc3c 287 //fhMi3[ic*fNPID+ipid]=(TH3D*)fhEtalon->Clone(key) ;
288 //fhMi3[ic*fNPID+ipid]->SetName(key) ;
289 //fhMi3[ic*fNPID+ipid]->SetTitle(title) ;
290 fhMi3[ic*fNPID+ipid] = new TH3D(key,title,nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax,nmassbins,massmin,massmax) ;
477d6cee 291 outputContainer->Add(fhMi3[ic*fNPID+ipid]) ;
1c5acb87 292 }
477d6cee 293 }
294
295
296 fhEvents=new TH3D("hEvents","Number of events",fNCentrBin,0.,1.*fNCentrBin,
297 fNZvertBin,0.,1.*fNZvertBin,fNrpBin,0.,1.*fNrpBin) ;
298 outputContainer->Add(fhEvents) ;
299
50f39b97 300
301 fhRealOpeningAngle = new TH2D
302 ("hRealOpeningAngle","Angle between all #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,200,0,0.5);
303 fhRealOpeningAngle->SetYTitle("#theta(rad)");
304 fhRealOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
305 outputContainer->Add(fhRealOpeningAngle) ;
306
307 fhRealCosOpeningAngle = new TH2D
308 ("hRealCosOpeningAngle","Cosinus of angle between all #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,200,-1,1);
309 fhRealCosOpeningAngle->SetYTitle("cos (#theta) ");
310 fhRealCosOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
311 outputContainer->Add(fhRealCosOpeningAngle) ;
312
313
477d6cee 314 //Histograms filled only if MC data is requested
315 if(IsDataMC() || (GetReader()->GetDataType() == AliCaloTrackReader::kMC) ){
5a2dbc3c 316 // if(fhEtalon->GetXaxis()->GetXbins() && fhEtalon->GetXaxis()->GetXbins()->GetSize()){ //Variable bin size
317 // fhPrimPt = new TH1D("hPrimPt","Primary pi0 pt",fhEtalon->GetXaxis()->GetNbins(),fhEtalon->GetXaxis()->GetXbins()->GetArray()) ;
318 // fhPrimAccPt = new TH1D("hPrimAccPt","Primary pi0 pt with both photons in acceptance",fhEtalon->GetXaxis()->GetNbins(),
319 // fhEtalon->GetXaxis()->GetXbins()->GetArray()) ;
320 // }
321 // else{
322 // fhPrimPt = new TH1D("hPrimPt","Primary pi0 pt",fhEtalon->GetXaxis()->GetNbins(),fhEtalon->GetXaxis()->GetXmin(),fhEtalon->GetXaxis()->GetXmax()) ;
323 // fhPrimAccPt = new TH1D("hPrimAccPt","Primary pi0 pt with both photons in acceptance",
324 // fhEtalon->GetXaxis()->GetNbins(),fhEtalon->GetXaxis()->GetXmin(),fhEtalon->GetXaxis()->GetXmax()) ;
325 // }
326
327 fhPrimPt = new TH1D("hPrimPt","Primary pi0 pt",nptbins,ptmin,ptmax) ;
328 fhPrimAccPt = new TH1D("hPrimAccPt","Primary pi0 pt with both photons in acceptance",nptbins,ptmin,ptmax) ;
477d6cee 329 outputContainer->Add(fhPrimPt) ;
330 outputContainer->Add(fhPrimAccPt) ;
331
5a2dbc3c 332 fhPrimY = new TH1D("hPrimaryRapidity","Rapidity of primary pi0",netabins,etamin,etamax) ;
477d6cee 333 outputContainer->Add(fhPrimY) ;
334
5a2dbc3c 335 fhPrimAccY = new TH1D("hPrimAccRapidity","Rapidity of primary pi0",netabins,etamin,etamax) ;
477d6cee 336 outputContainer->Add(fhPrimAccY) ;
337
5a2dbc3c 338 fhPrimPhi = new TH1D("hPrimaryPhi","Azimithal of primary pi0",nphibins,phimin*TMath::RadToDeg(),phimax*TMath::RadToDeg()) ;
477d6cee 339 outputContainer->Add(fhPrimPhi) ;
340
5a2dbc3c 341 fhPrimAccPhi = new TH1D("hPrimAccPhi","Azimithal of primary pi0 with accepted daughters",nphibins,phimin*TMath::RadToDeg(),phimax*TMath::RadToDeg()) ;
477d6cee 342 outputContainer->Add(fhPrimAccPhi) ;
50f39b97 343
344
345 fhPrimOpeningAngle = new TH2D
4df35693 346 ("hPrimOpeningAngle","Angle between all primary #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,100,0,0.5);
50f39b97 347 fhPrimOpeningAngle->SetYTitle("#theta(rad)");
348 fhPrimOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
349 outputContainer->Add(fhPrimOpeningAngle) ;
350
351 fhPrimCosOpeningAngle = new TH2D
4df35693 352 ("hPrimCosOpeningAngle","Cosinus of angle between all primary #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,100,-1,1);
50f39b97 353 fhPrimCosOpeningAngle->SetYTitle("cos (#theta) ");
354 fhPrimCosOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
355 outputContainer->Add(fhPrimCosOpeningAngle) ;
356
477d6cee 357 }
50f39b97 358
6921fa00 359 for(Int_t imod=0; imod<fNModules; imod++){
50f39b97 360 //Module dependent invariant mass
361 sprintf(key,"hReMod_%d",imod) ;
362 sprintf(title,"Real m_{#gamma#gamma} distr. for Module %d",imod) ;
363 //fhEtalon->Clone(key);
364 //fhReMod[imod]=(TH3D*)fhEtalon->Clone(key) ;
365 //fhReMod[imod]->SetName(key) ;
366 //fhReMod[imod]->SetTitle(title) ;
367 fhReMod[imod] = new TH3D(key,title,nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax,nmassbins,massmin,massmax) ;
368 outputContainer->Add(fhReMod[imod]) ;
6921fa00 369 }
50f39b97 370
477d6cee 371 return outputContainer;
1c5acb87 372}
373
374//_________________________________________________________________________________________________________________________________________________
375void AliAnaPi0::Print(const Option_t * /*opt*/) const
376{
477d6cee 377 //Print some relevant parameters set for the analysis
a3aebfff 378 printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
477d6cee 379 AliAnaPartCorrBaseClass::Print(" ");
a3aebfff 380
477d6cee 381 printf("Number of bins in Centrality: %d \n",fNCentrBin) ;
382 printf("Number of bins in Z vert. pos: %d \n",fNZvertBin) ;
383 printf("Number of bins in Reac. Plain: %d \n",fNrpBin) ;
384 printf("Depth of event buffer: %d \n",fNmaxMixEv) ;
385 printf("Number of different PID used: %d \n",fNPID) ;
386 printf("Cuts: \n") ;
387 printf("Z vertex position: -%2.3f < z < %2.3f \n",fZvtxCut,fZvtxCut) ;
50f39b97 388 printf("Number of modules: %d \n",fNModules) ;
389 printf("Select pairs with their angle: %d \n",fUseAngleCut) ;
477d6cee 390 printf("------------------------------------------------------\n") ;
1c5acb87 391}
392
393
394//____________________________________________________________________________________________________________________________________________________
6639984f 395void AliAnaPi0::MakeAnalysisFillHistograms()
1c5acb87 396{
477d6cee 397 //Process one event and extract photons from AOD branch
398 // filled with AliAnaPhoton and fill histos with invariant mass
399
400 //Apply some cuts on event: vertex position and centrality range
401 Int_t iRun=(GetReader()->GetInputEvent())->GetRunNumber() ;
402 if(IsBadRun(iRun)) return ;
403
404 Double_t vert[]={0,0,0} ; //vertex ;
405 GetReader()->GetVertex(vert);
406 if(vert[2]<-fZvtxCut || vert[2]> fZvtxCut) return ; //Event can not be used (vertex, centrality,... cuts not fulfilled)
407
408 //Get Centrality and calculate centrality bin
409 //Does not exist in ESD yet???????
410 Int_t curCentrBin=0 ;
411
412 //Get Reaction Plain position and calculate RP bin
413 //does not exist in ESD yet????
414 Int_t curRPBin=0 ;
415
416 Int_t curZvertBin=(Int_t)(0.5*fNZvertBin*(vert[2]+fZvtxCut)/fZvtxCut) ;
417
418 fhEvents->Fill(curCentrBin+0.5,curZvertBin+0.5,curRPBin+0.5) ;
419
420 Int_t nPhot = GetInputAODBranch()->GetEntriesFast() ;
a3aebfff 421 if(GetDebug() > 1) printf("AliAnaPi0::MakeAnalysisFillHistograms() - Photon entries %d\n", nPhot);
6921fa00 422 Int_t module1 = -1;
423 Int_t module2 = -1;
477d6cee 424 for(Int_t i1=0; i1<nPhot-1; i1++){
425 AliAODPWG4Particle * p1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i1)) ;
426 TLorentzVector photon1(p1->Px(),p1->Py(),p1->Pz(),p1->E());
59b6bd99 427 //Get Module number
428 module1 = GetModuleNumber(p1);
477d6cee 429 for(Int_t i2=i1+1; i2<nPhot; i2++){
430 AliAODPWG4Particle * p2 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i2)) ;
431 TLorentzVector photon2(p2->Px(),p2->Py(),p2->Pz(),p2->E());
59b6bd99 432 //Get module number
433 module2 = GetModuleNumber(p2);
477d6cee 434 Double_t m = (photon1 + photon2).M() ;
435 Double_t pt = (photon1 + photon2).Pt();
436 Double_t a = TMath::Abs(p1->E()-p2->E())/(p1->E()+p2->E()) ;
437 if(GetDebug() > 2)
a3aebfff 438 printf("AliAnaPi0::MakeAnalysisFillHistograms() - Current Event: pT: photon1 %2.2f, photon2 %2.2f; Pair: pT %2.2f, mass %2.3f, a %f2.3\n",
477d6cee 439 p1->Pt(), p2->Pt(), pt,m,a);
50f39b97 440
441 //Check if opening angle is too large or too small compared to what is expected
442 Double_t angle = photon1.Angle(photon2.Vect());
443 //if(fUseAngleCut && !GetNeutralMesonSelection()->IsAngleInWindow((photon1+photon2).E(),angle)) continue;
444 //printf("angle %f\n",angle);
445 if(fUseAngleCut && angle < 0.1) continue;
446 fhRealOpeningAngle ->Fill(pt,angle);
447 fhRealCosOpeningAngle->Fill(pt,TMath::Cos(angle));
448
59b6bd99 449 //Fill module dependent histograms
450 //if(module1==module2) printf("mod1 %d\n",module1);
451 if(module1==module2 && module1 >=0 && module1<fNModules)
452 fhReMod[module1]->Fill(pt,a,m) ;
453
477d6cee 454 for(Int_t ipid=0; ipid<fNPID; ipid++)
455 {
456 if((p1->IsPIDOK(ipid,AliCaloPID::kPhoton)) && (p2->IsPIDOK(ipid,AliCaloPID::kPhoton))){
457 fhRe1[curCentrBin*fNPID+ipid]->Fill(pt,a,m) ;
458 if(p1->DistToBad()>0 && p2->DistToBad()>0){
459 fhRe2[curCentrBin*fNPID+ipid]->Fill(pt,a,m) ;
460 if(p1->DistToBad()>1 && p2->DistToBad()>1){
461 fhRe3[curCentrBin*fNPID+ipid]->Fill(pt,a,m) ;
462 }
463 }
464 }
465 }
466 }
467 }
468
469 //Fill mixed
470
471 TList * evMixList=fEventsList[curCentrBin*fNZvertBin*fNrpBin+curZvertBin*fNrpBin+curRPBin] ;
472 Int_t nMixed = evMixList->GetSize() ;
473 for(Int_t ii=0; ii<nMixed; ii++){
474 TClonesArray* ev2= (TClonesArray*) (evMixList->At(ii));
475 Int_t nPhot2=ev2->GetEntriesFast() ;
476 Double_t m = -999;
a3aebfff 477 if(GetDebug() > 1) printf("AliAnaPi0::MakeAnalysisFillHistograms() - Mixed event %d photon entries %d\n", ii, nPhot);
477d6cee 478
479 for(Int_t i1=0; i1<nPhot; i1++){
480 AliAODPWG4Particle * p1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i1)) ;
481 TLorentzVector photon1(p1->Px(),p1->Py(),p1->Pz(),p1->E());
482 for(Int_t i2=0; i2<nPhot2; i2++){
483 AliAODPWG4Particle * p2 = (AliAODPWG4Particle*) (ev2->At(i2)) ;
1c5acb87 484
477d6cee 485 TLorentzVector photon2(p2->Px(),p2->Py(),p2->Pz(),p2->E());
486 m = (photon1+photon2).M() ;
487 Double_t pt = (photon1 + photon2).Pt();
488 Double_t a = TMath::Abs(p1->E()-p2->E())/(p1->E()+p2->E()) ;
50f39b97 489
490 //Check if opening angle is too large or too small compared to what is expected
491 Double_t angle = photon1.Angle(photon2.Vect());
492 //if(fUseAngleCut && !GetNeutralMesonSelection()->IsAngleInWindow((photon1+photon2).E(),angle)) continue;
493 if(fUseAngleCut && angle < 0.1) continue;
494
477d6cee 495 if(GetDebug() > 2)
a3aebfff 496 printf("AliAnaPi0::MakeAnalysisFillHistograms() - Mixed Event: pT: photon1 %2.2f, photon2 %2.2f; Pair: pT %2.2f, mass %2.3f, a %f2.3\n",
477d6cee 497 p1->Pt(), p2->Pt(), pt,m,a);
498 for(Int_t ipid=0; ipid<fNPID; ipid++){
499 if((p1->IsPIDOK(ipid,AliCaloPID::kPhoton)) && (p2->IsPIDOK(ipid,AliCaloPID::kPhoton))){
500 fhMi1[curCentrBin*fNPID+ipid]->Fill(pt,a,m) ;
501 if(p1->DistToBad()>0 && p2->DistToBad()>0){
502 fhMi2[curCentrBin*fNPID+ipid]->Fill(pt,a,m) ;
503 if(p1->DistToBad()>1 && p2->DistToBad()>1){
504 fhMi3[curCentrBin*fNPID+ipid]->Fill(pt,a,m) ;
505 }
506
507 }
508 }
1c5acb87 509 }
477d6cee 510 }
511 }
512 }
513
514 TClonesArray *currentEvent = new TClonesArray(*GetInputAODBranch());
515 //Add current event to buffer and Remove redandant events
516 if(currentEvent->GetEntriesFast()>0){
517 evMixList->AddFirst(currentEvent) ;
518 currentEvent=0 ; //Now list of particles belongs to buffer and it will be deleted with buffer
519 if(evMixList->GetSize()>=fNmaxMixEv)
520 {
521 TClonesArray * tmp = (TClonesArray*) (evMixList->Last()) ;
522 evMixList->RemoveLast() ;
523 delete tmp ;
524 }
525 }
526 else{ //empty event
527 delete currentEvent ;
528 currentEvent=0 ;
529 }
530
531 //Acceptance
de8a210f 532 if(IsDataMC() && GetReader()->ReadStack()){
50f39b97 533 AliStack * stack = GetMCStack();
534 if(stack && (IsDataMC() || (GetReader()->GetDataType() == AliCaloTrackReader::kMC)) ){
535 for(Int_t i=0 ; i<stack->GetNprimary(); i++){
536 TParticle * prim = stack->Particle(i) ;
537 if(prim->GetPdgCode() == 111){
538 Double_t pi0Pt = prim->Pt() ;
539 //printf("pi0, pt %2.2f\n",pi0Pt);
540 if(prim->Energy() == TMath::Abs(prim->Pz())) continue ; //Protection against floating point exception
541 Double_t pi0Y = 0.5*TMath::Log((prim->Energy()-prim->Pz())/(prim->Energy()+prim->Pz())) ;
542 Double_t phi = TMath::RadToDeg()*prim->Phi() ;
543 if(TMath::Abs(pi0Y) < 0.5){
544 fhPrimPt->Fill(pi0Pt) ;
545 }
546 fhPrimY ->Fill(pi0Y) ;
547 fhPrimPhi->Fill(phi) ;
548
549 //Check if both photons hit Calorimeter
550 Int_t iphot1=prim->GetFirstDaughter() ;
551 Int_t iphot2=prim->GetLastDaughter() ;
552 if(iphot1>-1 && iphot1<stack->GetNtrack() && iphot2>-1 && iphot2<stack->GetNtrack()){
553 TParticle * phot1 = stack->Particle(iphot1) ;
554 TParticle * phot2 = stack->Particle(iphot2) ;
555 if(phot1 && phot2 && phot1->GetPdgCode()==22 && phot2->GetPdgCode()==22){
556 //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",
557 // phot1->Pt(), phot1->Phi()*180./3.1415, phot1->Eta(), phot2->Pt(), phot2->Phi()*180./3.1415, phot2->Eta());
558
559 TLorentzVector lv1, lv2;
560 phot1->Momentum(lv1);
561 phot2->Momentum(lv2);
562
563 Bool_t inacceptance = kFALSE;
564 if(fCalorimeter == "PHOS"){
765d44e7 565 if(GetPHOSGeometry() && GetCaloUtils()->IsPHOSGeoMatrixSet()){
50f39b97 566 Int_t mod ;
567 Double_t x,z ;
765d44e7 568 if(GetPHOSGeometry()->ImpactOnEmc(phot1,mod,z,x) && GetPHOSGeometry()->ImpactOnEmc(phot2,mod,z,x))
50f39b97 569 inacceptance = kTRUE;
570 if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
571 }
572 else{
573
574 if(GetFiducialCut()->IsInFiducialCut(lv1,fCalorimeter) && GetFiducialCut()->IsInFiducialCut(lv2,fCalorimeter))
575 inacceptance = kTRUE ;
576 if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
577 }
591cc579 578
50f39b97 579 }
765d44e7 580 else if(fCalorimeter == "EMCAL" && GetCaloUtils()->IsEMCALGeoMatrixSet()){
581 if(GetEMCALGeometry()){
582 if(GetEMCALGeometry()->Impact(phot1) && GetEMCALGeometry()->Impact(phot2))
50f39b97 583 inacceptance = kTRUE;
584 if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
585 }
586 else{
587 if(GetFiducialCut()->IsInFiducialCut(lv1,fCalorimeter) && GetFiducialCut()->IsInFiducialCut(lv2,fCalorimeter))
588 inacceptance = kTRUE ;
589 if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
590 }
591 }
592
593 if(inacceptance){
594
595 fhPrimAccPt->Fill(pi0Pt) ;
596 fhPrimAccPhi->Fill(phi) ;
597 fhPrimAccY->Fill(pi0Y) ;
598 Double_t angle = lv1.Angle(lv2.Vect());
599 fhPrimOpeningAngle ->Fill(pi0Pt,angle);
600 fhPrimCosOpeningAngle->Fill(pi0Pt,TMath::Cos(angle));
601
602 }//Accepted
603 }// 2 photons
604 }//Check daughters exist
605 }// Primary pi0
606 }//loop on primaries
607 }//stack exists and data is MC
591cc579 608 }//read stack
609 else if(GetReader()->ReadAODMCParticles()){
50f39b97 610 if(GetDebug() >= 0) printf("AliAnaPi0::MakeAnalysisFillHistograms() - Acceptance calculation with MCParticles not implemented yet\n");
591cc579 611 }
477d6cee 612
1c5acb87 613}
614
a5cc4f03 615//________________________________________________________________________
616void AliAnaPi0::ReadHistograms(TList* outputList)
617{
50f39b97 618 // Needed when Terminate is executed in distributed environment
619 // Refill analysis histograms of this class with corresponding histograms in output list.
620
621 // Histograms of this analsys are kept in the same list as other analysis, recover the position of
622 // the first one and then add the next.
623 Int_t index = outputList->IndexOf(outputList->FindObject(GetAddedHistogramsStringToName()+"hRe_cen0_pid0_dist1"));
624
625 if(!fhRe1) fhRe1 = new TH3D*[fNCentrBin*fNPID] ;
626 if(!fhRe2) fhRe2 = new TH3D*[fNCentrBin*fNPID] ;
627 if(!fhRe3) fhRe3 = new TH3D*[fNCentrBin*fNPID] ;
628 if(!fhMi1) fhMi1 = new TH3D*[fNCentrBin*fNPID] ;
629 if(!fhMi2) fhMi2 = new TH3D*[fNCentrBin*fNPID] ;
630 if(!fhMi3) fhMi3 = new TH3D*[fNCentrBin*fNPID] ;
631 if(!fhReMod) fhReMod = new TH3D*[fNModules] ;
632
633 for(Int_t ic=0; ic<fNCentrBin; ic++){
634 for(Int_t ipid=0; ipid<fNPID; ipid++){
635 fhRe1[ic*fNPID+ipid] = (TH3D*) outputList->At(index++);
636 fhMi1[ic*fNPID+ipid] = (TH3D*) outputList->At(index++);
637 fhRe2[ic*fNPID+ipid] = (TH3D*) outputList->At(index++);
638 fhMi2[ic*fNPID+ipid] = (TH3D*) outputList->At(index++);
639 fhRe3[ic*fNPID+ipid] = (TH3D*) outputList->At(index++);
640 fhMi3[ic*fNPID+ipid] = (TH3D*) outputList->At(index++);
641 }
642 }
643
644 fhEvents = (TH3D *) outputList->At(index++);
645
646 //Histograms filled only if MC data is requested
647 if(IsDataMC() || (GetReader()->GetDataType() == AliCaloTrackReader::kMC) ){
648 fhPrimPt = (TH1D*) outputList->At(index++);
649 fhPrimAccPt = (TH1D*) outputList->At(index++);
650 fhPrimY = (TH1D*) outputList->At(index++);
651 fhPrimAccY = (TH1D*) outputList->At(index++);
652 fhPrimPhi = (TH1D*) outputList->At(index++);
653 fhPrimAccPhi = (TH1D*) outputList->At(index++);
654 }
655
656 for(Int_t imod=0; imod < fNModules; imod++)
657 fhReMod[imod] = (TH3D*) outputList->At(index++);
658
a5cc4f03 659}
660
661
6639984f 662//____________________________________________________________________________________________________________________________________________________
a5cc4f03 663void AliAnaPi0::Terminate(TList* outputList)
6639984f 664{
665 //Do some calculations and plots from the final histograms.
477d6cee 666
fbeaf916 667 printf(" *** %s Terminate:\n", GetName()) ;
50f39b97 668
a5cc4f03 669 //Recover histograms from output histograms list, needed for distributed analysis.
670 ReadHistograms(outputList);
50f39b97 671
2e557d1c 672 if (!fhRe1) {
50f39b97 673 printf("AliAnaPi0::Terminate() - Error: Remote output histograms not imported in AliAnaPi0 object");
674 return;
2e557d1c 675 }
50f39b97 676
a3aebfff 677 printf("AliAnaPi0::Terminate() Mgg Real : %5.3f , RMS : %5.3f \n", fhRe1[0]->GetMean(), fhRe1[0]->GetRMS() ) ;
50f39b97 678
71dd883b 679 char nameIM[128];
680 sprintf(nameIM,"AliAnaPi0_%s_cPt",fCalorimeter.Data());
681 TCanvas * cIM = new TCanvas(nameIM, "", 400, 10, 600, 700) ;
6639984f 682 cIM->Divide(2, 2);
50f39b97 683
6639984f 684 cIM->cd(1) ;
685 //gPad->SetLogy();
583c48f1 686 TH1D * hIMAllPt = (TH1D*) fhRe1[0]->ProjectionZ(Form("IMPtAll_%s",fCalorimeter.Data()));
6639984f 687 hIMAllPt->SetLineColor(2);
688 hIMAllPt->SetTitle("No cut on p_{T, #gamma#gamma} ");
689 hIMAllPt->Draw();
690
691 cIM->cd(2) ;
583c48f1 692 TH3F * hRe1Pt5 = (TH3F*)fhRe1[0]->Clone(Form("IMPt5_%s",fCalorimeter.Data()));
6639984f 693 hRe1Pt5->GetXaxis()->SetRangeUser(0,5);
583c48f1 694 TH1D * hIMPt5 = (TH1D*) hRe1Pt5->Project3D(Form("IMPt5_%s_pz",fCalorimeter.Data()));
6639984f 695 hIMPt5->SetLineColor(2);
696 hIMPt5->SetTitle("0 < p_{T, #gamma#gamma} < 5 GeV/c");
697 hIMPt5->Draw();
698
699 cIM->cd(3) ;
583c48f1 700 TH3F * hRe1Pt10 = (TH3F*)fhRe1[0]->Clone(Form("IMPt10_%s",fCalorimeter.Data()));
6639984f 701 hRe1Pt10->GetXaxis()->SetRangeUser(5,10);
583c48f1 702 TH1D * hIMPt10 = (TH1D*) hRe1Pt10->Project3D(Form("IMPt10_%s_pz",fCalorimeter.Data()));
6639984f 703 hIMPt10->SetLineColor(2);
704 hIMPt10->SetTitle("5 < p_{T, #gamma#gamma} < 10 GeV/c");
705 hIMPt10->Draw();
706
707 cIM->cd(4) ;
583c48f1 708 TH3F * hRe1Pt20 = (TH3F*)fhRe1[0]->Clone(Form("IMPt20_%s",fCalorimeter.Data()));
6639984f 709 hRe1Pt20->GetXaxis()->SetRangeUser(10,20);
583c48f1 710 TH1D * hIMPt20 = (TH1D*) hRe1Pt20->Project3D(Form("IMPt20_%s_pz",fCalorimeter.Data()));
6639984f 711 hIMPt20->SetLineColor(2);
712 hIMPt20->SetTitle("10 < p_{T, #gamma#gamma} < 20 GeV/c");
713 hIMPt20->Draw();
714
71dd883b 715 char nameIMF[128];
716 sprintf(nameIMF,"AliAnaPi0_%s_Mgg.eps",fCalorimeter.Data());
717 cIM->Print(nameIMF);
6639984f 718
71dd883b 719 char namePt[128];
720 sprintf(namePt,"AliAnaPi0_%s_cPt",fCalorimeter.Data());
721 TCanvas * cPt = new TCanvas(namePt, "", 400, 10, 600, 700) ;
6639984f 722 cPt->Divide(2, 2);
723
724 cPt->cd(1) ;
725 //gPad->SetLogy();
726 TH1D * hPt = (TH1D*) fhRe1[0]->Project3D("x");
727 hPt->SetLineColor(2);
728 hPt->SetTitle("No cut on M_{#gamma#gamma} ");
729 hPt->Draw();
730
731 cPt->cd(2) ;
583c48f1 732 TH3F * hRe1IM1 = (TH3F*)fhRe1[0]->Clone(Form("Pt1_%s",fCalorimeter.Data()));
6639984f 733 hRe1IM1->GetZaxis()->SetRangeUser(0.05,0.21);
734 TH1D * hPtIM1 = (TH1D*) hRe1IM1->Project3D("x");
735 hPtIM1->SetLineColor(2);
736 hPtIM1->SetTitle("0.05 < M_{#gamma#gamma} < 0.21 GeV/c^{2}");
737 hPtIM1->Draw();
738
739 cPt->cd(3) ;
583c48f1 740 TH3F * hRe1IM2 = (TH3F*)fhRe1[0]->Clone(Form("Pt2_%s",fCalorimeter.Data()));
6639984f 741 hRe1IM2->GetZaxis()->SetRangeUser(0.09,0.17);
742 TH1D * hPtIM2 = (TH1D*) hRe1IM2->Project3D("x");
743 hPtIM2->SetLineColor(2);
744 hPtIM2->SetTitle("0.09 < M_{#gamma#gamma} < 0.17 GeV/c^{2}");
745 hPtIM2->Draw();
746
747 cPt->cd(4) ;
583c48f1 748 TH3F * hRe1IM3 = (TH3F*)fhRe1[0]->Clone(Form("Pt3_%s",fCalorimeter.Data()));
6639984f 749 hRe1IM3->GetZaxis()->SetRangeUser(0.11,0.15);
750 TH1D * hPtIM3 = (TH1D*) hRe1IM1->Project3D("x");
751 hPtIM3->SetLineColor(2);
752 hPtIM3->SetTitle("0.11 < M_{#gamma#gamma} < 0.15 GeV/c^{2}");
753 hPtIM3->Draw();
754
71dd883b 755 char namePtF[128];
756 sprintf(namePtF,"AliAnaPi0_%s_Pt.eps",fCalorimeter.Data());
757 cPt->Print(namePtF);
1c5acb87 758
6639984f 759
760 char line[1024] ;
71dd883b 761 sprintf(line, ".!tar -zcf %s_%s.tar.gz *.eps", GetName(),fCalorimeter.Data()) ;
6639984f 762 gROOT->ProcessLine(line);
71dd883b 763 sprintf(line, ".!rm -fR AliAnaPi0_%s*.eps",fCalorimeter.Data());
6639984f 764 gROOT->ProcessLine(line);
765
71dd883b 766 printf(" AliAnaPi0::Terminate() - !! All the eps files are in %s_%s.tar.gz !!!\n", GetName(), fCalorimeter.Data());
1c5acb87 767
6639984f 768}
1c5acb87 769
770
771
772
773