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