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 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  = GetHistoPtBins();
 123   Int_t nphibins = GetHistoPhiBins();
 124   Int_t netabins = GetHistoEtaBins();
 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 = 0;
 268   Int_t tag2 = 0;
 269   Int_t tag  = 0;
 270
 271   if(!GetInputAODBranch()){
 272     printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - 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       AliAODPWG4Particle * photon2 =  (AliAODPWG4Particle*) (GetInputAODBranch()->At(jphoton));
 283       mom2 = *(photon2->Momentum());
 284           Int_t input = -1;     //if -1 photons come from different files, not a pi0
 285           if(photon1->GetInputFileIndex() == photon2->GetInputFileIndex()) input = photon1->GetInputFileIndex();
 286
 287       //Select good pair (good phi, pt cuts, aperture and invariant mass)
 288       if(GetNeutralMesonSelection()->SelectPair(mom1, mom2))
 289         {
 290           if(GetDebug()>1)
 291             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());
 292
 293           //Play with the MC stack if available
 294           if(IsDataMC()){
 295             //Check origin of the candidates
 296                 Int_t  label1 = photon1->GetLabel();
 297                 Int_t  label2 = photon2->GetLabel();
 298             tag1 = GetMCAnalysisUtils()->CheckOrigin(label1, GetReader(), photon1->GetInputFileIndex());
 299             tag2 = GetMCAnalysisUtils()->CheckOrigin(label2, GetReader(), photon2->GetInputFileIndex());
 300
 301             if(GetDebug() > 0) printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - Origin of: photon1 %d; photon2 %d \n",tag1, tag2);
 302             if(GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCPi0Decay) && GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCPi0Decay)){
 303
 304               //Check if pi0 mother is the same
 305                   if(GetReader()->ReadStack()){
 306                                 TParticle * mother1 = GetMCStack()->Particle(label1);//photon in kine tree
 307                                 label1 = mother1->GetFirstMother();
 308                                 //mother1 = GetMCStack()->Particle(label1);//pi0
 309
 310                                 TParticle * mother2 = GetMCStack()->Particle(label2);//photon in kine tree
 311                                 label2 = mother2->GetFirstMother();
 312                                 //mother2 = GetMCStack()->Particle(label2);//pi0
 313                   }
 314                   else if(GetReader()->ReadAODMCParticles() && (input > -1)){
 315                                 AliAODMCParticle * mother1 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles(photon1->GetInputFileIndex()))->At(label1);//photon in kine tree
 316                                 label1 = mother1->GetMother();
 317                                 //mother1 = GetMCStack()->Particle(label1);//pi0
 318                                 AliAODMCParticle * mother2 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles(photon2->GetInputFileIndex()))->At(label2);//photon in kine tree
 319                                 label2 = mother2->GetMother();
 320                                 //mother2 = GetMCStack()->Particle(label2);//pi0
 321                   }
 322
 323               //printf("mother1 %d, mother2 %d\n",label1,label2);
 324               if(label1 == label2)
 325                           GetMCAnalysisUtils()->SetTagBit(tag,AliMCAnalysisUtils::kMCPi0);
 326             }
 327           }//Work with stack also
 328
 329           //Create AOD for analysis
 330           mom = mom1+mom2;
 331           AliAODPWG4Particle pi0 = AliAODPWG4Particle(mom);
 332           //pi0.SetLabel(calo->GetLabel(0));
 333           pi0.SetPdg(AliCaloPID::kPi0);
 334           pi0.SetDetector(photon1->GetDetector());
 335           pi0.SetTag(tag);
 336           //Set the indeces of the original caloclusters
 337           pi0.SetCaloLabel(photon1->GetCaloLabel(0), photon2->GetCaloLabel(0));
 338           pi0.SetInputFileIndex(input);
 339           AddAODParticle(pi0);
 340         }//pi0
 341     }//2n photon loop
 342
 343   }//1st photon loop
 344
 345   if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - End fill AODs \n");
 346
 347 }
 348
 349 //__________________________________________________________________
 350 void  AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS()
 351 {
 352   //Do analysis and fill aods
 353   //Search for the photon decay in calorimeters
 354   //Read photon list from AOD, produced in class AliAnaPhoton and AliGammaConversion
 355   //Check if 2 photons have the mass of the pi0.
 356
 357   TLorentzVector mom1;
 358   TLorentzVector mom2;
 359   TLorentzVector mom ;
 360   Int_t tag1 = 0;
 361   Int_t tag2 = 0;
 362   Int_t tag  = 0;
 363
 364   if(!GetInputAODBranch()){
 365     printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - No input calo photons in AOD branch with name < %s > , STOP\n",GetInputAODName().Data());
 366     abort();
 367   }
 368   for(Int_t iphoton = 0; iphoton < GetInputAODBranch()->GetEntriesFast(); iphoton++){
 369     AliAODPWG4Particle * photon1 =  (AliAODPWG4Particle*) (GetInputAODBranch()->At(iphoton));
 370     mom1 = *(photon1->Momentum());
 371
 372     //Play with the MC stack if available
 373     fInputAODGammaConv = (TClonesArray *) GetReader()->GetOutputEvent()->FindListObject(fInputAODGammaConvName);
 374     if(!fInputAODGammaConv) {
 375       printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - No input gamma conversions in AOD branch with name < %s >, STOP\n",fInputAODGammaConvName.Data());
 376       abort();
 377     }
 378     for(Int_t jphoton = iphoton+1; jphoton < fInputAODGammaConv->GetEntriesFast()-1; jphoton++){
 379       AliAODPWG4Particle * photon2 =  (AliAODPWG4Particle*) (fInputAODGammaConv->At(jphoton));
 380       mom2 = *(photon2->Momentum());
 381
 382           Int_t input = -1;     //if -1 photons come from different files, not a pi0
 383           if(photon1->GetInputFileIndex() == photon2->GetInputFileIndex()) input = photon1->GetInputFileIndex();
 384
 385       //Select good pair (good phi, pt cuts, aperture and invariant mass)
 386       if(GetNeutralMesonSelection()->SelectPair(mom1, mom2)){
 387         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());
 388
 389         if(IsDataMC()){
 390           Int_t label1 = photon1->GetLabel();
 391           Int_t label2 = photon2->GetLabel();
 392           tag1 = GetMCAnalysisUtils()->CheckOrigin(label1, GetReader(), photon1->GetInputFileIndex());
 393           tag2 = GetMCAnalysisUtils()->CheckOrigin(label2, GetReader(), photon2->GetInputFileIndex());
 394           if(GetDebug() > 0) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - Origin of: photon1 %d; photon2 %d \n",tag1, tag2);
 395           if(GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCPi0Decay) && GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCPi0Decay)){
 396             //Check if pi0 mother is the same
 397
 398                 if(GetReader()->ReadStack()){
 399                         TParticle * mother1 = GetMCStack()->Particle(label1);//photon in kine tree
 400                         label1 = mother1->GetFirstMother();
 401                         //mother1 = GetMCStack()->Particle(label1);//pi0
 402
 403                         TParticle * mother2 = GetMCStack()->Particle(label2);//photon in kine tree
 404                         label2 = mother2->GetFirstMother();
 405                         //mother2 = GetMCStack()->Particle(label2);//pi0
 406             }
 407                 else if(GetReader()->ReadAODMCParticles() && (input > -1)){
 408                         AliAODMCParticle * mother1 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles(photon1->GetInputFileIndex()))->At(label1);//photon in kine tree
 409                         label1 = mother1->GetMother();
 410                         //mother1 = GetMCStack()->Particle(label1);//pi0
 411                         AliAODMCParticle * mother2 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles(photon2->GetInputFileIndex()))->At(label2);//photon in kine tree
 412                         label2 = mother2->GetMother();
 413                         //mother2 = GetMCStack()->Particle(label2);//pi0
 414                 }
 415
 416                 //printf("mother1 %d, mother2 %d\n",label1,label2);
 417             if(label1 == label2)
 418               GetMCAnalysisUtils()->SetTagBit(tag,AliMCAnalysisUtils::kMCPi0);
 419           }
 420         }//Work with stack also
 421
 422         //Create AOD for analysis
 423         mom = mom1+mom2;
 424         AliAODPWG4Particle pi0 = AliAODPWG4Particle(mom);
 425         //pi0.SetLabel(calo->GetLabel(0));
 426         pi0.SetPdg(AliCaloPID::kPi0);
 427         pi0.SetDetector(photon1->GetDetector());
 428         pi0.SetTag(tag);
 429         //Set the indeces of the original tracks or caloclusters
 430         pi0.SetCaloLabel(photon1->GetCaloLabel(0), -1);
 431         pi0.SetTrackLabel(photon2->GetTrackLabel(0), photon2->GetTrackLabel(1));
 432         pi0.SetInputFileIndex(input);
 433         AddAODParticle(pi0);
 434       }//pi0
 435     }//2n photon loop
 436
 437   }//1st photon loop
 438
 439   if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - End fill AODs \n");
 440
 441 }
 442
 443
 444 //__________________________________________________________________
 445 void  AliAnaPi0EbE::MakeShowerShapeIdentification()
 446 {
 447   //Search for pi0 in fCalorimeter with shower shape analysis
 448
 449   TObjArray * pl = new TObjArray;
 450
 451   //Get vertex for photon momentum calculation
 452   Double_t vertex[]  = {0,0,0} ; //vertex
 453   Double_t vertex2[] = {0,0,0} ; //vertex from second aod input
 454   if(GetReader()->GetDataType() != AliCaloTrackReader::kMC)
 455   {
 456           GetReader()->GetVertex(vertex);
 457           if(GetReader()->GetSecondInputAODTree()) GetReader()->GetSecondInputAODVertex(vertex2);
 458   }
 459
 460   //Select the Calorimeter of the photon
 461   if(fCalorimeter == "PHOS")
 462     pl = GetAODPHOS();
 463   else if (fCalorimeter == "EMCAL")
 464     pl = GetAODEMCAL();
 465   //Fill AODCaloClusters and AODParticle with PHOS aods
 466   TLorentzVector mom ;
 467   for(Int_t icalo = 0; icalo < pl->GetEntriesFast(); icalo++){
 468     AliAODCaloCluster * calo = (AliAODCaloCluster*) (pl->At(icalo));
 469
 470     //Cluster selection, not charged, with pi0 id and in fiducial cut
 471
 472         //Input from second AOD?
 473         Int_t input = 0;
 474         if     (fCalorimeter == "EMCAL" && GetReader()->GetAODEMCALNormalInputEntries() <= icalo) input = 1 ;
 475         else if(fCalorimeter == "PHOS"  && GetReader()->GetAODPHOSNormalInputEntries()  <= icalo) input = 1;
 476
 477         //Get Momentum vector,
 478         if     (input == 0) calo->GetMomentum(mom,vertex) ;//Assume that come from vertex in straight line
 479         else if(input == 1) calo->GetMomentum(mom,vertex2);//Assume that come from vertex in straight line
 480
 481         //If too small or big pt, skip it
 482     if(mom.Pt() < GetMinPt() || mom.Pt() > GetMaxPt() ) continue ;
 483     //Check acceptance selection
 484     if(IsFiducialCutOn()){
 485       Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter) ;
 486       if(! in ) continue ;
 487     }
 488
 489     //Create AOD for analysis
 490     AliAODPWG4Particle aodpi0 = AliAODPWG4Particle(mom);
 491     aodpi0.SetLabel(calo->GetLabel(0));
 492     //Set the indeces of the original caloclusters
 493     aodpi0.SetCaloLabel(calo->GetID(),-1);
 494     aodpi0.SetDetector(fCalorimeter);
 495     if(GetDebug() > 1)
 496       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());
 497
 498     //Check Distance to Bad channel, set bit.
 499     Double_t distBad=calo->GetDistToBadChannel() ; //Distance to bad channel
 500     if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;
 501     if(distBad < fMinDist) //In bad channel (PHOS cristal size 2.2x2.2 cm)
 502       continue ;
 503
 504     if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - FillAOD: Bad channel cut passed %4.2f\n",distBad);
 505
 506     if(distBad > fMinDist3) aodpi0.SetDistToBad(2) ;
 507     else if(distBad > fMinDist2) aodpi0.SetDistToBad(1) ;
 508     else aodpi0.SetDistToBad(0) ;
 509
 510     //Check PID
 511     //PID selection or bit setting
 512     if(GetReader()->GetDataType() == AliCaloTrackReader::kMC){
 513       //Get most probable PID, check PID weights (in MC this option is mandatory)
 514       aodpi0.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,calo->PID(),mom.E()));//PID with weights
 515       if(GetDebug() > 1)
 516         printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - FillAOD: PDG of identified particle %d\n",aodpi0.GetPdg());
 517       //If primary is not pi0, skip it.
 518       if(aodpi0.GetPdg() != AliCaloPID::kPi0) continue ;
 519     }
 520     else if(IsCaloPIDOn()){
 521       //Skip matched clusters with tracks
 522       if(calo->GetNTracksMatched() > 0) continue ;
 523
 524       //Get most probable PID, 2 options check PID weights
 525       //or redo PID, recommended option for EMCal.
 526       if(!IsCaloPIDRecalculationOn())
 527         aodpi0.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,calo->PID(),mom.E()));//PID with weights
 528       else
 529         aodpi0.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,mom,calo));//PID recalculated
 530
 531       if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - PDG of identified particle %d\n",aodpi0.GetPdg());
 532
 533       //If cluster does not pass pid, not pi0, skip it.
 534       if(aodpi0.GetPdg() != AliCaloPID::kPi0) continue ;
 535
 536     }
 537     else{
 538       //Set PID bits for later selection (AliAnaPi0 for example)
 539       //GetPDG already called in SetPIDBits.
 540       GetCaloPID()->SetPIDBits(fCalorimeter,calo,&aodpi0);
 541       if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - PID Bits set \n");
 542     }
 543
 544     if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - Pi0 selection cuts passed: pT %3.2f, pdg %d\n",aodpi0.Pt(), aodpi0.GetPdg());
 545
 546     //Play with the MC stack if available
 547     //Check origin of the candidates
 548     if(IsDataMC()){
 549       if((GetReader()->GetDataType() == AliCaloTrackReader::kMC && fAnaType!=kSSCalo) ||
 550                  GetReader()->GetDataType() != AliCaloTrackReader::kMC){
 551                   aodpi0.SetInputFileIndex(input);
 552                   Int_t tag     =0;
 553                   tag = GetMCAnalysisUtils()->CheckOrigin(calo->GetLabel(0),GetReader(), aodpi0.GetInputFileIndex());
 554                   //GetMCAnalysisUtils()->CheckMultipleOrigin(calo->GetLabels(),calo->GetNLabel(), GetReader(), aodpi0.GetInputFileIndex(), tag);
 555                   aodpi0.SetTag(tag);
 556                   if(GetDebug() > 0) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - Origin of candidate %d\n",aodpi0.GetTag());
 557       }
 558     }//Work with stack also
 559
 560     //Add AOD with pi0 object to aod branch
 561     AddAODParticle(aodpi0);
 562
 563   }//loop
 564
 565   if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - End fill AODs \n");
 566
 567 }
 568 //__________________________________________________________________
 569 void  AliAnaPi0EbE::MakeAnalysisFillHistograms()
 570 {
 571   //Do analysis and fill histograms
 572
 573   if(!GetOutputAODBranch()){
 574     printf("AliAnaPi0EbE::MakeAnalysisFillHistograms()  - No output pi0 in AOD branch with name < %s >,STOP \n",GetOutputAODName().Data());
 575     abort();
 576   }
 577   //Loop on stored AOD pi0
 578   Int_t naod = GetOutputAODBranch()->GetEntriesFast();
 579   if(GetDebug() > 0) printf("AliAnaPi0EbE::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod);
 580
 581   for(Int_t iaod = 0; iaod < naod ; iaod++){
 582
 583     AliAODPWG4Particle* pi0 =  (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
 584     Int_t pdg = pi0->GetPdg();
 585
 586     if(IsCaloPIDOn() && pdg != AliCaloPID::kPi0) continue;
 587
 588     //Fill pi0 histograms
 589     Float_t pt  = pi0->Pt();
 590     Float_t phi = pi0->Phi();
 591     Float_t eta = pi0->Eta();
 592
 593     fhPtPi0  ->Fill(pt);
 594     fhPhiPi0 ->Fill(pt,phi);
 595     fhEtaPi0 ->Fill(pt,eta);
 596
 597     if(IsDataMC()){
 598       if((GetReader()->GetDataType() == AliCaloTrackReader::kMC && fAnaType!=kSSCalo) ||
 599          GetReader()->GetDataType() != AliCaloTrackReader::kMC){
 600         if(GetMCAnalysisUtils()->CheckTagBit(pi0->GetTag(), AliMCAnalysisUtils::kMCPi0)){
 601           fhPtMCPi0  ->Fill(pt);
 602           fhPhiMCPi0 ->Fill(pt,phi);
 603           fhEtaMCPi0 ->Fill(pt,eta);
 604         }
 605         else{
 606           fhPtMCNoPi0  ->Fill(pt);
 607           fhPhiMCNoPi0 ->Fill(pt,phi);
 608           fhEtaMCNoPi0 ->Fill(pt,eta);
 609         }
 610       }
 611     }//Histograms with MC
 612
 613   }// aod loop
 614
 615 }
 616
 617
 618 //____________________________________________________________________________
 619 void AliAnaPi0EbE::Init()
 620 {
 621   //Init
 622   //Do some checks
 623   if(fCalorimeter == "PHOS" && !GetReader()->IsPHOSSwitchedOn()){
 624     printf("AliAnaPi0EbE::Init() - !!STOP: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!\n");
 625     abort();
 626   }
 627   else  if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn()){
 628     printf("AliAnaPi0EbE::Init() - !!STOP: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!\n");
 629     abort();
 630   }
 631
 632 }
 633
 634 //____________________________________________________________________________
 635 void AliAnaPi0EbE::InitParameters()
 636 {
 637   //Initialize the parameters of the analysis.
 638   SetOutputAODClassName("AliAODPWG4Particle");
 639   SetOutputAODName("Pi0");
 640
 641   AddToHistogramsName("AnaPi0EbE_");
 642
 643   fInputAODGammaConvName = "gammaconv" ;
 644   fAnaType = kIMCalo ;
 645   fCalorimeter = "EMCAL" ;
 646   fMinDist  = 2.;
 647   fMinDist2 = 4.;
 648   fMinDist3 = 5.;
 649
 650 }
 651
 652 //__________________________________________________________________
 653 void AliAnaPi0EbE::Print(const Option_t * opt) const
 654 {
 655   //Print some relevant parameters set for the analysis
 656   if(! opt)
 657     return;
 658
 659   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
 660   AliAnaPartCorrBaseClass::Print("");
 661   printf("Analysis Type = %d \n",  fAnaType) ;
 662   if(fAnaType == kSSCalo){
 663     printf("Calorimeter            =     %s\n", fCalorimeter.Data()) ;
 664     printf("Min Distance to Bad Channel   = %2.1f\n",fMinDist);
 665     printf("Min Distance to Bad Channel 2 = %2.1f\n",fMinDist2);
 666     printf("Min Distance to Bad Channel 3 = %2.1f\n",fMinDist3);
 667   }
 668   printf("    \n") ;
 669
 670 }