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