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