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