1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
17 //_________________________________________________________________________
18 // Class to collect two-photon invariant mass distributions for
19 // extractin raw pi0 yield.
21 //-- Author: Dmitri Peressounko (RRC "KI")
22 //-- Adapted to PartCorr frame by Lamia Benhabib (SUBATECH)
23 //-- and Gustavo Conesa (INFN-Frascati)
24 //_________________________________________________________________________
27 // --- ROOT system ---
30 //#include "Riostream.h"
34 #include "TClonesArray.h"
35 #include "TObjString.h"
36 #include "TGeoManager.h"
38 //---- AliRoot system ----
39 #include "AliAnaPi0.h"
40 #include "AliCaloTrackReader.h"
41 #include "AliCaloPID.h"
43 #include "AliFiducialCut.h"
44 #include "TParticle.h"
45 #include "AliAODCaloCluster.h"
46 #include "AliVEvent.h"
47 #include "AliESDCaloCluster.h"
48 #include "AliESDEvent.h"
49 #include "AliAODEvent.h"
50 #include "AliNeutralMesonSelection.h"
54 //________________________________________________________________________________________________________________________________________________
55 AliAnaPi0::AliAnaPi0() : AliAnaPartCorrBaseClass(),
56 fNCentrBin(0),fNZvertBin(0),fNrpBin(0),
57 fNPID(0),fNmaxMixEv(0), fZvtxCut(0.),fCalorimeter(""),
58 fNModules(12), fUseAngleCut(kFALSE), fEventsList(0x0), //fhEtalon(0x0),
59 fhReMod(0x0), fhRe1(0x0),fhMi1(0x0),fhRe2(0x0),fhMi2(0x0),fhRe3(0x0),fhMi3(0x0),fhEvents(0x0),
60 fhRealOpeningAngle(0x0),fhRealCosOpeningAngle(0x0),
61 fhPrimPt(0x0), fhPrimAccPt(0x0), fhPrimY(0x0), fhPrimAccY(0x0), fhPrimPhi(0x0), fhPrimAccPhi(0x0),
62 fhPrimOpeningAngle(0x0),fhPrimCosOpeningAngle(0x0)
69 //________________________________________________________________________________________________________________________________________________
70 AliAnaPi0::AliAnaPi0(const AliAnaPi0 & ex) : AliAnaPartCorrBaseClass(ex),
71 fNCentrBin(ex.fNCentrBin),fNZvertBin(ex.fNZvertBin),fNrpBin(ex.fNrpBin),
72 fNPID(ex.fNPID),fNmaxMixEv(ex.fNmaxMixEv),fZvtxCut(ex.fZvtxCut), fCalorimeter(ex.fCalorimeter),
73 fNModules(ex.fNModules), fUseAngleCut(ex.fUseAngleCut), fEventsList(ex.fEventsList), //fhEtalon(ex.fhEtalon),
74 fhReMod(ex.fhReMod), fhRe1(ex.fhRe1),fhMi1(ex.fhMi1),fhRe2(ex.fhRe2),fhMi2(ex.fhMi2),
75 fhRe3(ex.fhRe3),fhMi3(ex.fhMi3),fhEvents(ex.fhEvents),
76 fhRealOpeningAngle(ex.fhRealOpeningAngle),fhRealCosOpeningAngle(ex.fhRealCosOpeningAngle),
77 fhPrimPt(ex.fhPrimPt), fhPrimAccPt(ex.fhPrimAccPt), fhPrimY(ex.fhPrimY),
78 fhPrimAccY(ex.fhPrimAccY), fhPrimPhi(ex.fhPrimPhi), fhPrimAccPhi(ex.fhPrimAccPhi),
79 fhPrimOpeningAngle(ex.fhPrimOpeningAngle),fhPrimCosOpeningAngle(ex.fhPrimCosOpeningAngle)
85 //________________________________________________________________________________________________________________________________________________
86 AliAnaPi0 & AliAnaPi0::operator = (const AliAnaPi0 & ex)
88 // assignment operator
90 if(this == &ex)return *this;
91 ((AliAnaPartCorrBaseClass *)this)->operator=(ex);
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;
107 //________________________________________________________________________________________________________________________________________________
108 AliAnaPi0::~AliAnaPi0() {
109 // Remove event containers
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] ;
119 delete[] fEventsList;
125 //________________________________________________________________________________________________________________________________________________
126 void AliAnaPi0::InitParameters()
128 //Init parameters when first called the analysis
129 //Set default parameters
130 SetInputAODName("PWG4Particle");
132 AddToHistogramsName("AnaPi0_");
133 fNModules = 12; // set maximum to maximum number of EMCAL modules
140 fCalorimeter = "PHOS";
141 fUseAngleCut = kFALSE;
144 //________________________________________________________________________________________________________________________________________________
145 //void AliAnaPi0::Init()
147 //Init some data members needed in analysis
149 //Histograms binning and range
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)") ;
159 //________________________________________________________________________________________________________________________________________________
160 TList * AliAnaPi0::GetCreateOutputObjects()
162 // Create histograms to be saved in output file and
163 // store them in fOutputContainer
165 //create event containers
166 fEventsList = new TList*[fNCentrBin*fNZvertBin*fNrpBin] ;
168 for(Int_t ic=0; ic<fNCentrBin; ic++){
169 for(Int_t iz=0; iz<fNZvertBin; iz++){
170 for(Int_t irp=0; irp<fNrpBin; irp++){
171 fEventsList[ic*fNZvertBin*fNrpBin+iz*fNrpBin+irp] = new TList() ;
176 //If Geometry library loaded, do geometry selection during analysis.
177 if(fCalorimeter=="PHOS"){
178 if(!GetReader()->GetPHOSGeometry()) printf("AliAnaPi0::GetCreateOutputObjects() - Initialize PHOS geometry!\n");
179 GetReader()->InitPHOSGeometry();
182 else if(fCalorimeter=="EMCAL"){
183 if(!GetReader()->GetEMCALGeometry()) printf("AliAnaPi0::GetCreateOutputObjects() - Initialize EMCAL geometry!\n");
184 GetReader()->InitEMCALGeometry();
187 TList * outputContainer = new TList() ;
188 outputContainer->SetName(GetName());
190 fhReMod = new TH3D*[fNModules] ;
191 fhRe1 = new TH3D*[fNCentrBin*fNPID] ;
192 fhRe2 = new TH3D*[fNCentrBin*fNPID] ;
193 fhRe3 = new TH3D*[fNCentrBin*fNPID] ;
194 fhMi1 = new TH3D*[fNCentrBin*fNPID] ;
195 fhMi2 = new TH3D*[fNCentrBin*fNPID] ;
196 fhMi3 = new TH3D*[fNCentrBin*fNPID] ;
200 Int_t nptbins = GetHistoPtBins();
201 Int_t nphibins = GetHistoPhiBins();
202 Int_t netabins = GetHistoEtaBins();
203 Float_t ptmax = GetHistoPtMax();
204 Float_t phimax = GetHistoPhiMax();
205 Float_t etamax = GetHistoEtaMax();
206 Float_t ptmin = GetHistoPtMin();
207 Float_t phimin = GetHistoPhiMin();
208 Float_t etamin = GetHistoEtaMin();
210 Int_t nmassbins = GetHistoMassBins();
211 Int_t nasymbins = GetHistoAsymmetryBins();
212 Float_t massmax = GetHistoMassMax();
213 Float_t asymmax = GetHistoAsymmetryMax();
214 Float_t massmin = GetHistoMassMin();
215 Float_t asymmin = GetHistoAsymmetryMin();
217 for(Int_t ic=0; ic<fNCentrBin; ic++){
218 for(Int_t ipid=0; ipid<fNPID; ipid++){
220 //Distance to bad module 1
221 sprintf(key,"hRe_cen%d_pid%d_dist1",ic,ipid) ;
222 sprintf(title,"Real m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
224 //fhEtalon->Clone(key);
225 //fhRe1[ic*fNPID+ipid]=(TH3D*)fhEtalon->Clone(key) ;
226 //fhRe1[ic*fNPID+ipid]->SetName(key) ;
227 //fhRe1[ic*fNPID+ipid]->SetTitle(title) ;
228 fhRe1[ic*fNPID+ipid] = new TH3D(key,title,nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax,nmassbins,massmin,massmax) ;
229 outputContainer->Add(fhRe1[ic*fNPID+ipid]) ;
231 sprintf(key,"hMi_cen%d_pid%d_dist1",ic,ipid) ;
232 sprintf(title,"Mixed m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
233 //fhMi1[ic*fNPID+ipid]=(TH3D*)fhEtalon->Clone(key) ;
234 //fhMi1[ic*fNPID+ipid]->SetName(key) ;
235 //fhMi1[ic*fNPID+ipid]->SetTitle(title) ;
236 fhMi1[ic*fNPID+ipid] = new TH3D(key,title,nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax,nmassbins,massmin,massmax) ;
237 outputContainer->Add(fhMi1[ic*fNPID+ipid]) ;
239 //Distance to bad module 2
240 sprintf(key,"hRe_cen%d_pid%d_dist2",ic,ipid) ;
241 sprintf(title,"Real m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
242 //fhRe2[ic*fNPID+ipid]=(TH3D*)fhEtalon->Clone(key) ;
243 //fhRe2[ic*fNPID+ipid]->SetName(key) ;
244 //fhRe2[ic*fNPID+ipid]->SetTitle(title) ;
245 fhRe2[ic*fNPID+ipid] = new TH3D(key,title,nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax,nmassbins,massmin,massmax) ;
246 outputContainer->Add(fhRe2[ic*fNPID+ipid]) ;
248 sprintf(key,"hMi_cen%d_pid%d_dist2",ic,ipid) ;
249 sprintf(title,"Mixed m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
250 //fhMi2[ic*fNPID+ipid]=(TH3D*)fhEtalon->Clone(key) ;
251 //fhMi2[ic*fNPID+ipid]->SetName(key) ;
252 //fhMi2[ic*fNPID+ipid]->SetTitle(title) ;
253 fhMi2[ic*fNPID+ipid] = new TH3D(key,title,nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax,nmassbins,massmin,massmax) ;
254 outputContainer->Add(fhMi2[ic*fNPID+ipid]) ;
256 //Distance to bad module 3
257 sprintf(key,"hRe_cen%d_pid%d_dist3",ic,ipid) ;
258 sprintf(title,"Real m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
259 //fhRe3[ic*fNPID+ipid]=(TH3D*)fhEtalon->Clone(key) ;
260 //fhRe3[ic*fNPID+ipid]->SetName(key) ;
261 //fhRe3[ic*fNPID+ipid]->SetTitle(title) ;
262 fhRe3[ic*fNPID+ipid] = new TH3D(key,title,nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax,nmassbins,massmin,massmax) ;
263 outputContainer->Add(fhRe3[ic*fNPID+ipid]) ;
265 sprintf(key,"hMi_cen%d_pid%d_dist3",ic,ipid) ;
266 sprintf(title,"Mixed m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
267 //fhMi3[ic*fNPID+ipid]=(TH3D*)fhEtalon->Clone(key) ;
268 //fhMi3[ic*fNPID+ipid]->SetName(key) ;
269 //fhMi3[ic*fNPID+ipid]->SetTitle(title) ;
270 fhMi3[ic*fNPID+ipid] = new TH3D(key,title,nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax,nmassbins,massmin,massmax) ;
271 outputContainer->Add(fhMi3[ic*fNPID+ipid]) ;
276 fhEvents=new TH3D("hEvents","Number of events",fNCentrBin,0.,1.*fNCentrBin,
277 fNZvertBin,0.,1.*fNZvertBin,fNrpBin,0.,1.*fNrpBin) ;
278 outputContainer->Add(fhEvents) ;
281 fhRealOpeningAngle = new TH2D
282 ("hRealOpeningAngle","Angle between all #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,200,0,0.5);
283 fhRealOpeningAngle->SetYTitle("#theta(rad)");
284 fhRealOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
285 outputContainer->Add(fhRealOpeningAngle) ;
287 fhRealCosOpeningAngle = new TH2D
288 ("hRealCosOpeningAngle","Cosinus of angle between all #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,200,-1,1);
289 fhRealCosOpeningAngle->SetYTitle("cos (#theta) ");
290 fhRealCosOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
291 outputContainer->Add(fhRealCosOpeningAngle) ;
294 //Histograms filled only if MC data is requested
295 if(IsDataMC() || (GetReader()->GetDataType() == AliCaloTrackReader::kMC) ){
296 // if(fhEtalon->GetXaxis()->GetXbins() && fhEtalon->GetXaxis()->GetXbins()->GetSize()){ //Variable bin size
297 // fhPrimPt = new TH1D("hPrimPt","Primary pi0 pt",fhEtalon->GetXaxis()->GetNbins(),fhEtalon->GetXaxis()->GetXbins()->GetArray()) ;
298 // fhPrimAccPt = new TH1D("hPrimAccPt","Primary pi0 pt with both photons in acceptance",fhEtalon->GetXaxis()->GetNbins(),
299 // fhEtalon->GetXaxis()->GetXbins()->GetArray()) ;
302 // fhPrimPt = new TH1D("hPrimPt","Primary pi0 pt",fhEtalon->GetXaxis()->GetNbins(),fhEtalon->GetXaxis()->GetXmin(),fhEtalon->GetXaxis()->GetXmax()) ;
303 // fhPrimAccPt = new TH1D("hPrimAccPt","Primary pi0 pt with both photons in acceptance",
304 // fhEtalon->GetXaxis()->GetNbins(),fhEtalon->GetXaxis()->GetXmin(),fhEtalon->GetXaxis()->GetXmax()) ;
307 fhPrimPt = new TH1D("hPrimPt","Primary pi0 pt",nptbins,ptmin,ptmax) ;
308 fhPrimAccPt = new TH1D("hPrimAccPt","Primary pi0 pt with both photons in acceptance",nptbins,ptmin,ptmax) ;
309 outputContainer->Add(fhPrimPt) ;
310 outputContainer->Add(fhPrimAccPt) ;
312 fhPrimY = new TH1D("hPrimaryRapidity","Rapidity of primary pi0",netabins,etamin,etamax) ;
313 outputContainer->Add(fhPrimY) ;
315 fhPrimAccY = new TH1D("hPrimAccRapidity","Rapidity of primary pi0",netabins,etamin,etamax) ;
316 outputContainer->Add(fhPrimAccY) ;
318 fhPrimPhi = new TH1D("hPrimaryPhi","Azimithal of primary pi0",nphibins,phimin*TMath::RadToDeg(),phimax*TMath::RadToDeg()) ;
319 outputContainer->Add(fhPrimPhi) ;
321 fhPrimAccPhi = new TH1D("hPrimAccPhi","Azimithal of primary pi0 with accepted daughters",nphibins,phimin*TMath::RadToDeg(),phimax*TMath::RadToDeg()) ;
322 outputContainer->Add(fhPrimAccPhi) ;
325 fhPrimOpeningAngle = new TH2D
326 ("hPrimOpeningAngle","Angle between all primary #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,200,0,0.5);
327 fhPrimOpeningAngle->SetYTitle("#theta(rad)");
328 fhPrimOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
329 outputContainer->Add(fhPrimOpeningAngle) ;
331 fhPrimCosOpeningAngle = new TH2D
332 ("hPrimCosOpeningAngle","Cosinus of angle between all primary #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,200,-1,1);
333 fhPrimCosOpeningAngle->SetYTitle("cos (#theta) ");
334 fhPrimCosOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
335 outputContainer->Add(fhPrimCosOpeningAngle) ;
339 for(Int_t imod=0; imod<fNModules; imod++){
340 //Module dependent invariant mass
341 sprintf(key,"hReMod_%d",imod) ;
342 sprintf(title,"Real m_{#gamma#gamma} distr. for Module %d",imod) ;
343 //fhEtalon->Clone(key);
344 //fhReMod[imod]=(TH3D*)fhEtalon->Clone(key) ;
345 //fhReMod[imod]->SetName(key) ;
346 //fhReMod[imod]->SetTitle(title) ;
347 fhReMod[imod] = new TH3D(key,title,nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax,nmassbins,massmin,massmax) ;
348 outputContainer->Add(fhReMod[imod]) ;
352 //Save parameters used for analysis
353 TString parList ; //this will be list of parameters used for this analysis.
355 sprintf(onePar,"--- AliAnaPi0 ---\n") ;
357 sprintf(onePar,"Number of bins in Centrality: %d \n",fNCentrBin) ;
359 sprintf(onePar,"Number of bins in Z vert. pos: %d \n",fNZvertBin) ;
361 sprintf(onePar,"Number of bins in Reac. Plain: %d \n",fNrpBin) ;
363 sprintf(onePar,"Depth of event buffer: %d \n",fNmaxMixEv) ;
365 sprintf(onePar,"Number of different PID used: %d \n",fNPID) ;
367 sprintf(onePar,"Cuts: \n") ;
369 sprintf(onePar,"Z vertex position: -%f < z < %f \n",fZvtxCut,fZvtxCut) ;
371 sprintf(onePar,"Calorimeter: %s \n",fCalorimeter.Data()) ;
373 sprintf(onePar,"Number of modules: %d \n",fNModules) ;
376 TObjString *oString= new TObjString(parList) ;
377 outputContainer->Add(oString);
379 return outputContainer;
382 //_________________________________________________________________________________________________________________________________________________
383 void AliAnaPi0::Print(const Option_t * /*opt*/) const
385 //Print some relevant parameters set for the analysis
386 printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
387 AliAnaPartCorrBaseClass::Print(" ");
389 printf("Number of bins in Centrality: %d \n",fNCentrBin) ;
390 printf("Number of bins in Z vert. pos: %d \n",fNZvertBin) ;
391 printf("Number of bins in Reac. Plain: %d \n",fNrpBin) ;
392 printf("Depth of event buffer: %d \n",fNmaxMixEv) ;
393 printf("Number of different PID used: %d \n",fNPID) ;
395 printf("Z vertex position: -%2.3f < z < %2.3f \n",fZvtxCut,fZvtxCut) ;
396 printf("Number of modules: %d \n",fNModules) ;
397 printf("Select pairs with their angle: %d \n",fUseAngleCut) ;
398 printf("------------------------------------------------------\n") ;
401 //____________________________________________________________________________________________________________________________________________________
402 Int_t AliAnaPi0::GetModuleNumber(AliAODPWG4Particle * particle)
404 //Get the EMCAL/PHOS module number that corresponds to this particle
407 if(fCalorimeter=="EMCAL"){
408 GetReader()->GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(particle->Eta(),particle->Phi(), absId);
410 printf("AliAnaPi0::GetModuleNumber() - EMCAL: cluster eta %f, phi %f, absid %d, SuperModule %d\n",
411 particle->Eta(), particle->Phi()*TMath::RadToDeg(),absId, GetReader()->GetEMCALGeometry()->GetSuperModuleNumber(absId));
412 return GetReader()->GetEMCALGeometry()->GetSuperModuleNumber(absId) ;
416 if(!strcmp((GetReader()->GetInputEvent())->GetName(),"AliESDEvent")) {
417 AliESDCaloCluster *cluster = ((AliESDEvent*)GetReader()->GetInputEvent())->GetCaloCluster(particle->GetCaloLabel(0));
418 if ( cluster->GetNCells() > 0) {
419 absId = cluster->GetCellAbsId(0);
421 printf("AliAnaPi0::GetModuleNumber(ESD) - PHOS: cluster eta %f, phi %f, e %f, e cluster %f, absId %d\n",
422 particle->Eta(), particle->Phi()*TMath::RadToDeg(), particle->E(), cluster->E(), absId);
427 AliAODCaloCluster *cluster = ((AliAODEvent*)GetReader()->GetInputEvent())->GetCaloCluster(particle->GetCaloLabel(0));
428 if ( cluster->GetNCells() > 0) {
429 absId = cluster->GetCellAbsId(0);
431 printf("AliAnaPi0::GetModuleNumber(AOD) - PHOS: cluster eta %f, phi %f, e %f, e cluster %f, absId %d\n",
432 particle->Eta(), particle->Phi()*TMath::RadToDeg(), particle->E(), cluster->E(), absId);
438 GetReader()->GetPHOSGeometry()->AbsToRelNumbering(absId,relId);
440 printf("PHOS: Module %d\n",relId[0]-1);
449 //____________________________________________________________________________________________________________________________________________________
450 void AliAnaPi0::MakeAnalysisFillHistograms()
452 //Process one event and extract photons from AOD branch
453 // filled with AliAnaPhoton and fill histos with invariant mass
455 //Apply some cuts on event: vertex position and centrality range
456 Int_t iRun=(GetReader()->GetInputEvent())->GetRunNumber() ;
457 if(IsBadRun(iRun)) return ;
459 Double_t vert[]={0,0,0} ; //vertex ;
460 GetReader()->GetVertex(vert);
461 if(vert[2]<-fZvtxCut || vert[2]> fZvtxCut) return ; //Event can not be used (vertex, centrality,... cuts not fulfilled)
463 //Get Centrality and calculate centrality bin
464 //Does not exist in ESD yet???????
465 Int_t curCentrBin=0 ;
467 //Get Reaction Plain position and calculate RP bin
468 //does not exist in ESD yet????
471 Int_t curZvertBin=(Int_t)(0.5*fNZvertBin*(vert[2]+fZvtxCut)/fZvtxCut) ;
473 fhEvents->Fill(curCentrBin+0.5,curZvertBin+0.5,curRPBin+0.5) ;
475 Int_t nPhot = GetInputAODBranch()->GetEntriesFast() ;
476 if(GetDebug() > 1) printf("AliAnaPi0::MakeAnalysisFillHistograms() - Photon entries %d\n", nPhot);
479 for(Int_t i1=0; i1<nPhot-1; i1++){
480 AliAODPWG4Particle * p1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i1)) ;
481 TLorentzVector photon1(p1->Px(),p1->Py(),p1->Pz(),p1->E());
483 module1 = GetModuleNumber(p1);
484 for(Int_t i2=i1+1; i2<nPhot; i2++){
485 AliAODPWG4Particle * p2 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i2)) ;
486 TLorentzVector photon2(p2->Px(),p2->Py(),p2->Pz(),p2->E());
488 module2 = GetModuleNumber(p2);
489 Double_t m = (photon1 + photon2).M() ;
490 Double_t pt = (photon1 + photon2).Pt();
491 Double_t a = TMath::Abs(p1->E()-p2->E())/(p1->E()+p2->E()) ;
493 printf("AliAnaPi0::MakeAnalysisFillHistograms() - Current Event: pT: photon1 %2.2f, photon2 %2.2f; Pair: pT %2.2f, mass %2.3f, a %f2.3\n",
494 p1->Pt(), p2->Pt(), pt,m,a);
496 //Check if opening angle is too large or too small compared to what is expected
497 Double_t angle = photon1.Angle(photon2.Vect());
498 //if(fUseAngleCut && !GetNeutralMesonSelection()->IsAngleInWindow((photon1+photon2).E(),angle)) continue;
499 //printf("angle %f\n",angle);
500 if(fUseAngleCut && angle < 0.1) continue;
501 fhRealOpeningAngle ->Fill(pt,angle);
502 fhRealCosOpeningAngle->Fill(pt,TMath::Cos(angle));
504 //Fill module dependent histograms
505 //if(module1==module2) printf("mod1 %d\n",module1);
506 if(module1==module2 && module1 >=0 && module1<fNModules)
507 fhReMod[module1]->Fill(pt,a,m) ;
509 for(Int_t ipid=0; ipid<fNPID; ipid++)
511 if((p1->IsPIDOK(ipid,AliCaloPID::kPhoton)) && (p2->IsPIDOK(ipid,AliCaloPID::kPhoton))){
512 fhRe1[curCentrBin*fNPID+ipid]->Fill(pt,a,m) ;
513 if(p1->DistToBad()>0 && p2->DistToBad()>0){
514 fhRe2[curCentrBin*fNPID+ipid]->Fill(pt,a,m) ;
515 if(p1->DistToBad()>1 && p2->DistToBad()>1){
516 fhRe3[curCentrBin*fNPID+ipid]->Fill(pt,a,m) ;
526 TList * evMixList=fEventsList[curCentrBin*fNZvertBin*fNrpBin+curZvertBin*fNrpBin+curRPBin] ;
527 Int_t nMixed = evMixList->GetSize() ;
528 for(Int_t ii=0; ii<nMixed; ii++){
529 TClonesArray* ev2= (TClonesArray*) (evMixList->At(ii));
530 Int_t nPhot2=ev2->GetEntriesFast() ;
532 if(GetDebug() > 1) printf("AliAnaPi0::MakeAnalysisFillHistograms() - Mixed event %d photon entries %d\n", ii, nPhot);
534 for(Int_t i1=0; i1<nPhot; i1++){
535 AliAODPWG4Particle * p1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i1)) ;
536 TLorentzVector photon1(p1->Px(),p1->Py(),p1->Pz(),p1->E());
537 for(Int_t i2=0; i2<nPhot2; i2++){
538 AliAODPWG4Particle * p2 = (AliAODPWG4Particle*) (ev2->At(i2)) ;
540 TLorentzVector photon2(p2->Px(),p2->Py(),p2->Pz(),p2->E());
541 m = (photon1+photon2).M() ;
542 Double_t pt = (photon1 + photon2).Pt();
543 Double_t a = TMath::Abs(p1->E()-p2->E())/(p1->E()+p2->E()) ;
545 //Check if opening angle is too large or too small compared to what is expected
546 Double_t angle = photon1.Angle(photon2.Vect());
547 //if(fUseAngleCut && !GetNeutralMesonSelection()->IsAngleInWindow((photon1+photon2).E(),angle)) continue;
548 if(fUseAngleCut && angle < 0.1) continue;
551 printf("AliAnaPi0::MakeAnalysisFillHistograms() - Mixed Event: pT: photon1 %2.2f, photon2 %2.2f; Pair: pT %2.2f, mass %2.3f, a %f2.3\n",
552 p1->Pt(), p2->Pt(), pt,m,a);
553 for(Int_t ipid=0; ipid<fNPID; ipid++){
554 if((p1->IsPIDOK(ipid,AliCaloPID::kPhoton)) && (p2->IsPIDOK(ipid,AliCaloPID::kPhoton))){
555 fhMi1[curCentrBin*fNPID+ipid]->Fill(pt,a,m) ;
556 if(p1->DistToBad()>0 && p2->DistToBad()>0){
557 fhMi2[curCentrBin*fNPID+ipid]->Fill(pt,a,m) ;
558 if(p1->DistToBad()>1 && p2->DistToBad()>1){
559 fhMi3[curCentrBin*fNPID+ipid]->Fill(pt,a,m) ;
569 TClonesArray *currentEvent = new TClonesArray(*GetInputAODBranch());
570 //Add current event to buffer and Remove redandant events
571 if(currentEvent->GetEntriesFast()>0){
572 evMixList->AddFirst(currentEvent) ;
573 currentEvent=0 ; //Now list of particles belongs to buffer and it will be deleted with buffer
574 if(evMixList->GetSize()>=fNmaxMixEv)
576 TClonesArray * tmp = (TClonesArray*) (evMixList->Last()) ;
577 evMixList->RemoveLast() ;
582 delete currentEvent ;
587 if(IsDataMC() && GetReader()->ReadStack()){
588 AliStack * stack = GetMCStack();
589 if(stack && (IsDataMC() || (GetReader()->GetDataType() == AliCaloTrackReader::kMC)) ){
590 for(Int_t i=0 ; i<stack->GetNprimary(); i++){
591 TParticle * prim = stack->Particle(i) ;
592 if(prim->GetPdgCode() == 111){
593 Double_t pi0Pt = prim->Pt() ;
594 //printf("pi0, pt %2.2f\n",pi0Pt);
595 if(prim->Energy() == TMath::Abs(prim->Pz())) continue ; //Protection against floating point exception
596 Double_t pi0Y = 0.5*TMath::Log((prim->Energy()-prim->Pz())/(prim->Energy()+prim->Pz())) ;
597 Double_t phi = TMath::RadToDeg()*prim->Phi() ;
598 if(TMath::Abs(pi0Y) < 0.5){
599 fhPrimPt->Fill(pi0Pt) ;
601 fhPrimY ->Fill(pi0Y) ;
602 fhPrimPhi->Fill(phi) ;
604 //Check if both photons hit Calorimeter
605 Int_t iphot1=prim->GetFirstDaughter() ;
606 Int_t iphot2=prim->GetLastDaughter() ;
607 if(iphot1>-1 && iphot1<stack->GetNtrack() && iphot2>-1 && iphot2<stack->GetNtrack()){
608 TParticle * phot1 = stack->Particle(iphot1) ;
609 TParticle * phot2 = stack->Particle(iphot2) ;
610 if(phot1 && phot2 && phot1->GetPdgCode()==22 && phot2->GetPdgCode()==22){
611 //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",
612 // phot1->Pt(), phot1->Phi()*180./3.1415, phot1->Eta(), phot2->Pt(), phot2->Phi()*180./3.1415, phot2->Eta());
614 TLorentzVector lv1, lv2;
615 phot1->Momentum(lv1);
616 phot2->Momentum(lv2);
618 Bool_t inacceptance = kFALSE;
619 if(fCalorimeter == "PHOS"){
620 if(GetReader()->GetPHOSGeometry() && gGeoManager){
623 if(GetReader()->GetPHOSGeometry()->ImpactOnEmc(phot1,mod,z,x) && GetReader()->GetPHOSGeometry()->ImpactOnEmc(phot2,mod,z,x))
624 inacceptance = kTRUE;
625 if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
629 if(GetFiducialCut()->IsInFiducialCut(lv1,fCalorimeter) && GetFiducialCut()->IsInFiducialCut(lv2,fCalorimeter))
630 inacceptance = kTRUE ;
631 if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
635 else if(fCalorimeter == "EMCAL" && gGeoManager){
636 if(GetReader()->GetEMCALGeometry()){
637 if(GetReader()->GetEMCALGeometry()->Impact(phot1) && GetReader()->GetEMCALGeometry()->Impact(phot2))
638 inacceptance = kTRUE;
639 if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
642 if(GetFiducialCut()->IsInFiducialCut(lv1,fCalorimeter) && GetFiducialCut()->IsInFiducialCut(lv2,fCalorimeter))
643 inacceptance = kTRUE ;
644 if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
650 fhPrimAccPt->Fill(pi0Pt) ;
651 fhPrimAccPhi->Fill(phi) ;
652 fhPrimAccY->Fill(pi0Y) ;
653 Double_t angle = lv1.Angle(lv2.Vect());
654 fhPrimOpeningAngle ->Fill(pi0Pt,angle);
655 fhPrimCosOpeningAngle->Fill(pi0Pt,TMath::Cos(angle));
659 }//Check daughters exist
662 }//stack exists and data is MC
664 else if(GetReader()->ReadAODMCParticles()){
665 if(GetDebug() >= 0) printf("AliAnaPi0::MakeAnalysisFillHistograms() - Acceptance calculation with MCParticles not implemented yet\n");
670 //________________________________________________________________________
671 void AliAnaPi0::ReadHistograms(TList* outputList)
673 // Needed when Terminate is executed in distributed environment
674 // Refill analysis histograms of this class with corresponding histograms in output list.
676 // Histograms of this analsys are kept in the same list as other analysis, recover the position of
677 // the first one and then add the next.
678 Int_t index = outputList->IndexOf(outputList->FindObject(GetAddedHistogramsStringToName()+"hRe_cen0_pid0_dist1"));
680 if(!fhRe1) fhRe1 = new TH3D*[fNCentrBin*fNPID] ;
681 if(!fhRe2) fhRe2 = new TH3D*[fNCentrBin*fNPID] ;
682 if(!fhRe3) fhRe3 = new TH3D*[fNCentrBin*fNPID] ;
683 if(!fhMi1) fhMi1 = new TH3D*[fNCentrBin*fNPID] ;
684 if(!fhMi2) fhMi2 = new TH3D*[fNCentrBin*fNPID] ;
685 if(!fhMi3) fhMi3 = new TH3D*[fNCentrBin*fNPID] ;
686 if(!fhReMod) fhReMod = new TH3D*[fNModules] ;
688 for(Int_t ic=0; ic<fNCentrBin; ic++){
689 for(Int_t ipid=0; ipid<fNPID; ipid++){
690 fhRe1[ic*fNPID+ipid] = (TH3D*) outputList->At(index++);
691 fhMi1[ic*fNPID+ipid] = (TH3D*) outputList->At(index++);
692 fhRe2[ic*fNPID+ipid] = (TH3D*) outputList->At(index++);
693 fhMi2[ic*fNPID+ipid] = (TH3D*) outputList->At(index++);
694 fhRe3[ic*fNPID+ipid] = (TH3D*) outputList->At(index++);
695 fhMi3[ic*fNPID+ipid] = (TH3D*) outputList->At(index++);
699 fhEvents = (TH3D *) outputList->At(index++);
701 //Histograms filled only if MC data is requested
702 if(IsDataMC() || (GetReader()->GetDataType() == AliCaloTrackReader::kMC) ){
703 fhPrimPt = (TH1D*) outputList->At(index++);
704 fhPrimAccPt = (TH1D*) outputList->At(index++);
705 fhPrimY = (TH1D*) outputList->At(index++);
706 fhPrimAccY = (TH1D*) outputList->At(index++);
707 fhPrimPhi = (TH1D*) outputList->At(index++);
708 fhPrimAccPhi = (TH1D*) outputList->At(index++);
711 for(Int_t imod=0; imod < fNModules; imod++)
712 fhReMod[imod] = (TH3D*) outputList->At(index++);
717 //____________________________________________________________________________________________________________________________________________________
718 void AliAnaPi0::Terminate(TList* outputList)
720 //Do some calculations and plots from the final histograms.
722 printf(" *** %s Terminate:\n", GetName()) ;
724 //Recover histograms from output histograms list, needed for distributed analysis.
725 ReadHistograms(outputList);
728 printf("AliAnaPi0::Terminate() - Error: Remote output histograms not imported in AliAnaPi0 object");
732 printf("AliAnaPi0::Terminate() Mgg Real : %5.3f , RMS : %5.3f \n", fhRe1[0]->GetMean(), fhRe1[0]->GetRMS() ) ;
735 sprintf(nameIM,"AliAnaPi0_%s_cPt",fCalorimeter.Data());
736 TCanvas * cIM = new TCanvas(nameIM, "", 400, 10, 600, 700) ;
741 TH1D * hIMAllPt = (TH1D*) fhRe1[0]->ProjectionZ(Form("IMPtAll_%s",fCalorimeter.Data()));
742 hIMAllPt->SetLineColor(2);
743 hIMAllPt->SetTitle("No cut on p_{T, #gamma#gamma} ");
747 TH3F * hRe1Pt5 = (TH3F*)fhRe1[0]->Clone(Form("IMPt5_%s",fCalorimeter.Data()));
748 hRe1Pt5->GetXaxis()->SetRangeUser(0,5);
749 TH1D * hIMPt5 = (TH1D*) hRe1Pt5->Project3D(Form("IMPt5_%s_pz",fCalorimeter.Data()));
750 hIMPt5->SetLineColor(2);
751 hIMPt5->SetTitle("0 < p_{T, #gamma#gamma} < 5 GeV/c");
755 TH3F * hRe1Pt10 = (TH3F*)fhRe1[0]->Clone(Form("IMPt10_%s",fCalorimeter.Data()));
756 hRe1Pt10->GetXaxis()->SetRangeUser(5,10);
757 TH1D * hIMPt10 = (TH1D*) hRe1Pt10->Project3D(Form("IMPt10_%s_pz",fCalorimeter.Data()));
758 hIMPt10->SetLineColor(2);
759 hIMPt10->SetTitle("5 < p_{T, #gamma#gamma} < 10 GeV/c");
763 TH3F * hRe1Pt20 = (TH3F*)fhRe1[0]->Clone(Form("IMPt20_%s",fCalorimeter.Data()));
764 hRe1Pt20->GetXaxis()->SetRangeUser(10,20);
765 TH1D * hIMPt20 = (TH1D*) hRe1Pt20->Project3D(Form("IMPt20_%s_pz",fCalorimeter.Data()));
766 hIMPt20->SetLineColor(2);
767 hIMPt20->SetTitle("10 < p_{T, #gamma#gamma} < 20 GeV/c");
771 sprintf(nameIMF,"AliAnaPi0_%s_Mgg.eps",fCalorimeter.Data());
775 sprintf(namePt,"AliAnaPi0_%s_cPt",fCalorimeter.Data());
776 TCanvas * cPt = new TCanvas(namePt, "", 400, 10, 600, 700) ;
781 TH1D * hPt = (TH1D*) fhRe1[0]->Project3D("x");
782 hPt->SetLineColor(2);
783 hPt->SetTitle("No cut on M_{#gamma#gamma} ");
787 TH3F * hRe1IM1 = (TH3F*)fhRe1[0]->Clone(Form("Pt1_%s",fCalorimeter.Data()));
788 hRe1IM1->GetZaxis()->SetRangeUser(0.05,0.21);
789 TH1D * hPtIM1 = (TH1D*) hRe1IM1->Project3D("x");
790 hPtIM1->SetLineColor(2);
791 hPtIM1->SetTitle("0.05 < M_{#gamma#gamma} < 0.21 GeV/c^{2}");
795 TH3F * hRe1IM2 = (TH3F*)fhRe1[0]->Clone(Form("Pt2_%s",fCalorimeter.Data()));
796 hRe1IM2->GetZaxis()->SetRangeUser(0.09,0.17);
797 TH1D * hPtIM2 = (TH1D*) hRe1IM2->Project3D("x");
798 hPtIM2->SetLineColor(2);
799 hPtIM2->SetTitle("0.09 < M_{#gamma#gamma} < 0.17 GeV/c^{2}");
803 TH3F * hRe1IM3 = (TH3F*)fhRe1[0]->Clone(Form("Pt3_%s",fCalorimeter.Data()));
804 hRe1IM3->GetZaxis()->SetRangeUser(0.11,0.15);
805 TH1D * hPtIM3 = (TH1D*) hRe1IM1->Project3D("x");
806 hPtIM3->SetLineColor(2);
807 hPtIM3->SetTitle("0.11 < M_{#gamma#gamma} < 0.15 GeV/c^{2}");
811 sprintf(namePtF,"AliAnaPi0_%s_Pt.eps",fCalorimeter.Data());
816 sprintf(line, ".!tar -zcf %s_%s.tar.gz *.eps", GetName(),fCalorimeter.Data()) ;
817 gROOT->ProcessLine(line);
818 sprintf(line, ".!rm -fR AliAnaPi0_%s*.eps",fCalorimeter.Data());
819 gROOT->ProcessLine(line);
821 printf(" AliAnaPi0::Terminate() - !! All the eps files are in %s_%s.tar.gz !!!\n", GetName(), fCalorimeter.Data());