]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG4/PartCorrDep/AliAnaPi0EbE.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Several changes:
[u/mrichter/AliRoot.git] / PWG4 / PartCorrDep / AliAnaPi0EbE.cxx
CommitLineData
477d6cee 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 iGetEntriesFast(s 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: AliAnaPi0EbE.cxx 28688 2008-09-11 15:04:07Z gconesab $ */
16
17//_________________________________________________________________________
18// Class for the analysis of high pT pi0 event by event
19// Pi0 identified by one of the following:
20// -Invariant mass of 2 cluster in calorimeter
21// -Shower shape analysis in calorimeter
22// -Invariant mass of one cluster in calorimeter and one photon reconstructed in TPC (in near future)
23//
24// -- Author: Gustavo Conesa (LNF-INFN) & Raphaelle Ichou (SUBATECH)
25//////////////////////////////////////////////////////////////////////////////
26
27
28// --- ROOT system ---
29#include <TList.h>
30#include <TClonesArray.h>
31#include <TObjString.h>
32#include <TH2F.h>
33//#include "Riostream.h"
34
35// --- Analysis system ---
36#include "AliAnaPi0EbE.h"
37#include "AliCaloTrackReader.h"
38#include "AliIsolationCut.h"
39#include "AliNeutralMesonSelection.h"
40#include "AliCaloPID.h"
41#include "AliMCAnalysisUtils.h"
42#include "AliAODPWG4ParticleCorrelation.h"
43#include "AliStack.h"
44#include "AliFidutialCut.h"
45#include "TParticle.h"
46#include "AliAODCaloCluster.h"
47#include "AliAODEvent.h"
48
49ClassImp(AliAnaPi0EbE)
50
691bdd02 51//____________________________________________________________________________
477d6cee 52AliAnaPi0EbE::AliAnaPi0EbE() :
53AliAnaPartCorrBaseClass(), fAnaType(kIMCalo),fCalorimeter(""),
54fMinDist(0.),fMinDist2(0.),fMinDist3(0.),
55fInputAODGammaConv(0x0),fInputAODGammaConvName(""),
56fhPtPi0(0),fhPhiPi0(0),fhEtaPi0(0),
57fhPtMCNoPi0(0),fhPhiMCNoPi0(0),fhEtaMCNoPi0(0),
58fhPtMCPi0(0),fhPhiMCPi0(0),fhEtaMCPi0(0)
59{
60 //default ctor
61
62 //Initialize parameters
63 InitParameters();
64
65}
66
67//____________________________________________________________________________
68AliAnaPi0EbE::AliAnaPi0EbE(const AliAnaPi0EbE & p) :
69AliAnaPartCorrBaseClass(p), fAnaType(p.fAnaType), fCalorimeter(p.fCalorimeter),
70fMinDist(p.fMinDist),fMinDist2(p.fMinDist2), fMinDist3(p.fMinDist3),
71fInputAODGammaConv(new TClonesArray(*p.fInputAODGammaConv)),
72fInputAODGammaConvName(p.fInputAODGammaConvName),
73fhPtPi0(p.fhPtPi0),fhPhiPi0(p.fhPhiPi0),fhEtaPi0(p.fhEtaPi0),
74fhPtMCNoPi0(p.fhPtMCNoPi0),fhPhiMCNoPi0(p.fhPhiMCNoPi0),fhEtaMCNoPi0(p.fhEtaMCNoPi0),
75fhPtMCPi0(p.fhPtMCPi0),fhPhiMCPi0(p.fhPhiMCPi0),fhEtaMCPi0(p.fhEtaMCPi0)
76{
77 // cpy ctor
78}
79
80//_________________________________________________________________________
81AliAnaPi0EbE & AliAnaPi0EbE::operator = (const AliAnaPi0EbE & p)
82{
83 // assignment operator
84
85 if(&p == this) return *this;
86
87 fAnaType = p.fAnaType ;
88 fCalorimeter = p.fCalorimeter ;
89 fMinDist = p.fMinDist;
90 fMinDist2 = p.fMinDist2;
91 fMinDist3 = p.fMinDist3;
92
93 fInputAODGammaConv = new TClonesArray(*p.fInputAODGammaConv) ;
94 fInputAODGammaConvName = p.fInputAODGammaConvName ;
95
96 fhPtPi0 = p.fhPtPi0; fhPhiPi0 = p.fhPhiPi0; fhEtaPi0 = p.fhEtaPi0;
97 fhPtMCNoPi0 = p.fhPtMCNoPi0; fhPhiMCNoPi0 = p.fhPhiMCNoPi0; fhEtaMCNoPi0 = p.fhEtaMCPi0;
98 fhPtMCPi0 = p.fhPtMCPi0; fhPhiMCPi0 = p.fhPhiMCPi0; fhEtaMCPi0 = p.fhEtaMCPi0;
99
100 return *this;
101
102}
103
104//____________________________________________________________________________
105AliAnaPi0EbE::~AliAnaPi0EbE()
106{
107 //dtor
108 if(fInputAODGammaConv){
109 fInputAODGammaConv->Clear() ;
110 delete fInputAODGammaConv ;
111 }
112}
113
114//________________________________________________________________________
115TList * AliAnaPi0EbE::GetCreateOutputObjects()
116{
117 // Create histograms to be saved in output file and
118 // store them in outputContainer
119 TList * outputContainer = new TList() ;
120 outputContainer->SetName("Pi0EbEHistos") ;
121
122 Int_t nptbins = GetHistoNPtBins();
123 Int_t nphibins = GetHistoNPhiBins();
124 Int_t netabins = GetHistoNEtaBins();
125 Float_t ptmax = GetHistoPtMax();
126 Float_t phimax = GetHistoPhiMax();
127 Float_t etamax = GetHistoEtaMax();
128 Float_t ptmin = GetHistoPtMin();
129 Float_t phimin = GetHistoPhiMin();
130 Float_t etamin = GetHistoEtaMin();
131
132 fhPtPi0 = new TH1F("hPtPi0","Number of identified #pi^{0} decay",nptbins,ptmin,ptmax);
133 fhPtPi0->SetYTitle("N");
134 fhPtPi0->SetXTitle("p_{T #pi^{0}}(GeV/c)");
135 outputContainer->Add(fhPtPi0) ;
136
137 fhPhiPi0 = new TH2F
138 ("hPhiPi0","#phi_{#pi^{0}}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
139 fhPhiPi0->SetYTitle("#phi");
140 fhPhiPi0->SetXTitle("p_{T #pi^{0}} (GeV/c)");
141 outputContainer->Add(fhPhiPi0) ;
142
143 fhEtaPi0 = new TH2F
144 ("hEtaPi0","#phi_{#pi^{0}}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
145 fhEtaPi0->SetYTitle("#eta");
146 fhEtaPi0->SetXTitle("p_{T #pi^{0}} (GeV/c)");
147 outputContainer->Add(fhEtaPi0) ;
148
149 if(IsDataMC()) {
150 if((GetReader()->GetDataType() == AliCaloTrackReader::kMC && fAnaType!=kSSCalo) ||
151 GetReader()->GetDataType() != AliCaloTrackReader::kMC){
152
153 fhPtMCPi0 = new TH1F("hPtMCPi0","Identified pi0 from pi0",nptbins,ptmin,ptmax);
154 fhPtMCPi0->SetYTitle("N");
155 fhPtMCPi0->SetXTitle("p_{T #pi^{0}}(GeV/c)");
156 outputContainer->Add(fhPtMCPi0) ;
157
158 fhPhiMCPi0 = new TH2F
159 ("hPhiMCPi0","Identified pi0 from pi0",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
160 fhPhiMCPi0->SetYTitle("#phi");
161 fhPhiMCPi0->SetXTitle("p_{T #pi^{0}} (GeV/c)");
162 outputContainer->Add(fhPhiMCPi0) ;
163
164 fhEtaMCPi0 = new TH2F
165 ("hEtaMCPi0","Identified pi0 from pi0",nptbins,ptmin,ptmax,netabins,etamin,etamax);
166 fhEtaMCPi0->SetYTitle("#eta");
167 fhEtaMCPi0->SetXTitle("p_{T #pi^{0}} (GeV/c)");
168 outputContainer->Add(fhEtaMCPi0) ;
169
170 fhPtMCNoPi0 = new TH1F("hPtMCNoPi0","Identified pi0 not from pi0",nptbins,ptmin,ptmax);
171 fhPtMCNoPi0->SetYTitle("N");
172 fhPtMCNoPi0->SetXTitle("p_{T #pi^{0}}(GeV/c)");
173 outputContainer->Add(fhPtMCNoPi0) ;
174
175 fhPhiMCNoPi0 = new TH2F
176 ("hPhiMCNoPi0","Identified pi0 not from pi0",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
177 fhPhiMCNoPi0->SetYTitle("#phi");
178 fhPhiMCNoPi0->SetXTitle("p_{T #pi^{0}} (GeV/c)");
179 outputContainer->Add(fhPhiMCNoPi0) ;
180
181 fhEtaMCNoPi0 = new TH2F
182 ("hEtaMCNoPi0","Identified pi0 not from pi0",nptbins,ptmin,ptmax,netabins,etamin,etamax);
183 fhEtaMCNoPi0->SetYTitle("#eta");
184 fhEtaMCNoPi0->SetXTitle("p_{T #pi^{0}} (GeV/c)");
185 outputContainer->Add(fhEtaMCNoPi0) ;
186
187 }
188 }//Histos with MC
189
190
191 //Keep neutral meson selection histograms if requiered
192 //Setting done in AliNeutralMesonSelection
193
194 if(fAnaType!=kSSCalo && GetNeutralMesonSelection()){
195
196 TList * nmsHistos = GetNeutralMesonSelection()->GetCreateOutputObjects() ;
197 if(GetNeutralMesonSelection()->AreNeutralMesonSelectionHistosKept())
198 for(Int_t i = 0; i < nmsHistos->GetEntries(); i++) outputContainer->Add(nmsHistos->At(i)) ;
199 }
200
201 //Save parameters used for analysis
202 TString parList ; //this will be list of parameters used for this analysis.
203 char onePar[255] ;
204
205 sprintf(onePar,"--- AliAnaPi0EbE ---\n") ;
206 parList+=onePar ;
207 sprintf(onePar,"fAnaType=%d (Pi0 selection type) \n",fAnaType) ;
208 parList+=onePar ;
209
210 if(fAnaType == kSSCalo){
211 sprintf(onePar,"Calorimeter: %s\n",fCalorimeter.Data()) ;
212 parList+=onePar ;
213 sprintf(onePar,"fMinDist =%2.2f (Minimal distance to bad channel to accept cluster) \n",fMinDist) ;
214 parList+=onePar ;
215 sprintf(onePar,"fMinDist2=%2.2f (Cuts on Minimal distance to study acceptance evaluation) \n",fMinDist2) ;
216 parList+=onePar ;
217 sprintf(onePar,"fMinDist3=%2.2f (One more cut on distance used for acceptance-efficiency study) \n",fMinDist3) ;
218 parList+=onePar ;
219 }
220
221 //Get parameters set in base class.
222 parList += GetBaseParametersList() ;
223
224 //Get parameters set in PID class.
225 if(fAnaType == kSSCalo) parList += GetCaloPID()->GetPIDParametersList() ;
226
227 TObjString *oString= new TObjString(parList) ;
228 outputContainer->Add(oString);
229
230 return outputContainer ;
231
232}
233
234//__________________________________________________________________
235void AliAnaPi0EbE::MakeAnalysisFillAOD()
236{
237 //Do analysis and fill aods
238
239 switch(fAnaType)
240 {
241 case kIMCalo:
242 MakeInvMassInCalorimeter();
243 break;
244
245 case kSSCalo:
246 MakeShowerShapeIdentification();
247 break;
248
249 case kIMCaloTracks:
250 MakeInvMassInCalorimeterAndCTS();
251 break;
252
253 }
254}
255
256//__________________________________________________________________
257void AliAnaPi0EbE::MakeInvMassInCalorimeter()
258{
259 //Do analysis and fill aods
260 //Search for the photon decay in calorimeters
261 //Read photon list from AOD, produced in class AliAnaPhoton
262 //Check if 2 photons have the mass of the pi0.
263
264 TLorentzVector mom1;
265 TLorentzVector mom2;
266 TLorentzVector mom ;
267 Int_t tag1 =-1;
268 Int_t tag2 =-1;
269 Int_t tag = -1;
270
271 if(!GetInputAODBranch()){
272 printf("AliAnaPi0EbE::kIMCalo:: No input calo photons in AOD with name branch < %s >, STOP \n",GetInputAODName().Data());
273 abort();
274 }
275 for(Int_t iphoton = 0; iphoton < GetInputAODBranch()->GetEntriesFast(); iphoton++){
276 AliAODPWG4Particle * photon1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(iphoton));
277 mom1 = *(photon1->Momentum());
278
279
280 for(Int_t jphoton = iphoton+1; jphoton < GetInputAODBranch()->GetEntriesFast()-1; jphoton++){
281
282 AliAODPWG4ParticleCorrelation * photon2 = (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(jphoton));
283 mom2 = *(photon2->Momentum());
284
285 //Select good pair (good phi, pt cuts, aperture and invariant mass)
286 if(GetNeutralMesonSelection()->SelectPair(mom1, mom2))
287 {
288 if(GetDebug()>1)
289 printf("AliAnaPi0EbE::kIMCalo::Selected gamma pair: pt %f, phi %f, eta%f \n",(mom1+mom2).Pt(), (mom1+mom2).Phi()*180./3.1416, (mom1+mom2).Eta());
290
291 //Play with the MC stack if available
292 if(IsDataMC()){
293 //Check origin of the candidates
294 tag1 = GetMCAnalysisUtils()->CheckOrigin(photon1->GetLabel(), GetMCStack());
295 tag2 = GetMCAnalysisUtils()->CheckOrigin(photon2->GetLabel(), GetMCStack());
296
297 if(GetDebug() > 0) printf("AliAnaPi0EbE::kIMCalo::Origin of: photon1 %d; photon2 %d \n",tag1, tag2);
298 if(tag1 == AliMCAnalysisUtils::kMCPi0Decay && tag2 == AliMCAnalysisUtils::kMCPi0Decay){
299
300 //Check if pi0 mother is the same
301 Int_t label1 = photon1->GetLabel();
302 TParticle * mother1 = GetMCStack()->Particle(label1);//photon in kine tree
303 label1 = mother1->GetFirstMother();
304 //mother1 = GetMCStack()->Particle(label1);//pi0
305
306 Int_t label2 = photon2->GetLabel();
307 TParticle * mother2 = GetMCStack()->Particle(label2);//photon in kine tree
308 label2 = mother2->GetFirstMother();
309 //mother2 = GetMCStack()->Particle(label2);//pi0
310
311 //printf("mother1 %d, mother2 %d\n",label1,label2);
312 if(label1 == label2)
313 tag = AliMCAnalysisUtils::kMCPi0;
314 }
315 }//Work with stack also
316
317 //Create AOD for analysis
318 mom = mom1+mom2;
319 AliAODPWG4Particle pi0 = AliAODPWG4Particle(mom);
320 //pi0.SetLabel(calo->GetLabel(0));
321 pi0.SetPdg(AliCaloPID::kPi0);
322 pi0.SetDetector(photon1->GetDetector());
323 pi0.SetTag(tag);
324 //Set the indeces of the original caloclusters
325 pi0.SetCaloLabel(photon1->GetCaloLabel(0), photon2->GetCaloLabel(0));
326 AddAODParticle(pi0);
327 }//pi0
328 }//2n photon loop
329
330 }//1st photon loop
331
332 if(GetDebug() > 1) printf("AliAnaPi0EbE::kIMCalo::End fill AODs \n");
333
334}
335
336//__________________________________________________________________
337void AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS()
338{
339 //Do analysis and fill aods
340 //Search for the photon decay in calorimeters
341 //Read photon list from AOD, produced in class AliAnaPhoton and AliGammaConversion
342 //Check if 2 photons have the mass of the pi0.
343
344 TLorentzVector mom1;
345 TLorentzVector mom2;
346 TLorentzVector mom ;
347 Int_t tag1 =-1;
348 Int_t tag2 =-1;
349 Int_t tag = -1;
350
351 if(!GetInputAODBranch()){
352 printf("AliAnaPi0EbE::kIMCaloCTS: No input calo photons in AOD branch with name < %s > , STOP\n",GetInputAODName().Data());
353 abort();
354 }
355 for(Int_t iphoton = 0; iphoton < GetInputAODBranch()->GetEntriesFast(); iphoton++){
356 AliAODPWG4Particle * photon1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(iphoton));
357 mom1 = *(photon1->Momentum());
358
359 //Play with the MC stack if available
360 fInputAODGammaConv = (TClonesArray *) GetReader()->GetOutputEvent()->FindListObject(fInputAODGammaConvName);
361 if(!fInputAODGammaConv) {
362 printf("AliAnaPi0EbE::kIMCaloCTS: No input gamma conversions in AOD branch with name < %s >, STOP",fInputAODGammaConvName.Data());
363 abort();
364 }
365 for(Int_t jphoton = iphoton+1; jphoton < fInputAODGammaConv->GetEntriesFast()-1; jphoton++){
366 AliAODPWG4ParticleCorrelation * photon2 = (AliAODPWG4ParticleCorrelation*) (fInputAODGammaConv->At(jphoton));
367 mom2 = *(photon2->Momentum());
368 //Select good pair (good phi, pt cuts, aperture and invariant mass)
369 if(GetNeutralMesonSelection()->SelectPair(mom1, mom2)){
370 if(GetDebug() > 1) printf("AliAnaPi0EbE::kIMCaloCTS::Selected gamma pair: pt %f, phi %f, eta%f\n",(mom1+mom2).Pt(), (mom1+mom2).Phi()*180./3.1416, (mom1+mom2).Eta());
371
372 if(IsDataMC()){
373 //Check origin of the candidates
374 tag1 = GetMCAnalysisUtils()->CheckOrigin(photon1->GetLabel(), GetMCStack());
375 tag2 = GetMCAnalysisUtils()->CheckOrigin(photon2->GetLabel(), GetMCStack());
376 if(GetDebug() > 0) printf("AliAnaPi0EbE::kIMCaloCTS::Origin of: photon1 %d; photon2 %d \n",tag1, tag2);
377 if(tag1 == AliMCAnalysisUtils::kMCPi0Decay && tag2 == AliMCAnalysisUtils::kMCPi0Decay){
378 //Check if pi0 mother is the same
379 Int_t label1 = photon1->GetLabel();
380 TParticle * mother1 = GetMCStack()->Particle(label1);//photon in kine tree
381 label1 = mother1->GetFirstMother();
382 //mother1 = GetMCStack()->Particle(label1);//pi0
383
384 Int_t label2 = photon2->GetLabel();
385 TParticle * mother2 = GetMCStack()->Particle(label2);//photon in kine tree
386 label2 = mother2->GetFirstMother();
387 //mother2 = GetMCStack()->Particle(label2);//pi0
388
389 //printf("mother1 %d, mother2 %d\n",label1,label2);
390 if(label1 == label2)
391 tag = AliMCAnalysisUtils::kMCPi0;
392 }
393 }//Work with stack also
394
395 //Create AOD for analysis
396 mom = mom1+mom2;
397 AliAODPWG4ParticleCorrelation pi0 = AliAODPWG4ParticleCorrelation(mom);
398 //pi0.SetLabel(calo->GetLabel(0));
399 pi0.SetPdg(AliCaloPID::kPi0);
400 pi0.SetDetector(photon1->GetDetector());
401 pi0.SetTag(tag);
402 //Set the indeces of the original tracks or caloclusters
403 pi0.SetCaloLabel(photon1->GetCaloLabel(0), -1);
404 pi0.SetTrackLabel(photon2->GetTrackLabel(0), photon2->GetTrackLabel(1));
405 AddAODParticle(pi0);
406 }//pi0
407 }//2n photon loop
408
409 }//1st photon loop
410
411 if(GetDebug() > 1) printf("AliAnaPi0EbE::kIMCaloCTS::End fill AODs \n");
412
413}
414
415
416//__________________________________________________________________
417void AliAnaPi0EbE::MakeShowerShapeIdentification()
418{
419 //Search for pi0 in fCalorimeter with shower shape analysis
420
421 TRefArray * pl = new TRefArray;
422
423 //Get vertex for photon momentum calculation
424 Double_t vertex[]={0,0,0} ; //vertex ;
425 if(!GetReader()->GetDataType()== AliCaloTrackReader::kMC) GetReader()->GetVertex(vertex);
426
427 //Select the Calorimeter of the photon
428 if(fCalorimeter == "PHOS")
429 pl = GetAODPHOS();
430 else if (fCalorimeter == "EMCAL")
431 pl = GetAODEMCAL();
432 //Fill AODCaloClusters and AODParticle with PHOS aods
433 TLorentzVector mom ;
434 for(Int_t icalo = 0; icalo < pl->GetEntriesFast(); icalo++){
435 AliAODCaloCluster * calo = (AliAODCaloCluster*) (pl->At(icalo));
436
437 //Cluster selection, not charged, with pi0 id and in fidutial cut
438 //Get Momentum vector,
439 calo->GetMomentum(mom,vertex);//Assume that come from vertex in straight line
440 //If too small or big pt, skip it
441 if(mom.Pt() < GetMinPt() || mom.Pt() > GetMaxPt() ) continue ;
442 //Check acceptance selection
443 if(IsFidutialCutOn()){
444 Bool_t in = GetFidutialCut()->IsInFidutialCut(mom,fCalorimeter) ;
445 if(! in ) continue ;
446 }
447
448 //Create AOD for analysis
449 AliAODPWG4Particle aodpi0 = AliAODPWG4Particle(mom);
450 aodpi0.SetLabel(calo->GetLabel(0));
451 //Set the indeces of the original caloclusters
452 aodpi0.SetCaloLabel(calo->GetID(),-1);
453 aodpi0.SetDetector(fCalorimeter);
454 if(GetDebug() > 1)
455 printf("AliAnaPi0EbE::kSSCalo::::FillAOD: Min pt cut and fidutial cut passed: pt %3.2f, phi %2.2f, eta %1.2f\n",aodpi0.Pt(),aodpi0.Phi(),aodpi0.Eta());
456
457 //Check Distance to Bad channel, set bit.
458 Double_t distBad=calo->GetDistToBadChannel() ; //Distance to bad channel
459 if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;
460 if(distBad < fMinDist) //In bad channel (PHOS cristal size 2.2x2.2 cm)
461 continue ;
462
463 if(GetDebug() > 1) printf("AliAnaPi0EbE::kSSCalo::::FillAOD: Bad channel cut passed %4.2f\n",distBad);
464
465 if(distBad > fMinDist3) aodpi0.SetDistToBad(2) ;
466 else if(distBad > fMinDist2) aodpi0.SetDistToBad(1) ;
467 else aodpi0.SetDistToBad(0) ;
468
469 //Check PID
470 //PID selection or bit setting
471 if(GetReader()->GetDataType() == AliCaloTrackReader::kMC){
472 //Get most probable PID, check PID weights (in MC this option is mandatory)
473 aodpi0.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,calo->PID(),mom.E()));//PID with weights
474 if(GetDebug() > 1)
475 printf("AliAnaPi0EbE::kSSCalo::::FillAOD: PDG of identified particle %d\n",aodpi0.GetPdg());
476 //If primary is not pi0, skip it.
477 if(aodpi0.GetPdg() != AliCaloPID::kPi0) continue ;
478 }
479 else if(IsCaloPIDOn()){
480 //Skip matched clusters with tracks
481 if(calo->GetNTracksMatched() > 0) continue ;
482
483 //Get most probable PID, 2 options check PID weights
484 //or redo PID, recommended option for EMCal.
485 if(!IsCaloPIDRecalculationOn())
486 aodpi0.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,calo->PID(),mom.E()));//PID with weights
487 else
488 aodpi0.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,mom,calo));//PID recalculated
489
490 if(GetDebug() > 1) printf("AliAnaPi0EbE::kSSCalo::::FillAOD: PDG of identified particle %d\n",aodpi0.GetPdg());
491
492 //If cluster does not pass pid, not pi0, skip it.
493 if(aodpi0.GetPdg() != AliCaloPID::kPi0) continue ;
494
495 }
496 else{
497 //Set PID bits for later selection (AliAnaPi0 for example)
498 //GetPDG already called in SetPIDBits.
499 GetCaloPID()->SetPIDBits(fCalorimeter,calo,&aodpi0);
500 if(GetDebug() > 1) printf("AliAnaPi0EbE::kSSCalo::::FillAOD: PID Bits set \n");
501 }
502
503 if(GetDebug() > 1) printf("AliAnaPi0EbE::kSSCalo::::FillAOD: Pi0 selection cuts passed: pT %3.2f, pdg %d\n",aodpi0.Pt(), aodpi0.GetPdg());
504
505 //Play with the MC stack if available
506 //Check origin of the candidates
507 if(IsDataMC()){
508 if((GetReader()->GetDataType() == AliCaloTrackReader::kMC && fAnaType!=kSSCalo) ||
509 GetReader()->GetDataType() != AliCaloTrackReader::kMC){
510 aodpi0.SetTag(GetMCAnalysisUtils()->CheckOrigin(calo->GetLabel(0),GetMCStack()));
511 if(GetDebug() > 0) printf("AliAnaPi0EbE::kSSCalo::EbE::FillAOD: Origin of candidate %d\n",aodpi0.GetTag());
512 }
513 }//Work with stack also
514
515 //Add AOD with pi0 object to aod branch
516 AddAODParticle(aodpi0);
517
518 }//loop
519
520 if(GetDebug() > 1) printf("AliAnaPi0EbE::kSSCalo::::FillAOD: End fill AODs \n");
521
522}
523//__________________________________________________________________
524void AliAnaPi0EbE::MakeAnalysisFillHistograms()
691bdd02 525{
477d6cee 526 //Do analysis and fill histograms
691bdd02 527
477d6cee 528 if(!GetOutputAODBranch()){
529 printf("AliAnaPi0EbE::FillHistos - No output pi0 in AOD branch with name < %s >,STOP \n",GetOutputAODName().Data());
530 abort();
531 }
532 //Loop on stored AOD pi0
533 Int_t naod = GetOutputAODBranch()->GetEntriesFast();
534 if(GetDebug() > 0) printf("AliAnaPi0EbE::Histo::pi0 aod branch entries %d\n", naod);
535
536 for(Int_t iaod = 0; iaod < naod ; iaod++){
537
538 AliAODPWG4Particle* pi0 = (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
539 Int_t pdg = pi0->GetPdg();
540
541 if(pdg != AliCaloPID::kPi0) continue;
542
543 //Fill pi0 histograms
544 Float_t pt = pi0->Pt();
545 Float_t phi = pi0->Phi();
546 Float_t eta = pi0->Eta();
547
548
549 fhPtPi0 ->Fill(pt);
550 fhPhiPi0 ->Fill(pt,phi);
551 fhEtaPi0 ->Fill(pt,eta);
552
553 if(IsDataMC()){
554 if((GetReader()->GetDataType() == AliCaloTrackReader::kMC && fAnaType!=kSSCalo) ||
555 GetReader()->GetDataType() != AliCaloTrackReader::kMC){
556 if(pi0->GetTag()== AliMCAnalysisUtils::kMCPi0){
557 fhPtMCPi0 ->Fill(pt);
558 fhPhiMCPi0 ->Fill(pt,phi);
559 fhEtaMCPi0 ->Fill(pt,eta);
560 }
561 else{
562 fhPtMCNoPi0 ->Fill(pt);
563 fhPhiMCNoPi0 ->Fill(pt,phi);
564 fhEtaMCNoPi0 ->Fill(pt,eta);
565 }
566 }
567 }//Histograms with MC
568
569 }// aod loop
570
571}
572
573
574//____________________________________________________________________________
575void AliAnaPi0EbE::Init()
576{
691bdd02 577 //Init
578 //Do some checks
579 if(fCalorimeter == "PHOS" && !GetReader()->IsPHOSSwitchedOn()){
580 printf("!!ABORT: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!\n");
581 abort();
582 }
583 else if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn()){
584 printf("!!ABORT: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!\n");
585 abort();
586 }
587
588}
589
477d6cee 590//____________________________________________________________________________
591void AliAnaPi0EbE::InitParameters()
592{
593 //Initialize the parameters of the analysis.
594 SetOutputAODClassName("AliAODPWG4Particle");
595 SetOutputAODName("pi0s");
596 fInputAODGammaConvName = "gammaconv" ;
597 fAnaType = kIMCalo ;
598 fCalorimeter = "EMCAL" ;
599 fMinDist = 2.;
600 fMinDist2 = 4.;
601 fMinDist3 = 5.;
602
603}
604
605//__________________________________________________________________
606void AliAnaPi0EbE::Print(const Option_t * opt) const
607{
608 //Print some relevant parameters set for the analysis
609 if(! opt)
610 return;
611
612 printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
613 AliAnaPartCorrBaseClass::Print("");
614 printf("Analysis Type = %d \n", fAnaType) ;
615 if(fAnaType == kSSCalo){
616 printf("Calorimeter = %s\n", fCalorimeter.Data()) ;
617 printf("Min Distance to Bad Channel = %2.1f\n",fMinDist);
618 printf("Min Distance to Bad Channel 2 = %2.1f\n",fMinDist2);
619 printf("Min Distance to Bad Channel 3 = %2.1f\n",fMinDist3);
620 }
621 printf(" \n") ;
622
623}