1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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: AliAnaPi0EbE.cxx 28688 2008-09-11 15:04:07Z gconesab $ */
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 CTS
24 // -- Author: Gustavo Conesa (LNF-INFN) & Raphaelle Ichou (SUBATECH)
25 //////////////////////////////////////////////////////////////////////////////
28 // --- ROOT system ---
30 #include <TClonesArray.h>
31 #include <TObjString.h>
33 //#include "Riostream.h"
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"
43 #include "AliFiducialCut.h"
44 #include "TParticle.h"
45 #include "AliVCluster.h"
46 #include "AliAODEvent.h"
47 #include "AliAODMCParticle.h"
49 ClassImp(AliAnaPi0EbE)
51 //____________________________________________________________________________
52 AliAnaPi0EbE::AliAnaPi0EbE() :
53 AliAnaPartCorrBaseClass(), fAnaType(kIMCalo), fCalorimeter(""),
54 fMinDist(0.),fMinDist2(0.), fMinDist3(0.),
55 fInputAODGammaConv(0x0), fInputAODGammaConvName(""),
57 fhPtPi0(0), fhEPi0(0), fhEEtaPhiPi0(0),
59 fhEDispersion(0), fhELambda0(0), fhELambda1(0),
61 fhClusterPairDiffTimeE(0), fhClusterPairDiffTimeAsy(0),
63 fhPtMCNoPi0(0), fhPhiMCNoPi0(0), fhEtaMCNoPi0(0),
64 fhPtMCPi0(0), fhPhiMCPi0(0), fhEtaMCPi0(0)
68 for(Int_t i = 0; i < 5; i++){
71 fhEMCDispersion[i] = 0;
74 //Initialize parameters
79 //____________________________________________________________________________
80 AliAnaPi0EbE::~AliAnaPi0EbE()
83 if(fInputAODGammaConv){
84 fInputAODGammaConv->Clear() ;
85 delete fInputAODGammaConv ;
89 //________________________________________________________________________
90 TObjString * AliAnaPi0EbE::GetAnalysisCuts()
92 //Save parameters used for analysis
93 TString parList ; //this will be list of parameters used for this analysis.
94 const Int_t buffersize = 255;
95 char onePar[buffersize] ;
97 snprintf(onePar,buffersize,"--- AliAnaPi0EbE ---\n") ;
99 snprintf(onePar,buffersize,"fAnaType=%d (Pi0 selection type) \n",fAnaType) ;
102 if(fAnaType == kSSCalo){
103 snprintf(onePar,buffersize,"Calorimeter: %s\n",fCalorimeter.Data()) ;
105 snprintf(onePar,buffersize,"fMinDist =%2.2f (Minimal distance to bad channel to accept cluster) \n",fMinDist) ;
107 snprintf(onePar,buffersize,"fMinDist2=%2.2f (Cuts on Minimal distance to study acceptance evaluation) \n",fMinDist2) ;
109 snprintf(onePar,buffersize,"fMinDist3=%2.2f (One more cut on distance used for acceptance-efficiency study) \n",fMinDist3) ;
113 //Get parameters set in base class.
114 parList += GetBaseParametersList() ;
116 //Get parameters set in PID class.
117 if(fAnaType == kSSCalo) parList += GetCaloPID()->GetPIDParametersList() ;
119 return new TObjString(parList) ;
122 //________________________________________________________________________
123 TList * AliAnaPi0EbE::GetCreateOutputObjects()
125 // Create histograms to be saved in output file and
126 // store them in outputContainer
127 TList * outputContainer = new TList() ;
128 outputContainer->SetName("Pi0EbEHistos") ;
130 Int_t nptbins = GetHistoPtBins(); Float_t ptmax = GetHistoPtMax(); Float_t ptmin = GetHistoPtMin();
131 Int_t nphibins = GetHistoPhiBins(); Float_t phimax = GetHistoPhiMax(); Float_t phimin = GetHistoPhiMin();
132 Int_t netabins = GetHistoEtaBins(); Float_t etamax = GetHistoEtaMax(); Float_t etamin = GetHistoEtaMin();
133 Int_t ssbins = GetHistoShowerShapeBins(); Float_t ssmax = GetHistoShowerShapeMax(); Float_t ssmin = GetHistoShowerShapeMin();
134 Int_t tdbins = GetHistoDiffTimeBins() ; Float_t tdmax = GetHistoDiffTimeMax(); Float_t tdmin = GetHistoDiffTimeMin();
136 fhPtPi0 = new TH1F("hPtPi0","Number of identified #pi^{0} decay",nptbins,ptmin,ptmax);
137 fhPtPi0->SetYTitle("N");
138 fhPtPi0->SetXTitle("p_{T #pi^{0}}(GeV/c)");
139 outputContainer->Add(fhPtPi0) ;
141 fhEPi0 = new TH1F("hEPi0","Number of identified #pi^{0} decay",nptbins,ptmin,ptmax);
142 fhEPi0->SetYTitle("N");
143 fhEPi0->SetXTitle("E #pi^{0}(GeV)");
144 outputContainer->Add(fhEPi0) ;
146 fhEEtaPhiPi0 = new TH3F
147 ("hEEtaPhiPi0","Selected #pi^{0} pairs: E vs #eta vs #phi",nptbins,ptmin,ptmax,netabins,etamin,etamax, nphibins,phimin,phimax);
148 fhEEtaPhiPi0->SetZTitle("#phi");
149 fhEEtaPhiPi0->SetYTitle("#eta");
150 fhEEtaPhiPi0->SetXTitle("E (GeV)");
151 outputContainer->Add(fhEEtaPhiPi0) ;
155 if(fAnaType == kIMCalo){
157 fhEDispersion = new TH2F
158 ("hEDispersion","Selected #pi^{0} pairs: E vs dispersion",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
159 fhEDispersion->SetYTitle("D^{2}");
160 fhEDispersion->SetXTitle("E (GeV)");
161 outputContainer->Add(fhEDispersion) ;
163 fhELambda0 = new TH2F
164 ("hELambda0","Selected #pi^{0} pairs: E vs #lambda_{0}",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
165 fhELambda0->SetYTitle("#lambda_{0}^{2}");
166 fhELambda0->SetXTitle("E (GeV)");
167 outputContainer->Add(fhELambda0) ;
169 fhELambda1 = new TH2F
170 ("hELambda1","Selected #pi^{0} pairs: E vs #lambda_{1}",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
171 fhELambda1->SetYTitle("#lambda_{1}^{2}");
172 fhELambda1->SetXTitle("E (GeV)");
173 outputContainer->Add(fhELambda1) ;
175 fhClusterPairDiffTimeE = new TH2F("hClusterPairDiffTimeE","cluster pair time difference vs E",nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
176 fhClusterPairDiffTimeE->SetXTitle("E_{pair} (GeV)");
177 fhClusterPairDiffTimeE->SetYTitle("#Delta t (ns)");
178 outputContainer->Add(fhClusterPairDiffTimeE);
180 fhClusterPairDiffTimeAsy = new TH2F("hClusterPairDiffTime","cluster pair time difference vs pair asymmetry",100,0,1, tdbins,tdmin,tdmax);
181 fhClusterPairDiffTimeAsy->SetXTitle("Asymmetry");
182 fhClusterPairDiffTimeAsy->SetYTitle("#Delta t (ns)");
183 outputContainer->Add(fhClusterPairDiffTimeAsy);
185 }// Invariant mass analysis in calorimeters only
188 if((GetReader()->GetDataType() == AliCaloTrackReader::kMC && fAnaType!=kSSCalo) ||
189 GetReader()->GetDataType() != AliCaloTrackReader::kMC){
191 fhPtMCPi0 = new TH1F("hPtMCPi0","Identified pi0 from pi0",nptbins,ptmin,ptmax);
192 fhPtMCPi0->SetYTitle("N");
193 fhPtMCPi0->SetXTitle("p_{T #pi^{0}}(GeV/c)");
194 outputContainer->Add(fhPtMCPi0) ;
196 fhPhiMCPi0 = new TH2F
197 ("hPhiMCPi0","Identified pi0 from pi0",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
198 fhPhiMCPi0->SetYTitle("#phi");
199 fhPhiMCPi0->SetXTitle("p_{T #pi^{0}} (GeV/c)");
200 outputContainer->Add(fhPhiMCPi0) ;
202 fhEtaMCPi0 = new TH2F
203 ("hEtaMCPi0","Identified pi0 from pi0",nptbins,ptmin,ptmax,netabins,etamin,etamax);
204 fhEtaMCPi0->SetYTitle("#eta");
205 fhEtaMCPi0->SetXTitle("p_{T #pi^{0}} (GeV/c)");
206 outputContainer->Add(fhEtaMCPi0) ;
208 fhPtMCNoPi0 = new TH1F("hPtMCNoPi0","Identified pi0 not from pi0",nptbins,ptmin,ptmax);
209 fhPtMCNoPi0->SetYTitle("N");
210 fhPtMCNoPi0->SetXTitle("p_{T #pi^{0}}(GeV/c)");
211 outputContainer->Add(fhPtMCNoPi0) ;
213 fhPhiMCNoPi0 = new TH2F
214 ("hPhiMCNoPi0","Identified pi0 not from pi0",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
215 fhPhiMCNoPi0->SetYTitle("#phi");
216 fhPhiMCNoPi0->SetXTitle("p_{T #pi^{0}} (GeV/c)");
217 outputContainer->Add(fhPhiMCNoPi0) ;
219 fhEtaMCNoPi0 = new TH2F
220 ("hEtaMCNoPi0","Identified pi0 not from pi0",nptbins,ptmin,ptmax,netabins,etamin,etamax);
221 fhEtaMCNoPi0->SetYTitle("#eta");
222 fhEtaMCNoPi0->SetXTitle("p_{T #pi^{0}} (GeV/c)");
223 outputContainer->Add(fhEtaMCNoPi0) ;
225 if(fAnaType == kIMCalo){
226 TString ptype[] ={"#gamma","#gamma->e^{#pm}","#pi^{0}","e^{#pm}", "hadron"};
227 TString pname[] ={"Photon","Conversion", "Pi0", "Electron","Hadron"};
228 for(Int_t i = 0; i < 5; i++){
230 fhEMCLambda0[i] = new TH2F(Form("hELambda0_MC%s",pname[i].Data()),
231 Form("Selected pair, cluster from %s : E vs #lambda_{0}^{2}",ptype[i].Data()),
232 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
233 fhEMCLambda0[i]->SetYTitle("#lambda_{0}^{2}");
234 fhEMCLambda0[i]->SetXTitle("E (GeV)");
235 outputContainer->Add(fhEMCLambda0[i]) ;
237 fhEMCLambda1[i] = new TH2F(Form("hELambda1_MC%s",pname[i].Data()),
238 Form("Selected pair, cluster from %s : E vs #lambda_{1}^{2}",ptype[i].Data()),
239 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
240 fhEMCLambda1[i]->SetYTitle("#lambda_{1}^{2}");
241 fhEMCLambda1[i]->SetXTitle("E (GeV)");
242 outputContainer->Add(fhEMCLambda1[i]) ;
244 fhEMCDispersion[i] = new TH2F(Form("hEDispersion_MC%s",pname[i].Data()),
245 Form("Selected pair, cluster from %s : E vs dispersion^{2}",ptype[i].Data()),
246 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
247 fhEMCDispersion[i]->SetYTitle("D^{2}");
248 fhEMCDispersion[i]->SetXTitle("E (GeV)");
249 outputContainer->Add(fhEMCDispersion[i]) ;
258 //Keep neutral meson selection histograms if requiered
259 //Setting done in AliNeutralMesonSelection
261 if(fAnaType!=kSSCalo && GetNeutralMesonSelection()){
263 TList * nmsHistos = GetNeutralMesonSelection()->GetCreateOutputObjects() ;
264 if(GetNeutralMesonSelection()->AreNeutralMesonSelectionHistosKept())
265 for(Int_t i = 0; i < nmsHistos->GetEntries(); i++) outputContainer->Add(nmsHistos->At(i)) ;
270 return outputContainer ;
274 //____________________________________________________________________________
275 void AliAnaPi0EbE::Init()
279 if(fCalorimeter == "PHOS" && !GetReader()->IsPHOSSwitchedOn() && NewOutputAOD()){
280 printf("AliAnaPi0EbE::Init() - !!STOP: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!\n");
283 else if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn() && NewOutputAOD()){
284 printf("AliAnaPi0EbE::Init() - !!STOP: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!\n");
290 //____________________________________________________________________________
291 void AliAnaPi0EbE::InitParameters()
293 //Initialize the parameters of the analysis.
294 AddToHistogramsName("AnaPi0EbE_");
296 fInputAODGammaConvName = "gammaconv" ;
298 fCalorimeter = "EMCAL" ;
305 //__________________________________________________________________
306 void AliAnaPi0EbE::MakeAnalysisFillAOD()
308 //Do analysis and fill aods
313 MakeInvMassInCalorimeter();
317 MakeShowerShapeIdentification();
321 MakeInvMassInCalorimeterAndCTS();
327 //__________________________________________________________________
328 void AliAnaPi0EbE::MakeInvMassInCalorimeter()
330 //Do analysis and fill aods
331 //Search for the photon decay in calorimeters
332 //Read photon list from AOD, produced in class AliAnaPhoton
333 //Check if 2 photons have the mass of the pi0.
342 if(!GetInputAODBranch()){
343 printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - No input calo photons in AOD with name branch < %s >, STOP \n",GetInputAODName().Data());
347 for(Int_t iphoton = 0; iphoton < GetInputAODBranch()->GetEntriesFast()-1; iphoton++){
348 AliAODPWG4Particle * photon1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(iphoton));
350 //Vertex cut in case of mixed events
351 Int_t evtIndex1 = 0 ;
353 evtIndex1 = GetMixedEvent()->EventIndexForCaloCluster(photon1->GetCaloLabel(0)) ;
354 if(TMath::Abs(GetVertex(evtIndex1)[2]) > GetZvertexCut()) continue ; //vertex cut
355 mom1 = *(photon1->Momentum());
357 //Get shower shape information of clusters
358 TObjArray *clusters = 0;
359 if (!strcmp(fCalorimeter,"EMCAL")) clusters = GetEMCALClusters();
360 else if(!strcmp(fCalorimeter,"PHOS")) clusters = GetPHOSClusters() ;
362 Bool_t bFound1 = kFALSE;
363 Int_t caloLabel1 = photon1->GetCaloLabel(0);
365 AliVCluster *cluster1 = 0;
367 //Get original cluster, to recover some information
368 for(Int_t iclus = 0; iclus < clusters->GetEntriesFast(); iclus++){
369 AliVCluster *cluster= dynamic_cast<AliVCluster*> (clusters->At(iclus));
371 if (cluster->GetID()==caloLabel1) {
378 }// calorimeter clusters loop
381 for(Int_t jphoton = iphoton+1; jphoton < GetInputAODBranch()->GetEntriesFast(); jphoton++){
383 AliAODPWG4Particle * photon2 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(jphoton));
384 Int_t evtIndex2 = 0 ;
386 evtIndex2 = GetMixedEvent()->EventIndexForCaloCluster(photon2->GetCaloLabel(0)) ;
387 if(GetMixedEvent() && (evtIndex1 == evtIndex2))
389 if(TMath::Abs(GetVertex(evtIndex2)[2]) > GetZvertexCut()) continue ; //vertex cut
390 mom2 = *(photon2->Momentum());
393 Float_t e1 = photon1->E();
400 Float_t e2 = photon2->E();
405 Bool_t bFound2 = kFALSE;
406 Int_t caloLabel2 = photon2->GetCaloLabel(0);
407 AliVCluster *cluster2 = 0;
409 for(Int_t iclus = iclus1+1; iclus < clusters->GetEntriesFast(); iclus++){
410 AliVCluster *cluster= dynamic_cast<AliVCluster*> (clusters->At(iclus));
412 if(cluster->GetID()==caloLabel2) {
418 }// calorimeter clusters loop
421 if(cluster1 && bFound1){
422 disp1 = cluster1->GetDispersion()*cluster1->GetDispersion();
423 l11 = cluster1->GetM20();
424 l01 = cluster1->GetM02();
425 tof1 = cluster1->GetTOF()*1e9;
427 // else printf("cluster1 not available: calo label %d / %d, cluster ID %d\n",
428 // photon2->GetCaloLabel(0),(GetReader()->GetInputEvent())->GetNumberOfCaloClusters()-1,cluster1->GetID());
431 if(cluster2 && bFound2){
432 disp2 = cluster2->GetDispersion()*cluster2->GetDispersion();
433 l12 = cluster2->GetM20();
434 l02 = cluster2->GetM02();
435 tof2 = cluster2->GetTOF()*1e9;
437 // else printf("cluster2 not available: calo label %d / %d, cluster ID %d\n",
438 // photon2->GetCaloLabel(0),(GetReader()->GetInputEvent())->GetNumberOfCaloClusters()-1,cluster2->GetID());
440 //Select clusters with good time window difference
441 Double_t t12diff = tof1-tof2;
442 Float_t asymmetry = TMath::Abs(e1-e2)/(e1+e2);
443 fhClusterPairDiffTimeE ->Fill(e1+e2, t12diff);
444 fhClusterPairDiffTimeAsy->Fill(asymmetry,t12diff);
445 if(TMath::Abs(t12diff) > GetPairTimeCut()) continue;
448 //Select good pair (good phi, pt cuts, aperture and invariant mass)
449 if(GetNeutralMesonSelection()->SelectPair(mom1, mom2))
452 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());
454 //Play with the MC stack if available
456 //Check origin of the candidates
457 Int_t label1 = photon1->GetLabel();
458 Int_t label2 = photon2->GetLabel();
459 if(label1>=0)tag1 = GetMCAnalysisUtils()->CheckOrigin(label1, GetReader(), photon1->GetInputFileIndex());
460 if(label2>=0)tag2 = GetMCAnalysisUtils()->CheckOrigin(label2, GetReader(), photon2->GetInputFileIndex());
462 if(GetDebug() > 0) printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - Origin of: photon1 %d; photon2 %d \n",tag1, tag2);
463 if(GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCPi0Decay) &&
464 GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCPi0Decay)){
466 //Check if pi0 mother is the same
467 if(GetReader()->ReadStack()){
469 TParticle * mother1 = GetMCStack()->Particle(label1);//photon in kine tree
470 label1 = mother1->GetFirstMother();
471 //mother1 = GetMCStack()->Particle(label1);//pi0
474 TParticle * mother2 = GetMCStack()->Particle(label2);//photon in kine tree
475 label2 = mother2->GetFirstMother();
476 //mother2 = GetMCStack()->Particle(label2);//pi0
479 else if(GetReader()->ReadAODMCParticles()){//&& (input > -1)){
481 AliAODMCParticle * mother1 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles(photon1->GetInputFileIndex()))->At(label1);//photon in kine tree
482 label1 = mother1->GetMother();
483 //mother1 = GetMCStack()->Particle(label1);//pi0
486 AliAODMCParticle * mother2 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles(photon2->GetInputFileIndex()))->At(label2);//photon in kine tree
487 label2 = mother2->GetMother();
488 //mother2 = GetMCStack()->Particle(label2);//pi0
492 //printf("mother1 %d, mother2 %d\n",label1,label2);
493 if(label1 == label2 && label1>=0)
494 GetMCAnalysisUtils()->SetTagBit(tag,AliMCAnalysisUtils::kMCPi0);
496 }//Work with stack also
498 //Fill some histograms about shower shape
499 if(clusters && GetReader()->GetDataType()!=AliCaloTrackReader::kMC){
502 //printf("Signal Cl1: e %f, pt %f, disp %f, l1 %f, l0 %f, eta %f, phi %f \n",
503 // e,pt,disp,l1,l0,photon2->Eta(),photon2->Phi());
505 fhEDispersion->Fill(e1, disp1);
506 fhELambda0 ->Fill(e1, l01 );
507 fhELambda1 ->Fill(e1, l11 );
510 //printf("Signal Cl2: e %f, pt %f, disp %f, l1 %f, l0 %f, eta %f, phi %f \n",e
511 // ,pt,disp,l1,l0,photon2->Eta(),photon2->Phi());
513 fhEDispersion->Fill(e2, disp2);
514 fhELambda0 ->Fill(e2, l02 );
515 fhELambda1 ->Fill(e2, l12 );
519 if( GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCPhoton) &&
520 !GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCConversion) ){
521 fhEMCLambda0[0] ->Fill(e1, l01);
522 fhEMCLambda1[0] ->Fill(e1, l11);
523 fhEMCDispersion[0] ->Fill(e1, disp1);
524 }//photon no conversion
525 else if ( GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCElectron)){
526 fhEMCLambda0[3] ->Fill(e1, l01);
527 fhEMCLambda1[3] ->Fill(e1, l11);
528 fhEMCDispersion[3] ->Fill(e1, disp1);
530 else if ( GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCConversion) &&
531 GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCConversion) ){
532 fhEMCLambda0[1] ->Fill(e1, l01);
533 fhEMCLambda1[1] ->Fill(e1, l11);
534 fhEMCDispersion[1] ->Fill(e1, disp1);
536 else if ( GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCPi0) ){
537 fhEMCLambda0[2] ->Fill(e1, l01);
538 fhEMCLambda1[2] ->Fill(e1, l11);
539 fhEMCDispersion[2] ->Fill(e1, disp1);
542 fhEMCLambda0[4] ->Fill(e1, l01);
543 fhEMCLambda1[4] ->Fill(e1, l11);
544 fhEMCDispersion[4] ->Fill(e1, disp1);
548 if( GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCPhoton) &&
549 !GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCConversion) ){
550 fhEMCLambda0[0] ->Fill(e2, l02);
551 fhEMCLambda1[0] ->Fill(e2, l12);
552 fhEMCDispersion[0] ->Fill(e2, disp2);
553 }//photon no conversion
554 else if ( GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCElectron)){
555 fhEMCLambda0[3] ->Fill(e2, l02);
556 fhEMCLambda1[3] ->Fill(e2, l12);
557 fhEMCDispersion[3] ->Fill(e2, disp2);
559 else if ( GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCConversion) &&
560 GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCConversion) ){
561 fhEMCLambda0[1] ->Fill(e2, l02);
562 fhEMCLambda1[1] ->Fill(e2, l12);
563 fhEMCDispersion[1] ->Fill(e2, disp2);
565 else if ( GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCPi0) ){
566 fhEMCLambda0[2] ->Fill(e2, l02);
567 fhEMCLambda1[2] ->Fill(e2, l12);
568 fhEMCDispersion[2] ->Fill(e2, disp2);
571 fhEMCLambda0[4] ->Fill(e2, l02);
572 fhEMCLambda1[4] ->Fill(e2, l12);
573 fhEMCDispersion[4] ->Fill(e2, disp2);
578 //Create AOD for analysis
580 AliAODPWG4Particle pi0 = AliAODPWG4Particle(mom);
581 //pi0.SetLabel(calo->GetLabel());
582 pi0.SetIdentifiedParticleType(AliCaloPID::kPi0);
583 pi0.SetDetector(photon1->GetDetector());
585 //Set the indeces of the original caloclusters
586 pi0.SetCaloLabel(photon1->GetCaloLabel(0), photon2->GetCaloLabel(0));
587 //pi0.SetInputFileIndex(input);
595 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - End fill AODs \n");
599 //__________________________________________________________________
600 void AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS()
602 //Do analysis and fill aods
603 //Search for the photon decay in calorimeters
604 //Read photon list from AOD, produced in class AliAnaPhoton and AliGammaConversion
605 //Check if 2 photons have the mass of the pi0.
614 if(!GetInputAODBranch()){
615 printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - No input calo photons in AOD branch with name < %s > , STOP\n",GetInputAODName().Data());
619 for(Int_t iphoton = 0; iphoton < GetInputAODBranch()->GetEntriesFast(); iphoton++){
620 AliAODPWG4Particle * photon1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(iphoton));
621 mom1 = *(photon1->Momentum());
623 //Play with the MC stack if available
624 fInputAODGammaConv = (TClonesArray *) GetReader()->GetOutputEvent()->FindListObject(fInputAODGammaConvName);
625 if(!fInputAODGammaConv) {
626 printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - No input gamma conversions in AOD branch with name < %s >, STOP\n",fInputAODGammaConvName.Data());
629 for(Int_t jphoton = 0; jphoton < fInputAODGammaConv->GetEntriesFast(); jphoton++){
630 AliAODPWG4Particle * photon2 = (AliAODPWG4Particle*) (fInputAODGammaConv->At(jphoton));
632 evtIndex = GetMixedEvent()->EventIndexForCaloCluster(photon2->GetCaloLabel(0)) ;
633 if(TMath::Abs(GetVertex(evtIndex)[2]) > GetZvertexCut()) continue ; //vertex cut
635 mom2 = *(photon2->Momentum());
637 //Int_t input = -1; //if -1 photons come from different files, not a pi0
638 //if(photon1->GetInputFileIndex() == photon2->GetInputFileIndex()) input = photon1->GetInputFileIndex();
640 //Select good pair (good phi, pt cuts, aperture and invariant mass)
641 if(GetNeutralMesonSelection()->SelectPair(mom1, mom2)){
642 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());
645 Int_t label1 = photon1->GetLabel();
646 Int_t label2 = photon2->GetLabel();
647 if(label1>=0)tag1 = GetMCAnalysisUtils()->CheckOrigin(label1, GetReader(), photon1->GetInputFileIndex());
648 if(label2>=0)tag2 = GetMCAnalysisUtils()->CheckOrigin(label2, GetReader(), photon2->GetInputFileIndex());
649 if(GetDebug() > 0) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - Origin of: photon1 %d; photon2 %d \n",tag1, tag2);
650 if(GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCPi0Decay) &&
651 GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCPi0Decay)){
652 //Check if pi0 mother is the same
654 if(GetReader()->ReadStack()){
656 TParticle * mother1 = GetMCStack()->Particle(label1);//photon in kine tree
657 label1 = mother1->GetFirstMother();
658 //mother1 = GetMCStack()->Particle(label1);//pi0
661 TParticle * mother2 = GetMCStack()->Particle(label2);//photon in kine tree
662 label2 = mother2->GetFirstMother();
663 //mother2 = GetMCStack()->Particle(label2);//pi0
666 else if(GetReader()->ReadAODMCParticles()&& label1>=0 && label2>=0){ //&& (input > -1)){
668 AliAODMCParticle * mother1 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles(photon1->GetInputFileIndex()))->At(label1);//photon in kine tree
669 label1 = mother1->GetMother();
670 //mother1 = GetMCStack()->Particle(label1);//pi0
673 AliAODMCParticle * mother2 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles(photon2->GetInputFileIndex()))->At(label2);//photon in kine tree
674 label2 = mother2->GetMother();
675 //mother2 = GetMCStack()->Particle(label2);//pi0
679 //printf("mother1 %d, mother2 %d\n",label1,label2);
680 if(label1 == label2 && label1>=0)
681 GetMCAnalysisUtils()->SetTagBit(tag,AliMCAnalysisUtils::kMCPi0);
683 }//Work with stack also
685 //Create AOD for analysis
687 AliAODPWG4Particle pi0 = AliAODPWG4Particle(mom);
688 //pi0.SetLabel(calo->GetLabel());
689 pi0.SetIdentifiedParticleType(AliCaloPID::kPi0);
690 pi0.SetDetector(photon1->GetDetector());
692 //Set the indeces of the original tracks or caloclusters
693 pi0.SetCaloLabel(photon1->GetCaloLabel(0), -1);
694 pi0.SetTrackLabel(photon2->GetTrackLabel(0), photon2->GetTrackLabel(1));
695 //pi0.SetInputFileIndex(input);
702 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - End fill AODs \n");
707 //__________________________________________________________________
708 void AliAnaPi0EbE::MakeShowerShapeIdentification()
710 //Search for pi0 in fCalorimeter with shower shape analysis
712 TObjArray * pl = 0x0;
713 //Select the Calorimeter of the photon
714 if(fCalorimeter == "PHOS")
715 pl = GetPHOSClusters();
716 else if (fCalorimeter == "EMCAL")
717 pl = GetEMCALClusters();
720 Info("MakeShowerShapeIdentification","TObjArray with %s clusters is NULL!\n",fCalorimeter.Data());
725 for(Int_t icalo = 0; icalo < pl->GetEntriesFast(); icalo++){
726 AliVCluster * calo = (AliVCluster*) (pl->At(icalo));
729 if (GetMixedEvent()) {
730 evtIndex=GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ;
732 if(TMath::Abs(GetVertex(evtIndex)[2]) > GetZvertexCut()) continue ; //vertex cut
734 //Get Momentum vector,
735 if(GetReader()->GetDataType() != AliCaloTrackReader::kMC){
736 calo->GetMomentum(mom,GetVertex(evtIndex)) ;}//Assume that come from vertex in straight line
738 Double_t vertex[]={0,0,0};
739 calo->GetMomentum(mom,vertex) ;
742 //If too small or big pt, skip it
743 if(mom.Pt() < GetMinPt() || mom.Pt() > GetMaxPt() ) continue ;
744 //Check acceptance selection
745 if(IsFiducialCutOn()){
746 Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter) ;
750 //Create AOD for analysis
751 AliAODPWG4Particle aodpi0 = AliAODPWG4Particle(mom);
752 aodpi0.SetLabel(calo->GetLabel());
753 //Set the indeces of the original caloclusters
754 aodpi0.SetCaloLabel(calo->GetID(),-1);
755 aodpi0.SetDetector(fCalorimeter);
757 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());
759 //Check Distance to Bad channel, set bit.
760 Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
761 if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;
762 if(distBad < fMinDist) //In bad channel (PHOS cristal size 2.2x2.2 cm)
765 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - FillAOD: Bad channel cut passed %4.2f\n",distBad);
767 if(distBad > fMinDist3) aodpi0.SetDistToBad(2) ;
768 else if(distBad > fMinDist2) aodpi0.SetDistToBad(1) ;
769 else aodpi0.SetDistToBad(0) ;
772 //PID selection or bit setting
773 if(GetReader()->GetDataType() == AliCaloTrackReader::kMC){
774 //Get most probable PID, check PID weights (in MC this option is mandatory)
775 aodpi0.SetIdentifiedParticleType(GetCaloPID()->GetIdentifiedParticleType(fCalorimeter,calo->GetPID(),mom.E()));//PID with weights
777 printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - FillAOD: PDG of identified particle %d\n",aodpi0.GetIdentifiedParticleType());
778 //If primary is not pi0, skip it.
779 if(aodpi0.GetIdentifiedParticleType() != AliCaloPID::kPi0) continue ;
781 else if(IsCaloPIDOn()){
782 //Skip matched clusters with tracks
783 if(IsTrackMatched(calo)) continue ;
785 //Get most probable PID, 2 options check PID weights
786 //or redo PID, recommended option for EMCal.
787 if(!IsCaloPIDRecalculationOn())
788 aodpi0.SetIdentifiedParticleType(GetCaloPID()->GetIdentifiedParticleType(fCalorimeter,calo->GetPID(),mom.E()));//PID with weights
790 aodpi0.SetIdentifiedParticleType(GetCaloPID()->GetIdentifiedParticleType(fCalorimeter,mom,calo));//PID recalculated
792 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - PDG of identified particle %d\n",aodpi0.GetIdentifiedParticleType());
794 //If cluster does not pass pid, not pi0, skip it.
795 if(aodpi0.GetIdentifiedParticleType() != AliCaloPID::kPi0) continue ;
799 //Set PID bits for later selection (AliAnaPi0 for example)
800 //GetPDG already called in SetPIDBits.
801 GetCaloPID()->SetPIDBits(fCalorimeter,calo,&aodpi0, GetCaloUtils());
802 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - PID Bits set \n");
805 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - Pi0 selection cuts passed: pT %3.2f, pdg %d\n",aodpi0.Pt(), aodpi0.GetIdentifiedParticleType());
807 //Play with the MC stack if available
808 //Check origin of the candidates
810 if((GetReader()->GetDataType() == AliCaloTrackReader::kMC && fAnaType!=kSSCalo) ||
811 GetReader()->GetDataType() != AliCaloTrackReader::kMC){
812 //aodpi0.SetInputFileIndex(input);
814 tag = GetMCAnalysisUtils()->CheckOrigin(calo->GetLabel(),GetReader(), aodpi0.GetInputFileIndex());
815 //GetMCAnalysisUtils()->CheckMultipleOrigin(calo->GetLabels(),calo->GetNLabels(), GetReader(), aodpi0.GetInputFileIndex(), tag);
817 if(GetDebug() > 0) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - Origin of candidate %d\n",aodpi0.GetTag());
819 }//Work with stack also
821 //Add AOD with pi0 object to aod branch
822 AddAODParticle(aodpi0);
826 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - End fill AODs \n");
829 //__________________________________________________________________
830 void AliAnaPi0EbE::MakeAnalysisFillHistograms()
832 //Do analysis and fill histograms
834 if(!GetOutputAODBranch()){
835 printf("AliAnaPi0EbE::MakeAnalysisFillHistograms() - No output pi0 in AOD branch with name < %s >,STOP \n",GetOutputAODName().Data());
838 //Loop on stored AOD pi0
839 Int_t naod = GetOutputAODBranch()->GetEntriesFast();
840 if(GetDebug() > 0) printf("AliAnaPi0EbE::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod);
842 for(Int_t iaod = 0; iaod < naod ; iaod++){
844 AliAODPWG4Particle* pi0 = (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
845 Int_t pdg = pi0->GetIdentifiedParticleType();
847 if(IsCaloPIDOn() && pdg != AliCaloPID::kPi0) continue;
849 //Fill pi0 histograms
850 Float_t ener = pi0->E();
851 Float_t pt = pi0->Pt();
852 Float_t phi = pi0->Phi();
853 if(phi < 0) phi+=TMath::TwoPi();
854 Float_t eta = pi0->Eta();
858 fhEEtaPhiPi0 ->Fill(ener,eta,phi);
861 if((GetReader()->GetDataType() == AliCaloTrackReader::kMC && fAnaType!=kSSCalo) ||
862 GetReader()->GetDataType() != AliCaloTrackReader::kMC){
863 if(GetMCAnalysisUtils()->CheckTagBit(pi0->GetTag(), AliMCAnalysisUtils::kMCPi0)){
864 fhPtMCPi0 ->Fill(pt);
865 fhPhiMCPi0 ->Fill(pt,phi);
866 fhEtaMCPi0 ->Fill(pt,eta);
869 fhPtMCNoPi0 ->Fill(pt);
870 fhPhiMCNoPi0 ->Fill(pt,phi);
871 fhEtaMCNoPi0 ->Fill(pt,eta);
874 }//Histograms with MC
880 //__________________________________________________________________
881 void AliAnaPi0EbE::Print(const Option_t * opt) const
883 //Print some relevant parameters set for the analysis
887 printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
888 AliAnaPartCorrBaseClass::Print("");
889 printf("Analysis Type = %d \n", fAnaType) ;
890 if(fAnaType == kSSCalo){
891 printf("Calorimeter = %s\n", fCalorimeter.Data()) ;
892 printf("Min Distance to Bad Channel = %2.1f\n",fMinDist);
893 printf("Min Distance to Bad Channel 2 = %2.1f\n",fMinDist2);
894 printf("Min Distance to Bad Channel 3 = %2.1f\n",fMinDist3);