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