]> 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
array bound warnings are removed
[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>
0c1383b5 31#include <TObjString.h>
477d6cee 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"
477d6cee 42#include "AliStack.h"
ff45398a 43#include "AliFiducialCut.h"
477d6cee 44#include "TParticle.h"
0ae57829 45#include "AliVCluster.h"
477d6cee 46#include "AliAODEvent.h"
591cc579 47#include "AliAODMCParticle.h"
477d6cee 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}
78219bac 66/*
477d6cee 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}
78219bac 103*/
477d6cee 104//____________________________________________________________________________
105AliAnaPi0EbE::~AliAnaPi0EbE()
106{
107 //dtor
108 if(fInputAODGammaConv){
109 fInputAODGammaConv->Clear() ;
110 delete fInputAODGammaConv ;
111 }
112}
113
0c1383b5 114//________________________________________________________________________
115TObjString * AliAnaPi0EbE::GetAnalysisCuts()
116{
117 //Save parameters used for analysis
118 TString parList ; //this will be list of parameters used for this analysis.
5ae09196 119 const Int_t buffersize = 255;
120 char onePar[buffersize] ;
0c1383b5 121
5ae09196 122 snprintf(onePar,buffersize,"--- AliAnaPi0EbE ---\n") ;
0c1383b5 123 parList+=onePar ;
5ae09196 124 snprintf(onePar,buffersize,"fAnaType=%d (Pi0 selection type) \n",fAnaType) ;
0c1383b5 125 parList+=onePar ;
126
127 if(fAnaType == kSSCalo){
5ae09196 128 snprintf(onePar,buffersize,"Calorimeter: %s\n",fCalorimeter.Data()) ;
0c1383b5 129 parList+=onePar ;
5ae09196 130 snprintf(onePar,buffersize,"fMinDist =%2.2f (Minimal distance to bad channel to accept cluster) \n",fMinDist) ;
0c1383b5 131 parList+=onePar ;
5ae09196 132 snprintf(onePar,buffersize,"fMinDist2=%2.2f (Cuts on Minimal distance to study acceptance evaluation) \n",fMinDist2) ;
0c1383b5 133 parList+=onePar ;
5ae09196 134 snprintf(onePar,buffersize,"fMinDist3=%2.2f (One more cut on distance used for acceptance-efficiency study) \n",fMinDist3) ;
0c1383b5 135 parList+=onePar ;
136 }
137
138 //Get parameters set in base class.
139 parList += GetBaseParametersList() ;
140
141 //Get parameters set in PID class.
142 if(fAnaType == kSSCalo) parList += GetCaloPID()->GetPIDParametersList() ;
143
144 return new TObjString(parList) ;
145}
146
477d6cee 147//________________________________________________________________________
148TList * AliAnaPi0EbE::GetCreateOutputObjects()
149{
150 // Create histograms to be saved in output file and
151 // store them in outputContainer
152 TList * outputContainer = new TList() ;
153 outputContainer->SetName("Pi0EbEHistos") ;
154
5a2dbc3c 155 Int_t nptbins = GetHistoPtBins();
156 Int_t nphibins = GetHistoPhiBins();
157 Int_t netabins = GetHistoEtaBins();
477d6cee 158 Float_t ptmax = GetHistoPtMax();
159 Float_t phimax = GetHistoPhiMax();
160 Float_t etamax = GetHistoEtaMax();
161 Float_t ptmin = GetHistoPtMin();
162 Float_t phimin = GetHistoPhiMin();
163 Float_t etamin = GetHistoEtaMin();
164
165 fhPtPi0 = new TH1F("hPtPi0","Number of identified #pi^{0} decay",nptbins,ptmin,ptmax);
166 fhPtPi0->SetYTitle("N");
167 fhPtPi0->SetXTitle("p_{T #pi^{0}}(GeV/c)");
168 outputContainer->Add(fhPtPi0) ;
169
170 fhPhiPi0 = new TH2F
5ae09196 171 ("hPhiPi0","#phi_{#pi^{0}}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
477d6cee 172 fhPhiPi0->SetYTitle("#phi");
173 fhPhiPi0->SetXTitle("p_{T #pi^{0}} (GeV/c)");
174 outputContainer->Add(fhPhiPi0) ;
175
176 fhEtaPi0 = new TH2F
5ae09196 177 ("hEtaPi0","#phi_{#pi^{0}}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
477d6cee 178 fhEtaPi0->SetYTitle("#eta");
179 fhEtaPi0->SetXTitle("p_{T #pi^{0}} (GeV/c)");
180 outputContainer->Add(fhEtaPi0) ;
181
182 if(IsDataMC()) {
183 if((GetReader()->GetDataType() == AliCaloTrackReader::kMC && fAnaType!=kSSCalo) ||
184 GetReader()->GetDataType() != AliCaloTrackReader::kMC){
185
186 fhPtMCPi0 = new TH1F("hPtMCPi0","Identified pi0 from pi0",nptbins,ptmin,ptmax);
187 fhPtMCPi0->SetYTitle("N");
188 fhPtMCPi0->SetXTitle("p_{T #pi^{0}}(GeV/c)");
189 outputContainer->Add(fhPtMCPi0) ;
190
191 fhPhiMCPi0 = new TH2F
5ae09196 192 ("hPhiMCPi0","Identified pi0 from pi0",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
477d6cee 193 fhPhiMCPi0->SetYTitle("#phi");
194 fhPhiMCPi0->SetXTitle("p_{T #pi^{0}} (GeV/c)");
195 outputContainer->Add(fhPhiMCPi0) ;
196
197 fhEtaMCPi0 = new TH2F
5ae09196 198 ("hEtaMCPi0","Identified pi0 from pi0",nptbins,ptmin,ptmax,netabins,etamin,etamax);
477d6cee 199 fhEtaMCPi0->SetYTitle("#eta");
200 fhEtaMCPi0->SetXTitle("p_{T #pi^{0}} (GeV/c)");
201 outputContainer->Add(fhEtaMCPi0) ;
202
203 fhPtMCNoPi0 = new TH1F("hPtMCNoPi0","Identified pi0 not from pi0",nptbins,ptmin,ptmax);
204 fhPtMCNoPi0->SetYTitle("N");
205 fhPtMCNoPi0->SetXTitle("p_{T #pi^{0}}(GeV/c)");
206 outputContainer->Add(fhPtMCNoPi0) ;
207
208 fhPhiMCNoPi0 = new TH2F
5ae09196 209 ("hPhiMCNoPi0","Identified pi0 not from pi0",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
477d6cee 210 fhPhiMCNoPi0->SetYTitle("#phi");
211 fhPhiMCNoPi0->SetXTitle("p_{T #pi^{0}} (GeV/c)");
212 outputContainer->Add(fhPhiMCNoPi0) ;
213
214 fhEtaMCNoPi0 = new TH2F
5ae09196 215 ("hEtaMCNoPi0","Identified pi0 not from pi0",nptbins,ptmin,ptmax,netabins,etamin,etamax);
477d6cee 216 fhEtaMCNoPi0->SetYTitle("#eta");
217 fhEtaMCNoPi0->SetXTitle("p_{T #pi^{0}} (GeV/c)");
218 outputContainer->Add(fhEtaMCNoPi0) ;
219
220 }
221 }//Histos with MC
222
223
224 //Keep neutral meson selection histograms if requiered
225 //Setting done in AliNeutralMesonSelection
226
227 if(fAnaType!=kSSCalo && GetNeutralMesonSelection()){
228
229 TList * nmsHistos = GetNeutralMesonSelection()->GetCreateOutputObjects() ;
230 if(GetNeutralMesonSelection()->AreNeutralMesonSelectionHistosKept())
231 for(Int_t i = 0; i < nmsHistos->GetEntries(); i++) outputContainer->Add(nmsHistos->At(i)) ;
5ae09196 232 delete nmsHistos;
a14fd526 233
477d6cee 234 }
235
477d6cee 236 return outputContainer ;
237
238}
239
240//__________________________________________________________________
241void AliAnaPi0EbE::MakeAnalysisFillAOD()
242{
243 //Do analysis and fill aods
244
245 switch(fAnaType)
246 {
247 case kIMCalo:
248 MakeInvMassInCalorimeter();
249 break;
250
251 case kSSCalo:
252 MakeShowerShapeIdentification();
253 break;
254
255 case kIMCaloTracks:
256 MakeInvMassInCalorimeterAndCTS();
257 break;
258
259 }
260}
261
262//__________________________________________________________________
263void AliAnaPi0EbE::MakeInvMassInCalorimeter()
264{
c8fe2783 265 //Do analysis and fill aods
266 //Search for the photon decay in calorimeters
267 //Read photon list from AOD, produced in class AliAnaPhoton
268 //Check if 2 photons have the mass of the pi0.
477d6cee 269
270 TLorentzVector mom1;
271 TLorentzVector mom2;
272 TLorentzVector mom ;
591cc579 273 Int_t tag1 = 0;
274 Int_t tag2 = 0;
275 Int_t tag = 0;
477d6cee 276
277 if(!GetInputAODBranch()){
a3aebfff 278 printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - No input calo photons in AOD with name branch < %s >, STOP \n",GetInputAODName().Data());
477d6cee 279 abort();
280 }
281 for(Int_t iphoton = 0; iphoton < GetInputAODBranch()->GetEntriesFast(); iphoton++){
282 AliAODPWG4Particle * photon1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(iphoton));
c8fe2783 283
284 Int_t evtIndex1 = 0 ;
285 if(GetMixedEvent())
286 evtIndex1 = GetMixedEvent()->EventIndexForCaloCluster(photon1->GetCaloLabel(0)) ;
287
477d6cee 288 mom1 = *(photon1->Momentum());
289
290
291 for(Int_t jphoton = iphoton+1; jphoton < GetInputAODBranch()->GetEntriesFast()-1; jphoton++){
292
a3aebfff 293 AliAODPWG4Particle * photon2 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(jphoton));
c8fe2783 294 Int_t evtIndex2 = 0 ;
295 if(GetMixedEvent())
296 evtIndex2 = GetMixedEvent()->EventIndexForCaloCluster(photon2->GetCaloLabel(0)) ;
297 if(GetMixedEvent() && (evtIndex1 == evtIndex2))
298 continue ;
299
477d6cee 300 mom2 = *(photon2->Momentum());
c8fe2783 301 Int_t input = -1; //if -1 photons come from different files, not a pi0
302 if(photon1->GetInputFileIndex() == photon2->GetInputFileIndex())
303 input = photon1->GetInputFileIndex();
304
305 //Select good pair (good phi, pt cuts, aperture and invariant mass)
477d6cee 306 if(GetNeutralMesonSelection()->SelectPair(mom1, mom2))
c8fe2783 307 {
308 if(GetDebug()>1)
309 printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - Selected gamma pair: pt %f, phi %f, eta%f \n",(mom1+mom2).Pt(), (mom1+mom2).Phi()*180./3.1416, (mom1+mom2).Eta());
310
311 //Play with the MC stack if available
312 if(IsDataMC()){
313 //Check origin of the candidates
314 Int_t label1 = photon1->GetLabel();
315 Int_t label2 = photon2->GetLabel();
316 tag1 = GetMCAnalysisUtils()->CheckOrigin(label1, GetReader(), photon1->GetInputFileIndex());
317 tag2 = GetMCAnalysisUtils()->CheckOrigin(label2, GetReader(), photon2->GetInputFileIndex());
318
319 if(GetDebug() > 0) printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - Origin of: photon1 %d; photon2 %d \n",tag1, tag2);
320 if(GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCPi0Decay) && GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCPi0Decay)){
321
322 //Check if pi0 mother is the same
323 if(GetReader()->ReadStack()){
324 TParticle * mother1 = GetMCStack()->Particle(label1);//photon in kine tree
325 label1 = mother1->GetFirstMother();
326 //mother1 = GetMCStack()->Particle(label1);//pi0
327
328 TParticle * mother2 = GetMCStack()->Particle(label2);//photon in kine tree
329 label2 = mother2->GetFirstMother();
330 //mother2 = GetMCStack()->Particle(label2);//pi0
331 }
332 else if(GetReader()->ReadAODMCParticles() && (input > -1)){
333 AliAODMCParticle * mother1 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles(photon1->GetInputFileIndex()))->At(label1);//photon in kine tree
334 label1 = mother1->GetMother();
335 //mother1 = GetMCStack()->Particle(label1);//pi0
336 AliAODMCParticle * mother2 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles(photon2->GetInputFileIndex()))->At(label2);//photon in kine tree
337 label2 = mother2->GetMother();
338 //mother2 = GetMCStack()->Particle(label2);//pi0
339 }
340
341 //printf("mother1 %d, mother2 %d\n",label1,label2);
342 if(label1 == label2)
343 GetMCAnalysisUtils()->SetTagBit(tag,AliMCAnalysisUtils::kMCPi0);
344 }
345 }//Work with stack also
346
347 //Create AOD for analysis
348 mom = mom1+mom2;
349 AliAODPWG4Particle pi0 = AliAODPWG4Particle(mom);
350 //pi0.SetLabel(calo->GetLabel());
351 pi0.SetPdg(AliCaloPID::kPi0);
352 pi0.SetDetector(photon1->GetDetector());
353 pi0.SetTag(tag);
354 //Set the indeces of the original caloclusters
355 pi0.SetCaloLabel(photon1->GetCaloLabel(0), photon2->GetCaloLabel(0));
356 pi0.SetInputFileIndex(input);
357 AddAODParticle(pi0);
358 }//pi0
477d6cee 359 }//2n photon loop
360
361 }//1st photon loop
362
a3aebfff 363 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - End fill AODs \n");
477d6cee 364
365}
366
367//__________________________________________________________________
368void AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS()
369{
370 //Do analysis and fill aods
371 //Search for the photon decay in calorimeters
372 //Read photon list from AOD, produced in class AliAnaPhoton and AliGammaConversion
373 //Check if 2 photons have the mass of the pi0.
374
375 TLorentzVector mom1;
376 TLorentzVector mom2;
377 TLorentzVector mom ;
591cc579 378 Int_t tag1 = 0;
379 Int_t tag2 = 0;
380 Int_t tag = 0;
477d6cee 381
382 if(!GetInputAODBranch()){
a3aebfff 383 printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - No input calo photons in AOD branch with name < %s > , STOP\n",GetInputAODName().Data());
477d6cee 384 abort();
385 }
386 for(Int_t iphoton = 0; iphoton < GetInputAODBranch()->GetEntriesFast(); iphoton++){
387 AliAODPWG4Particle * photon1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(iphoton));
388 mom1 = *(photon1->Momentum());
389
390 //Play with the MC stack if available
391 fInputAODGammaConv = (TClonesArray *) GetReader()->GetOutputEvent()->FindListObject(fInputAODGammaConvName);
392 if(!fInputAODGammaConv) {
a3aebfff 393 printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - No input gamma conversions in AOD branch with name < %s >, STOP\n",fInputAODGammaConvName.Data());
477d6cee 394 abort();
395 }
396 for(Int_t jphoton = iphoton+1; jphoton < fInputAODGammaConv->GetEntriesFast()-1; jphoton++){
a3aebfff 397 AliAODPWG4Particle * photon2 = (AliAODPWG4Particle*) (fInputAODGammaConv->At(jphoton));
477d6cee 398 mom2 = *(photon2->Momentum());
7cd4e982 399
400 Int_t input = -1; //if -1 photons come from different files, not a pi0
401 if(photon1->GetInputFileIndex() == photon2->GetInputFileIndex()) input = photon1->GetInputFileIndex();
402
477d6cee 403 //Select good pair (good phi, pt cuts, aperture and invariant mass)
404 if(GetNeutralMesonSelection()->SelectPair(mom1, mom2)){
a3aebfff 405 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - Selected gamma pair: pt %f, phi %f, eta%f\n",(mom1+mom2).Pt(), (mom1+mom2).Phi()*180./3.1416, (mom1+mom2).Eta());
477d6cee 406
407 if(IsDataMC()){
7cd4e982 408 Int_t label1 = photon1->GetLabel();
409 Int_t label2 = photon2->GetLabel();
410 tag1 = GetMCAnalysisUtils()->CheckOrigin(label1, GetReader(), photon1->GetInputFileIndex());
411 tag2 = GetMCAnalysisUtils()->CheckOrigin(label2, GetReader(), photon2->GetInputFileIndex());
a3aebfff 412 if(GetDebug() > 0) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - Origin of: photon1 %d; photon2 %d \n",tag1, tag2);
591cc579 413 if(GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCPi0Decay) && GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCPi0Decay)){
477d6cee 414 //Check if pi0 mother is the same
591cc579 415
416 if(GetReader()->ReadStack()){
417 TParticle * mother1 = GetMCStack()->Particle(label1);//photon in kine tree
418 label1 = mother1->GetFirstMother();
419 //mother1 = GetMCStack()->Particle(label1);//pi0
477d6cee 420
591cc579 421 TParticle * mother2 = GetMCStack()->Particle(label2);//photon in kine tree
422 label2 = mother2->GetFirstMother();
423 //mother2 = GetMCStack()->Particle(label2);//pi0
424 }
7cd4e982 425 else if(GetReader()->ReadAODMCParticles() && (input > -1)){
426 AliAODMCParticle * mother1 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles(photon1->GetInputFileIndex()))->At(label1);//photon in kine tree
591cc579 427 label1 = mother1->GetMother();
428 //mother1 = GetMCStack()->Particle(label1);//pi0
7cd4e982 429 AliAODMCParticle * mother2 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles(photon2->GetInputFileIndex()))->At(label2);//photon in kine tree
591cc579 430 label2 = mother2->GetMother();
431 //mother2 = GetMCStack()->Particle(label2);//pi0
432 }
433
434 //printf("mother1 %d, mother2 %d\n",label1,label2);
477d6cee 435 if(label1 == label2)
591cc579 436 GetMCAnalysisUtils()->SetTagBit(tag,AliMCAnalysisUtils::kMCPi0);
477d6cee 437 }
438 }//Work with stack also
439
440 //Create AOD for analysis
441 mom = mom1+mom2;
a3aebfff 442 AliAODPWG4Particle pi0 = AliAODPWG4Particle(mom);
c8fe2783 443 //pi0.SetLabel(calo->GetLabel());
477d6cee 444 pi0.SetPdg(AliCaloPID::kPi0);
445 pi0.SetDetector(photon1->GetDetector());
446 pi0.SetTag(tag);
447 //Set the indeces of the original tracks or caloclusters
448 pi0.SetCaloLabel(photon1->GetCaloLabel(0), -1);
449 pi0.SetTrackLabel(photon2->GetTrackLabel(0), photon2->GetTrackLabel(1));
7cd4e982 450 pi0.SetInputFileIndex(input);
477d6cee 451 AddAODParticle(pi0);
452 }//pi0
453 }//2n photon loop
454
455 }//1st photon loop
456
a3aebfff 457 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - End fill AODs \n");
477d6cee 458
459}
460
461
462//__________________________________________________________________
463void AliAnaPi0EbE::MakeShowerShapeIdentification()
464{
465 //Search for pi0 in fCalorimeter with shower shape analysis
466
4a745797 467 TObjArray * pl = 0x0;
5ae09196 468 //Select the Calorimeter of the photon
469 if(fCalorimeter == "PHOS")
470 pl = GetAODPHOS();
471 else if (fCalorimeter == "EMCAL")
472 pl = GetAODEMCAL();
473
474 if(!pl) {
475 Info("MakeShowerShapeIdentification","TObjArray with %s clusters is NULL!\n",fCalorimeter.Data());
476 return;
477 }
477d6cee 478
479 //Get vertex for photon momentum calculation
233e0df8 480 Double_t vertex[] = {0,0,0} ; //vertex
f37fa8d2 481 //Double_t vertex2[] = {0,0,0} ; //vertex from second aod input
edd59991 482 if(GetReader()->GetDataType() != AliCaloTrackReader::kMC)
233e0df8 483 {
484 GetReader()->GetVertex(vertex);
1e68a3f4 485 //if(GetReader()->GetSecondInputAODTree()) GetReader()->GetSecondInputAODVertex(vertex2);
233e0df8 486 }
487
0ae57829 488
477d6cee 489 TLorentzVector mom ;
490 for(Int_t icalo = 0; icalo < pl->GetEntriesFast(); icalo++){
0ae57829 491 AliVCluster * calo = (AliVCluster*) (pl->At(icalo));
477d6cee 492
ff45398a 493 //Cluster selection, not charged, with pi0 id and in fiducial cut
233e0df8 494
495 //Input from second AOD?
496 Int_t input = 0;
1e68a3f4 497// if (fCalorimeter == "EMCAL" && GetReader()->GetAODEMCALNormalInputEntries() <= icalo) input = 1 ;
498// else if(fCalorimeter == "PHOS" && GetReader()->GetAODPHOSNormalInputEntries() <= icalo) input = 1;
233e0df8 499
500 //Get Momentum vector,
501 if (input == 0) calo->GetMomentum(mom,vertex) ;//Assume that come from vertex in straight line
1e68a3f4 502 //else if(input == 1) calo->GetMomentum(mom,vertex2);//Assume that come from vertex in straight line
233e0df8 503
504 //If too small or big pt, skip it
477d6cee 505 if(mom.Pt() < GetMinPt() || mom.Pt() > GetMaxPt() ) continue ;
506 //Check acceptance selection
ff45398a 507 if(IsFiducialCutOn()){
508 Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter) ;
477d6cee 509 if(! in ) continue ;
510 }
511
512 //Create AOD for analysis
513 AliAODPWG4Particle aodpi0 = AliAODPWG4Particle(mom);
c8fe2783 514 aodpi0.SetLabel(calo->GetLabel());
477d6cee 515 //Set the indeces of the original caloclusters
516 aodpi0.SetCaloLabel(calo->GetID(),-1);
517 aodpi0.SetDetector(fCalorimeter);
518 if(GetDebug() > 1)
ff45398a 519 printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - FillAOD: Min pt cut and fiducial cut passed: pt %3.2f, phi %2.2f, eta %1.2f\n",aodpi0.Pt(),aodpi0.Phi(),aodpi0.Eta());
477d6cee 520
521 //Check Distance to Bad channel, set bit.
c8fe2783 522 Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
477d6cee 523 if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;
524 if(distBad < fMinDist) //In bad channel (PHOS cristal size 2.2x2.2 cm)
525 continue ;
526
a3aebfff 527 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - FillAOD: Bad channel cut passed %4.2f\n",distBad);
477d6cee 528
529 if(distBad > fMinDist3) aodpi0.SetDistToBad(2) ;
530 else if(distBad > fMinDist2) aodpi0.SetDistToBad(1) ;
531 else aodpi0.SetDistToBad(0) ;
532
533 //Check PID
534 //PID selection or bit setting
535 if(GetReader()->GetDataType() == AliCaloTrackReader::kMC){
536 //Get most probable PID, check PID weights (in MC this option is mandatory)
c8fe2783 537 aodpi0.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,calo->GetPID(),mom.E()));//PID with weights
477d6cee 538 if(GetDebug() > 1)
a3aebfff 539 printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - FillAOD: PDG of identified particle %d\n",aodpi0.GetPdg());
477d6cee 540 //If primary is not pi0, skip it.
541 if(aodpi0.GetPdg() != AliCaloPID::kPi0) continue ;
542 }
543 else if(IsCaloPIDOn()){
544 //Skip matched clusters with tracks
f37fa8d2 545 if(IsTrackMatched(calo)) continue ;
477d6cee 546
547 //Get most probable PID, 2 options check PID weights
548 //or redo PID, recommended option for EMCal.
549 if(!IsCaloPIDRecalculationOn())
c8fe2783 550 aodpi0.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,calo->GetPID(),mom.E()));//PID with weights
477d6cee 551 else
552 aodpi0.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,mom,calo));//PID recalculated
553
a3aebfff 554 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - PDG of identified particle %d\n",aodpi0.GetPdg());
477d6cee 555
556 //If cluster does not pass pid, not pi0, skip it.
557 if(aodpi0.GetPdg() != AliCaloPID::kPi0) continue ;
558
559 }
560 else{
561 //Set PID bits for later selection (AliAnaPi0 for example)
562 //GetPDG already called in SetPIDBits.
563 GetCaloPID()->SetPIDBits(fCalorimeter,calo,&aodpi0);
a3aebfff 564 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - PID Bits set \n");
477d6cee 565 }
566
a3aebfff 567 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - Pi0 selection cuts passed: pT %3.2f, pdg %d\n",aodpi0.Pt(), aodpi0.GetPdg());
477d6cee 568
569 //Play with the MC stack if available
570 //Check origin of the candidates
571 if(IsDataMC()){
572 if((GetReader()->GetDataType() == AliCaloTrackReader::kMC && fAnaType!=kSSCalo) ||
233e0df8 573 GetReader()->GetDataType() != AliCaloTrackReader::kMC){
574 aodpi0.SetInputFileIndex(input);
8dc8c40b 575 Int_t tag =0;
c8fe2783 576 tag = GetMCAnalysisUtils()->CheckOrigin(calo->GetLabel(),GetReader(), aodpi0.GetInputFileIndex());
577 //GetMCAnalysisUtils()->CheckMultipleOrigin(calo->GetLabels(),calo->GetNLabels(), GetReader(), aodpi0.GetInputFileIndex(), tag);
8dc8c40b 578 aodpi0.SetTag(tag);
233e0df8 579 if(GetDebug() > 0) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - Origin of candidate %d\n",aodpi0.GetTag());
477d6cee 580 }
581 }//Work with stack also
582
583 //Add AOD with pi0 object to aod branch
584 AddAODParticle(aodpi0);
585
586 }//loop
587
a3aebfff 588 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - End fill AODs \n");
477d6cee 589
590}
591//__________________________________________________________________
592void AliAnaPi0EbE::MakeAnalysisFillHistograms()
691bdd02 593{
477d6cee 594 //Do analysis and fill histograms
691bdd02 595
477d6cee 596 if(!GetOutputAODBranch()){
a3aebfff 597 printf("AliAnaPi0EbE::MakeAnalysisFillHistograms() - No output pi0 in AOD branch with name < %s >,STOP \n",GetOutputAODName().Data());
477d6cee 598 abort();
599 }
600 //Loop on stored AOD pi0
601 Int_t naod = GetOutputAODBranch()->GetEntriesFast();
a3aebfff 602 if(GetDebug() > 0) printf("AliAnaPi0EbE::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod);
477d6cee 603
604 for(Int_t iaod = 0; iaod < naod ; iaod++){
605
606 AliAODPWG4Particle* pi0 = (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
607 Int_t pdg = pi0->GetPdg();
9415d854 608
609 if(IsCaloPIDOn() && pdg != AliCaloPID::kPi0) continue;
477d6cee 610
611 //Fill pi0 histograms
7cd4e982 612 Float_t pt = pi0->Pt();
477d6cee 613 Float_t phi = pi0->Phi();
614 Float_t eta = pi0->Eta();
615
477d6cee 616 fhPtPi0 ->Fill(pt);
617 fhPhiPi0 ->Fill(pt,phi);
618 fhEtaPi0 ->Fill(pt,eta);
619
620 if(IsDataMC()){
621 if((GetReader()->GetDataType() == AliCaloTrackReader::kMC && fAnaType!=kSSCalo) ||
622 GetReader()->GetDataType() != AliCaloTrackReader::kMC){
591cc579 623 if(GetMCAnalysisUtils()->CheckTagBit(pi0->GetTag(), AliMCAnalysisUtils::kMCPi0)){
477d6cee 624 fhPtMCPi0 ->Fill(pt);
625 fhPhiMCPi0 ->Fill(pt,phi);
626 fhEtaMCPi0 ->Fill(pt,eta);
627 }
628 else{
629 fhPtMCNoPi0 ->Fill(pt);
630 fhPhiMCNoPi0 ->Fill(pt,phi);
631 fhEtaMCNoPi0 ->Fill(pt,eta);
632 }
633 }
634 }//Histograms with MC
635
636 }// aod loop
637
638}
639
640
641//____________________________________________________________________________
642void AliAnaPi0EbE::Init()
643{
691bdd02 644 //Init
645 //Do some checks
1e86c71e 646 if(fCalorimeter == "PHOS" && !GetReader()->IsPHOSSwitchedOn() && NewOutputAOD()){
591cc579 647 printf("AliAnaPi0EbE::Init() - !!STOP: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!\n");
691bdd02 648 abort();
649 }
1e86c71e 650 else if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn() && NewOutputAOD()){
591cc579 651 printf("AliAnaPi0EbE::Init() - !!STOP: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!\n");
691bdd02 652 abort();
653 }
654
655}
656
477d6cee 657//____________________________________________________________________________
658void AliAnaPi0EbE::InitParameters()
659{
1e86c71e 660 //Initialize the parameters of the analysis.
a3aebfff 661 AddToHistogramsName("AnaPi0EbE_");
662
477d6cee 663 fInputAODGammaConvName = "gammaconv" ;
664 fAnaType = kIMCalo ;
665 fCalorimeter = "EMCAL" ;
666 fMinDist = 2.;
667 fMinDist2 = 4.;
668 fMinDist3 = 5.;
669
670}
671
672//__________________________________________________________________
673void AliAnaPi0EbE::Print(const Option_t * opt) const
674{
675 //Print some relevant parameters set for the analysis
676 if(! opt)
677 return;
678
679 printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
680 AliAnaPartCorrBaseClass::Print("");
681 printf("Analysis Type = %d \n", fAnaType) ;
682 if(fAnaType == kSSCalo){
683 printf("Calorimeter = %s\n", fCalorimeter.Data()) ;
684 printf("Min Distance to Bad Channel = %2.1f\n",fMinDist);
685 printf("Min Distance to Bad Channel 2 = %2.1f\n",fMinDist2);
686 printf("Min Distance to Bad Channel 3 = %2.1f\n",fMinDist3);
687 }
688 printf(" \n") ;
689
690}