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     if(TMath::Abs(GetVertex(evtIndex1)[2]) > GetZvertexCut()) continue ;  //vertex cut
 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       if(TMath::Abs(GetVertex(evtIndex2)[2]) > GetZvertexCut()) continue ;  //vertex cut
 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   Int_t evtIndex = 0;
 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       if(GetMixedEvent())
 429         evtIndex = GetMixedEvent()->EventIndexForCaloCluster(photon2->GetCaloLabel(0)) ;
 430       if(TMath::Abs(GetVertex(evtIndex)[2]) > GetZvertexCut()) continue ;  //vertex cut
 431
 432       mom2 = *(photon2->Momentum());
 433
 434       //Int_t input = -1;       //if -1 photons come from different files, not a pi0
 435       //if(photon1->GetInputFileIndex() == photon2->GetInputFileIndex()) input = photon1->GetInputFileIndex();
 436
 437       //Select good pair (good phi, pt cuts, aperture and invariant mass)
 438       if(GetNeutralMesonSelection()->SelectPair(mom1, mom2)){
 439         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());
 440
 441         if(IsDataMC()){
 442           Int_t label1 = photon1->GetLabel();
 443           Int_t label2 = photon2->GetLabel();
 444           if(label1>=0)tag1 = GetMCAnalysisUtils()->CheckOrigin(label1, GetReader(), photon1->GetInputFileIndex());
 445           if(label2>=0)tag2 = GetMCAnalysisUtils()->CheckOrigin(label2, GetReader(), photon2->GetInputFileIndex());
 446           if(GetDebug() > 0) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - Origin of: photon1 %d; photon2 %d \n",tag1, tag2);
 447           if(GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCPi0Decay) &&
 448              GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCPi0Decay)){
 449             //Check if pi0 mother is the same
 450
 451             if(GetReader()->ReadStack()){
 452               if(label1>=0){
 453                 TParticle * mother1 = GetMCStack()->Particle(label1);//photon in kine tree
 454                 label1 = mother1->GetFirstMother();
 455                 //mother1 = GetMCStack()->Particle(label1);//pi0
 456               }
 457               if(label2>=0){
 458                 TParticle * mother2 = GetMCStack()->Particle(label2);//photon in kine tree
 459                 label2 = mother2->GetFirstMother();
 460                 //mother2 = GetMCStack()->Particle(label2);//pi0
 461               }
 462             }
 463             else if(GetReader()->ReadAODMCParticles()&& label1>=0 && label2>=0){ //&& (input > -1)){
 464               if(label1>=0){
 465                 AliAODMCParticle * mother1 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles(photon1->GetInputFileIndex()))->At(label1);//photon in kine tree
 466                 label1 = mother1->GetMother();
 467                 //mother1 = GetMCStack()->Particle(label1);//pi0
 468               }
 469               if(label2>=0){
 470                 AliAODMCParticle * mother2 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles(photon2->GetInputFileIndex()))->At(label2);//photon in kine tree
 471                 label2 = mother2->GetMother();
 472                 //mother2 = GetMCStack()->Particle(label2);//pi0
 473               }
 474             }
 475
 476             //printf("mother1 %d, mother2 %d\n",label1,label2);
 477             if(label1 == label2 && label1>=0)
 478               GetMCAnalysisUtils()->SetTagBit(tag,AliMCAnalysisUtils::kMCPi0);
 479           }
 480         }//Work with stack also
 481
 482         //Create AOD for analysis
 483         mom = mom1+mom2;
 484         AliAODPWG4Particle pi0 = AliAODPWG4Particle(mom);
 485         //pi0.SetLabel(calo->GetLabel());
 486         pi0.SetPdg(AliCaloPID::kPi0);
 487         pi0.SetDetector(photon1->GetDetector());
 488         pi0.SetTag(tag);
 489         //Set the indeces of the original tracks or caloclusters
 490         pi0.SetCaloLabel(photon1->GetCaloLabel(0), -1);
 491         pi0.SetTrackLabel(photon2->GetTrackLabel(0), photon2->GetTrackLabel(1));
 492         //pi0.SetInputFileIndex(input);
 493         AddAODParticle(pi0);
 494       }//pi0
 495     }//2n photon loop
 496
 497   }//1st photon loop
 498
 499   if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - End fill AODs \n");
 500
 501 }
 502
 503
 504 //__________________________________________________________________
 505 void  AliAnaPi0EbE::MakeShowerShapeIdentification()
 506 {
 507   //Search for pi0 in fCalorimeter with shower shape analysis
 508
 509   TObjArray * pl = 0x0;
 510   //Select the Calorimeter of the photon
 511   if(fCalorimeter == "PHOS")
 512     pl = GetAODPHOS();
 513   else if (fCalorimeter == "EMCAL")
 514     pl = GetAODEMCAL();
 515
 516   if(!pl) {
 517     Info("MakeShowerShapeIdentification","TObjArray with %s clusters is NULL!\n",fCalorimeter.Data());
 518     return;
 519   }
 520
 521   //Get vertex for photon momentum calculation
 522   //Double_t vertex2[] = {0,0,0} ; //vertex from second aod input
 523   //if(GetReader()->GetDataType() != AliCaloTrackReader::kMC)
 524   //{
 525   //if(GetReader()->GetSecondInputAODTree()) GetReader()->GetSecondInputAODVertex(vertex2);
 526   //}
 527
 528
 529   TLorentzVector mom ;
 530   for(Int_t icalo = 0; icalo < pl->GetEntriesFast(); icalo++){
 531     AliVCluster * calo = (AliVCluster*) (pl->At(icalo));
 532
 533     Int_t evtIndex = 0 ;
 534     if (GetMixedEvent()) {
 535       evtIndex=GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ;
 536     }
 537     if(TMath::Abs(GetVertex(evtIndex)[2]) > GetZvertexCut()) continue ;  //vertex cut
 538
 539     //Cluster selection, not charged, with pi0 id and in fiducial cut
 540
 541     //Input from second AOD?
 542     //Int_t input = 0;
 543     //  if     (fCalorimeter == "EMCAL" && GetReader()->GetAODEMCALNormalInputEntries() <= icalo) input = 1 ;
 544     //  else if(fCalorimeter == "PHOS"  && GetReader()->GetAODPHOSNormalInputEntries()  <= icalo) input = 1;
 545
 546     //Get Momentum vector,
 547     //if     (input == 0)
 548     if(GetReader()->GetDataType() != AliCaloTrackReader::kMC){
 549       calo->GetMomentum(mom,GetVertex(evtIndex)) ;}//Assume that come from vertex in straight line
 550     else{
 551       Double_t vertex[]={0,0,0};
 552       calo->GetMomentum(mom,vertex) ;
 553     }
 554     //else if(input == 1) calo->GetMomentum(mom,vertex2);//Assume that come from vertex in straight line
 555
 556     //If too small or big pt, skip it
 557     if(mom.Pt() < GetMinPt() || mom.Pt() > GetMaxPt() ) continue ;
 558     //Check acceptance selection
 559     if(IsFiducialCutOn()){
 560       Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter) ;
 561       if(! in ) continue ;
 562     }
 563
 564     //Create AOD for analysis
 565     AliAODPWG4Particle aodpi0 = AliAODPWG4Particle(mom);
 566     aodpi0.SetLabel(calo->GetLabel());
 567     //Set the indeces of the original caloclusters
 568     aodpi0.SetCaloLabel(calo->GetID(),-1);
 569     aodpi0.SetDetector(fCalorimeter);
 570     if(GetDebug() > 1)
 571       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());
 572
 573     //Check Distance to Bad channel, set bit.
 574     Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
 575     if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;
 576     if(distBad < fMinDist) //In bad channel (PHOS cristal size 2.2x2.2 cm)
 577       continue ;
 578
 579     if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - FillAOD: Bad channel cut passed %4.2f\n",distBad);
 580
 581     if(distBad > fMinDist3) aodpi0.SetDistToBad(2) ;
 582     else if(distBad > fMinDist2) aodpi0.SetDistToBad(1) ;
 583     else aodpi0.SetDistToBad(0) ;
 584
 585     //Check PID
 586     //PID selection or bit setting
 587     if(GetReader()->GetDataType() == AliCaloTrackReader::kMC){
 588       //Get most probable PID, check PID weights (in MC this option is mandatory)
 589       aodpi0.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,calo->GetPID(),mom.E()));//PID with weights
 590       if(GetDebug() > 1)
 591         printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - FillAOD: PDG of identified particle %d\n",aodpi0.GetPdg());
 592       //If primary is not pi0, skip it.
 593       if(aodpi0.GetPdg() != AliCaloPID::kPi0) continue ;
 594     }
 595     else if(IsCaloPIDOn()){
 596       //Skip matched clusters with tracks
 597       if(IsTrackMatched(calo)) continue ;
 598
 599       //Get most probable PID, 2 options check PID weights
 600       //or redo PID, recommended option for EMCal.
 601       if(!IsCaloPIDRecalculationOn())
 602         aodpi0.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,calo->GetPID(),mom.E()));//PID with weights
 603       else
 604         aodpi0.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,mom,calo));//PID recalculated
 605
 606       if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - PDG of identified particle %d\n",aodpi0.GetPdg());
 607
 608       //If cluster does not pass pid, not pi0, skip it.
 609       if(aodpi0.GetPdg() != AliCaloPID::kPi0) continue ;
 610
 611     }
 612     else{
 613       //Set PID bits for later selection (AliAnaPi0 for example)
 614       //GetPDG already called in SetPIDBits.
 615       GetCaloPID()->SetPIDBits(fCalorimeter,calo,&aodpi0, GetCaloUtils());
 616       if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - PID Bits set \n");
 617     }
 618
 619     if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - Pi0 selection cuts passed: pT %3.2f, pdg %d\n",aodpi0.Pt(), aodpi0.GetPdg());
 620
 621     //Play with the MC stack if available
 622     //Check origin of the candidates
 623     if(IsDataMC()){
 624       if((GetReader()->GetDataType() == AliCaloTrackReader::kMC && fAnaType!=kSSCalo) ||
 625          GetReader()->GetDataType() != AliCaloTrackReader::kMC){
 626         //aodpi0.SetInputFileIndex(input);
 627         Int_t tag       =0;
 628         tag = GetMCAnalysisUtils()->CheckOrigin(calo->GetLabel(),GetReader(), aodpi0.GetInputFileIndex());
 629         //GetMCAnalysisUtils()->CheckMultipleOrigin(calo->GetLabels(),calo->GetNLabels(), GetReader(), aodpi0.GetInputFileIndex(), tag);
 630         aodpi0.SetTag(tag);
 631         if(GetDebug() > 0) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - Origin of candidate %d\n",aodpi0.GetTag());
 632       }
 633     }//Work with stack also
 634
 635     //Add AOD with pi0 object to aod branch
 636     AddAODParticle(aodpi0);
 637
 638   }//loop
 639
 640   if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - End fill AODs \n");
 641
 642 }
 643 //__________________________________________________________________
 644 void  AliAnaPi0EbE::MakeAnalysisFillHistograms()
 645 {
 646   //Do analysis and fill histograms
 647
 648   if(!GetOutputAODBranch()){
 649     printf("AliAnaPi0EbE::MakeAnalysisFillHistograms()  - No output pi0 in AOD branch with name < %s >,STOP \n",GetOutputAODName().Data());
 650     abort();
 651   }
 652   //Loop on stored AOD pi0
 653   Int_t naod = GetOutputAODBranch()->GetEntriesFast();
 654   if(GetDebug() > 0) printf("AliAnaPi0EbE::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod);
 655
 656   for(Int_t iaod = 0; iaod < naod ; iaod++){
 657
 658     AliAODPWG4Particle* pi0 =  (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
 659     Int_t pdg = pi0->GetPdg();
 660
 661     if(IsCaloPIDOn() && pdg != AliCaloPID::kPi0) continue;
 662
 663     //Fill pi0 histograms
 664     Float_t pt  = pi0->Pt();
 665     Float_t phi = pi0->Phi();
 666     if(phi < 0) phi+=TMath::TwoPi();
 667     Float_t eta = pi0->Eta();
 668
 669     fhPtPi0       ->Fill(pt);
 670     fhPtEtaPhiPi0 ->Fill(pt,eta,phi);
 671
 672     if(IsDataMC()){
 673       if((GetReader()->GetDataType() == AliCaloTrackReader::kMC && fAnaType!=kSSCalo) ||
 674          GetReader()->GetDataType() != AliCaloTrackReader::kMC){
 675         if(GetMCAnalysisUtils()->CheckTagBit(pi0->GetTag(), AliMCAnalysisUtils::kMCPi0)){
 676           fhPtMCPi0  ->Fill(pt);
 677           fhPhiMCPi0 ->Fill(pt,phi);
 678           fhEtaMCPi0 ->Fill(pt,eta);
 679         }
 680         else{
 681           fhPtMCNoPi0  ->Fill(pt);
 682           fhPhiMCNoPi0 ->Fill(pt,phi);
 683           fhEtaMCNoPi0 ->Fill(pt,eta);
 684         }
 685       }
 686     }//Histograms with MC
 687
 688   }// aod loop
 689
 690 }
 691
 692
 693 //____________________________________________________________________________
 694 void AliAnaPi0EbE::Init()
 695 {
 696   //Init
 697   //Do some checks
 698   if(fCalorimeter == "PHOS" && !GetReader()->IsPHOSSwitchedOn() && NewOutputAOD()){
 699     printf("AliAnaPi0EbE::Init() - !!STOP: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!\n");
 700     abort();
 701   }
 702   else  if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn() && NewOutputAOD()){
 703     printf("AliAnaPi0EbE::Init() - !!STOP: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!\n");
 704     abort();
 705   }
 706
 707 }
 708
 709 //____________________________________________________________________________
 710 void AliAnaPi0EbE::InitParameters()
 711 {
 712   //Initialize the parameters of the analysis.
 713   AddToHistogramsName("AnaPi0EbE_");
 714
 715   fInputAODGammaConvName = "gammaconv" ;
 716   fAnaType = kIMCalo ;
 717   fCalorimeter = "EMCAL" ;
 718   fMinDist  = 2.;
 719   fMinDist2 = 4.;
 720   fMinDist3 = 5.;
 721
 722 }
 723
 724 //__________________________________________________________________
 725 void AliAnaPi0EbE::Print(const Option_t * opt) const
 726 {
 727   //Print some relevant parameters set for the analysis
 728   if(! opt)
 729     return;
 730
 731   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
 732   AliAnaPartCorrBaseClass::Print("");
 733   printf("Analysis Type = %d \n",  fAnaType) ;
 734   if(fAnaType == kSSCalo){
 735     printf("Calorimeter            =     %s\n", fCalorimeter.Data()) ;
 736     printf("Min Distance to Bad Channel   = %2.1f\n",fMinDist);
 737     printf("Min Distance to Bad Channel 2 = %2.1f\n",fMinDist2);
 738     printf("Min Distance to Bad Channel 3 = %2.1f\n",fMinDist3);
 739   }
 740   printf("    \n") ;
 741
 742 }