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