PWG4/PartCorrDep/AliAnaPi0EbE.cxx

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