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 is 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 CTS
  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), fhEPi0(0), fhEEtaPhiPi0(0),
  58 fhEDispPi0(0), fhEDispBkg(0),
  59 fhELambda0Pi0(0), fhELambda0Bkg(0),
  60 fhELambda1Pi0(0), fhELambda1Bkg(0),
  61 fhELambdaPi0EtaCen(0),fhELambdaPi0EtaMid(0),fhELambdaPi0EtaBor(0),
  62 fhClusterPairDiffTimeE(0),fhClusterPairDiffTimeAsy(0),
  63 //MC histos
  64 fhELambdaPhNoConv(0),fhELambdaConvPhotons(0),fhELambdaElectrons(0),
  65 fhELambdaPi0NoPh(0),fhELambdaOtherParts(0),
  66 fhPtMCNoPi0(0),fhPhiMCNoPi0(0),fhEtaMCNoPi0(0),
  67 fhPtMCPi0(0),fhPhiMCPi0(0),fhEtaMCPi0(0)
  68 {
  69   //default ctor
  70
  71   //Initialize parameters
  72   InitParameters();
  73
  74 }
  75
  76 //____________________________________________________________________________
  77 AliAnaPi0EbE::~AliAnaPi0EbE()
  78 {
  79   //dtor
  80   if(fInputAODGammaConv){
  81     fInputAODGammaConv->Clear() ;
  82     delete fInputAODGammaConv ;
  83   }
  84 }
  85
  86 //________________________________________________________________________
  87 TObjString *  AliAnaPi0EbE::GetAnalysisCuts()
  88 {
  89         //Save parameters used for analysis
  90          TString parList ; //this will be list of parameters used for this analysis.
  91    const Int_t buffersize = 255;
  92          char onePar[buffersize] ;
  93
  94          snprintf(onePar,buffersize,"--- AliAnaPi0EbE ---\n") ;
  95          parList+=onePar ;
  96          snprintf(onePar,buffersize,"fAnaType=%d (Pi0 selection type) \n",fAnaType) ;
  97          parList+=onePar ;
  98
  99          if(fAnaType == kSSCalo){
 100            snprintf(onePar,buffersize,"Calorimeter: %s\n",fCalorimeter.Data()) ;
 101            parList+=onePar ;
 102            snprintf(onePar,buffersize,"fMinDist =%2.2f (Minimal distance to bad channel to accept cluster) \n",fMinDist) ;
 103            parList+=onePar ;
 104            snprintf(onePar,buffersize,"fMinDist2=%2.2f (Cuts on Minimal distance to study acceptance evaluation) \n",fMinDist2) ;
 105            parList+=onePar ;
 106            snprintf(onePar,buffersize,"fMinDist3=%2.2f (One more cut on distance used for acceptance-efficiency study) \n",fMinDist3) ;
 107            parList+=onePar ;
 108          }
 109
 110          //Get parameters set in base class.
 111          parList += GetBaseParametersList() ;
 112
 113          //Get parameters set in PID class.
 114          if(fAnaType == kSSCalo) parList += GetCaloPID()->GetPIDParametersList() ;
 115
 116          return new TObjString(parList) ;
 117 }
 118
 119 //________________________________________________________________________
 120 TList *  AliAnaPi0EbE::GetCreateOutputObjects()
 121 {
 122   // Create histograms to be saved in output file and
 123   // store them in outputContainer
 124   TList * outputContainer = new TList() ;
 125   outputContainer->SetName("Pi0EbEHistos") ;
 126
 127   Int_t nptbins  = GetHistoPtBins();           Float_t ptmax  = GetHistoPtMax();           Float_t ptmin  = GetHistoPtMin();
 128   Int_t nphibins = GetHistoPhiBins();          Float_t phimax = GetHistoPhiMax();          Float_t phimin = GetHistoPhiMin();
 129   Int_t netabins = GetHistoEtaBins();          Float_t etamax = GetHistoEtaMax();          Float_t etamin = GetHistoEtaMin();
 130   Int_t ssbins   = GetHistoShowerShapeBins();  Float_t ssmax  = GetHistoShowerShapeMax();  Float_t ssmin  = GetHistoShowerShapeMin();
 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   fhEPi0  = new TH1F("hEPi0","Number of identified  #pi^{0} decay",nptbins,ptmin,ptmax);
 138   fhEPi0->SetYTitle("N");
 139   fhEPi0->SetXTitle("E  #pi^{0}(GeV)");
 140   outputContainer->Add(fhEPi0) ;
 141
 142   fhEEtaPhiPi0  = new TH3F
 143   ("hEEtaPhiPi0","Selected #pi^{0} pairs: E vs #eta vs #phi",nptbins,ptmin,ptmax,netabins,etamin,etamax, nphibins,phimin,phimax);
 144   fhEEtaPhiPi0->SetZTitle("#phi");
 145   fhEEtaPhiPi0->SetYTitle("#eta");
 146   fhEEtaPhiPi0->SetXTitle("E (GeV)");
 147   outputContainer->Add(fhEEtaPhiPi0) ;
 148
 149   ////////
 150
 151   if(fAnaType == kIMCalo){
 152
 153     fhEDispPi0  = new TH2F
 154     ("hEDispPi0","Selected #pi^{0} pairs: E vs dispersion",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 155     fhEDispPi0->SetYTitle("dispersion");
 156     fhEDispPi0->SetXTitle("E (GeV)");
 157     outputContainer->Add(fhEDispPi0) ;
 158
 159     fhEDispBkg  = new TH2F
 160     ("hEDispBkg","Rejected #pi^{0} pairs: E vs dispersion",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 161     fhEDispBkg->SetYTitle("dispersion");
 162     fhEDispBkg->SetXTitle("E (GeV)");
 163     outputContainer->Add(fhEDispBkg) ;
 164
 165     fhELambda0Pi0  = new TH2F
 166     ("hELambda0Pi0","Selected #pi^{0} pairs: E vs #lambda_{0}",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 167     fhELambda0Pi0->SetYTitle("#lambda_{0}");
 168     fhELambda0Pi0->SetXTitle("E (GeV)");
 169     outputContainer->Add(fhELambda0Pi0) ;
 170
 171     fhELambda0Bkg  = new TH2F
 172     ("hELambda0Bkg","Rejected #pi^{0} pairs: E vs #lambda_{0}",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 173     fhELambda0Bkg->SetYTitle("#lambda_{0}");
 174     fhELambda0Bkg->SetXTitle("E (GeV)");
 175     outputContainer->Add(fhELambda0Bkg) ;
 176
 177     fhELambda1Pi0  = new TH2F
 178     ("hELambda1Pi0","Selected #pi^{0} pairs: E vs #lambda_{1}",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 179     fhELambda1Pi0->SetYTitle("#lambda_{1}");
 180     fhELambda1Pi0->SetXTitle("E (GeV)");
 181     outputContainer->Add(fhELambda1Pi0) ;
 182
 183     fhELambda1Bkg  = new TH2F
 184     ("hELambda1Bkg","Rejected #pi^{0} pairs: E vs #lambda_{1}",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 185     fhELambda1Bkg->SetYTitle("#lambda_{1}");
 186     fhELambda1Bkg->SetXTitle("E (GeV)");
 187     outputContainer->Add(fhELambda1Bkg) ;
 188
 189     fhELambdaPi0EtaCen  = new TH2F
 190     ("hELambdaPi0EtaCen","Selected #pi^{0} pairs: E vs #lambda_{0}, |#eta < 0.2|",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 191     fhELambdaPi0EtaCen->SetYTitle("#lambda_{0}");
 192     fhELambdaPi0EtaCen->SetXTitle("E (GeV)");
 193     outputContainer->Add(fhELambdaPi0EtaCen) ;
 194
 195     fhELambdaPi0EtaMid  = new TH2F
 196     ("hELambdaPi0EtaMid","Selected #pi^{0} pairs: E vs #lambda_{0}, |0.2< #eta < 0.5|",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 197     fhELambdaPi0EtaMid->SetYTitle("#lambda_{0}");
 198     fhELambdaPi0EtaMid->SetXTitle("E (GeV)");
 199     outputContainer->Add(fhELambdaPi0EtaMid) ;
 200
 201     fhELambdaPi0EtaBor  = new TH2F
 202     ("hELambdaPi0EtaBor","Selected #pi^{0} pairs: E vs #lambda_{0}, |#eta > 0.5|",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 203     fhELambdaPi0EtaBor->SetYTitle("#lambda_{0}");
 204     fhELambdaPi0EtaBor->SetXTitle("E (GeV)");
 205     outputContainer->Add(fhELambdaPi0EtaBor) ;
 206
 207     fhClusterPairDiffTimeE = new TH2F("hClusterPairDiffTimeE","cluster pair time difference vs E",nptbins,ptmin,ptmax, 200,-100,100);
 208     fhClusterPairDiffTimeE->SetXTitle("E_{pair} (GeV)");
 209     fhClusterPairDiffTimeE->SetYTitle("#Delta t (ns)");
 210     outputContainer->Add(fhClusterPairDiffTimeE);
 211
 212     fhClusterPairDiffTimeAsy = new TH2F("hClusterPairDiffTime","cluster pair time difference vs pair asymmetry",100,0,1, 200,-100,100);
 213     fhClusterPairDiffTimeAsy->SetXTitle("Asymmetry");
 214     fhClusterPairDiffTimeAsy->SetYTitle("#Delta t (ns)");
 215     outputContainer->Add(fhClusterPairDiffTimeAsy);
 216
 217   }// Invariant mass analysis in calorimeters only
 218
 219   if(IsDataMC()) {
 220     if((GetReader()->GetDataType() == AliCaloTrackReader::kMC && fAnaType!=kSSCalo) ||
 221        GetReader()->GetDataType() != AliCaloTrackReader::kMC){
 222
 223       fhPtMCPi0  = new TH1F("hPtMCPi0","Identified pi0 from pi0",nptbins,ptmin,ptmax);
 224       fhPtMCPi0->SetYTitle("N");
 225       fhPtMCPi0->SetXTitle("p_{T #pi^{0}}(GeV/c)");
 226       outputContainer->Add(fhPtMCPi0) ;
 227
 228       fhPhiMCPi0  = new TH2F
 229       ("hPhiMCPi0","Identified pi0 from pi0",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 230       fhPhiMCPi0->SetYTitle("#phi");
 231       fhPhiMCPi0->SetXTitle("p_{T #pi^{0}} (GeV/c)");
 232       outputContainer->Add(fhPhiMCPi0) ;
 233
 234       fhEtaMCPi0  = new TH2F
 235       ("hEtaMCPi0","Identified pi0 from pi0",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 236       fhEtaMCPi0->SetYTitle("#eta");
 237       fhEtaMCPi0->SetXTitle("p_{T #pi^{0}} (GeV/c)");
 238       outputContainer->Add(fhEtaMCPi0) ;
 239
 240       fhPtMCNoPi0  = new TH1F("hPtMCNoPi0","Identified pi0 not from pi0",nptbins,ptmin,ptmax);
 241       fhPtMCNoPi0->SetYTitle("N");
 242       fhPtMCNoPi0->SetXTitle("p_{T #pi^{0}}(GeV/c)");
 243       outputContainer->Add(fhPtMCNoPi0) ;
 244
 245       fhPhiMCNoPi0  = new TH2F
 246       ("hPhiMCNoPi0","Identified pi0 not from pi0",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 247       fhPhiMCNoPi0->SetYTitle("#phi");
 248       fhPhiMCNoPi0->SetXTitle("p_{T #pi^{0}} (GeV/c)");
 249       outputContainer->Add(fhPhiMCNoPi0) ;
 250
 251       fhEtaMCNoPi0  = new TH2F
 252       ("hEtaMCNoPi0","Identified pi0 not from pi0",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 253       fhEtaMCNoPi0->SetYTitle("#eta");
 254       fhEtaMCNoPi0->SetXTitle("p_{T #pi^{0}} (GeV/c)");
 255       outputContainer->Add(fhEtaMCNoPi0) ;
 256
 257       if(fAnaType == kIMCalo){
 258
 259         fhELambdaPhNoConv  = new TH2F
 260         ("hELambdaPhNoConv","Selected #pi^{0} pairs (Really photons and not conversion): E vs #lambda_{0} vs #lambda_{1}",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 261         fhELambdaPhNoConv->SetZTitle("#lambda_{1}");
 262         fhELambdaPhNoConv->SetYTitle("#lambda_{0}");
 263         fhELambdaPhNoConv->SetXTitle("E (GeV)");
 264         outputContainer->Add(fhELambdaPhNoConv) ;
 265
 266
 267         fhELambdaConvPhotons  = new TH2F
 268         ("hELambdaConvPhotons","Selected #pi^{0} pairs (Converted photons): E vs #lambda_{0} vs #lambda_{1}",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 269         fhELambdaConvPhotons->SetZTitle("#lambda_{1}");
 270         fhELambdaConvPhotons->SetYTitle("#lambda_{0}");
 271         fhELambdaConvPhotons->SetXTitle("E (GeV)");
 272         outputContainer->Add(fhELambdaConvPhotons) ;
 273
 274         fhELambdaElectrons  = new TH2F
 275         ("hELambdaElectrons","Selected #pi^{0} pairs(Electrons ): E vs #lambda_{0} vs #lambda_{1}",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 276         fhELambdaElectrons->SetZTitle("#lambda_{1}");
 277         fhELambdaElectrons->SetYTitle("#lambda_{0}");
 278         fhELambdaElectrons->SetXTitle("E (GeV)");
 279         outputContainer->Add(fhELambdaElectrons) ;
 280
 281         fhELambdaPi0NoPh  = new TH2F
 282         ("hELambdaPi0NoPh","Selected #pi^{0} pairs(Pi0s instead of photons): E vs #lambda_{0} vs #lambda_{1}",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 283         fhELambdaPi0NoPh->SetZTitle("#lambda_{1}");
 284         fhELambdaPi0NoPh->SetYTitle("#lambda_{0}");
 285         fhELambdaPi0NoPh->SetXTitle("E (GeV)");
 286         outputContainer->Add(fhELambdaPi0NoPh) ;
 287
 288         fhELambdaOtherParts  = new TH2F
 289         ("hELambdaOtherParts","Selected #pi^{0} pairs (Other parts): E vs #lambda_{0} vs #lambda_{1}",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 290         fhELambdaOtherParts->SetZTitle("#lambda_{1}");
 291         fhELambdaOtherParts->SetYTitle("#lambda_{0}");
 292         fhELambdaOtherParts->SetXTitle("E (GeV)");
 293         outputContainer->Add(fhELambdaOtherParts) ;
 294
 295       }//kIMCalo
 296
 297     }
 298   }//Histos with MC
 299
 300
 301   //Keep neutral meson selection histograms if requiered
 302   //Setting done in AliNeutralMesonSelection
 303
 304   if(fAnaType!=kSSCalo && GetNeutralMesonSelection()){
 305
 306     TList * nmsHistos = GetNeutralMesonSelection()->GetCreateOutputObjects() ;
 307     if(GetNeutralMesonSelection()->AreNeutralMesonSelectionHistosKept())
 308       for(Int_t i = 0; i < nmsHistos->GetEntries(); i++) outputContainer->Add(nmsHistos->At(i)) ;
 309     delete nmsHistos;
 310
 311   }
 312
 313   return outputContainer ;
 314
 315 }
 316
 317 //__________________________________________________________________
 318 void  AliAnaPi0EbE::MakeAnalysisFillAOD()
 319 {
 320   //Do analysis and fill aods
 321
 322   switch(fAnaType)
 323     {
 324     case kIMCalo:
 325       MakeInvMassInCalorimeter();
 326       break;
 327
 328     case kSSCalo:
 329       MakeShowerShapeIdentification();
 330       break;
 331
 332     case kIMCaloTracks:
 333       MakeInvMassInCalorimeterAndCTS();
 334       break;
 335
 336     }
 337 }
 338
 339 //__________________________________________________________________
 340 void  AliAnaPi0EbE::MakeInvMassInCalorimeter()
 341 {
 342   //Do analysis and fill aods
 343   //Search for the photon decay in calorimeters
 344   //Read photon list from AOD, produced in class AliAnaPhoton
 345   //Check if 2 photons have the mass of the pi0.
 346
 347   TLorentzVector mom1;
 348   TLorentzVector mom2;
 349   TLorentzVector mom ;
 350   Int_t tag1 = 0;
 351   Int_t tag2 = 0;
 352   Int_t tag  = 0;
 353
 354   if(!GetInputAODBranch()){
 355     printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - No input calo photons in AOD with name branch < %s >, STOP \n",GetInputAODName().Data());
 356     abort();
 357   }
 358
 359   for(Int_t iphoton = 0; iphoton < GetInputAODBranch()->GetEntriesFast()-1; iphoton++){
 360     AliAODPWG4Particle * photon1 =  (AliAODPWG4Particle*) (GetInputAODBranch()->At(iphoton));
 361
 362     //Vertex cut in case of mixed events
 363     Int_t evtIndex1 = 0 ;
 364     if(GetMixedEvent())
 365       evtIndex1 = GetMixedEvent()->EventIndexForCaloCluster(photon1->GetCaloLabel(0)) ;
 366     if(TMath::Abs(GetVertex(evtIndex1)[2]) > GetZvertexCut()) continue ;  //vertex cut
 367     mom1 = *(photon1->Momentum());
 368
 369     //Get shower shape information of clusters
 370     TObjArray *clusters = 0;
 371     if     (fCalorimeter="EMCAL") clusters = GetEMCALClusters();
 372     else if(fCalorimeter="PHOS" ) clusters = GetPHOSClusters() ;
 373
 374     Bool_t bFound1        = kFALSE;
 375     Int_t  caloLabel1     = photon1->GetCaloLabel(0);
 376     Bool_t iclus1         = -1;
 377     AliVCluster *cluster1 = 0;
 378     if(clusters) {
 379       //Get original cluster, to recover some information
 380       for(Int_t iclus = 0; iclus < clusters->GetEntriesFast(); iclus++){
 381         AliVCluster *cluster= dynamic_cast<AliVCluster*> (clusters->At(iclus));
 382         if(cluster){
 383           if     (cluster->GetID()==caloLabel1) {
 384             bFound1  = kTRUE  ;
 385             cluster1 = cluster;
 386             iclus1   = iclus;
 387           }
 388         }
 389         if(bFound1) break;
 390       }// calorimeter clusters loop
 391     }
 392
 393     for(Int_t jphoton = iphoton+1; jphoton < GetInputAODBranch()->GetEntriesFast(); jphoton++){
 394
 395       AliAODPWG4Particle * photon2 =  (AliAODPWG4Particle*) (GetInputAODBranch()->At(jphoton));
 396       Int_t evtIndex2 = 0 ;
 397       if(GetMixedEvent())
 398         evtIndex2 = GetMixedEvent()->EventIndexForCaloCluster(photon2->GetCaloLabel(0)) ;
 399       if(GetMixedEvent() && (evtIndex1 == evtIndex2))
 400         continue ;
 401       if(TMath::Abs(GetVertex(evtIndex2)[2]) > GetZvertexCut()) continue ;  //vertex cut
 402       mom2 = *(photon2->Momentum());
 403
 404       //Photon1
 405       Float_t e1    = photon1->E();
 406       Float_t eta1  = photon1->Eta();
 407       Float_t tof1  = -1;
 408       Float_t disp1 = -1;
 409       Float_t l01   = -1;
 410       Float_t l11   = -1;
 411
 412       //Photon2
 413       Float_t e2    = photon2->E();
 414       Float_t eta2  = photon2->Eta();
 415       Float_t disp2 = -1;
 416       Float_t tof2  = -1;
 417       Float_t l02   = -1;
 418       Float_t l12   = -1;
 419
 420       Bool_t bFound2        = kFALSE;
 421       Int_t  caloLabel2     = photon2->GetCaloLabel(0);
 422       AliVCluster *cluster2 = 0;
 423       if(clusters){
 424         for(Int_t iclus = iclus1+1; iclus < clusters->GetEntriesFast(); iclus++){
 425           AliVCluster *cluster= dynamic_cast<AliVCluster*> (clusters->At(iclus));
 426           if(cluster){
 427             if(cluster->GetID()==caloLabel2) {
 428               bFound2  = kTRUE  ;
 429               cluster2 = cluster;
 430             }
 431           }
 432           if(bFound2) break;
 433         }// calorimeter clusters loop
 434
 435         //Photon/Cluster 1
 436         if(cluster1 && bFound1){
 437           disp1 = cluster1->GetDispersion();
 438           l11   = cluster1->GetM20();
 439           l01   = cluster1->GetM02();
 440           tof1  = cluster1->GetTOF()*1e9;
 441         }
 442 //        else printf("cluster1 not available: calo label %d / %d, cluster ID %d\n",
 443 //                     photon2->GetCaloLabel(0),(GetReader()->GetInputEvent())->GetNumberOfCaloClusters()-1,cluster1->GetID());
 444
 445         //Photon/Cluster 2
 446         if(cluster2 && bFound2){
 447           disp2 = cluster2->GetDispersion();
 448           l12   = cluster2->GetM20();
 449           l02   = cluster2->GetM02();
 450           tof2  = cluster2->GetTOF()*1e9;
 451         }
 452 //        else printf("cluster2 not available: calo label %d / %d, cluster ID %d\n",
 453 //                     photon2->GetCaloLabel(0),(GetReader()->GetInputEvent())->GetNumberOfCaloClusters()-1,cluster2->GetID());
 454
 455         //Select clusters with good time window difference
 456         Double_t t12diff = tof1-tof2;
 457         Float_t asymmetry = TMath::Abs(e1-e2)/(e1+e2);
 458         fhClusterPairDiffTimeE  ->Fill(e1+e2,    t12diff);
 459         fhClusterPairDiffTimeAsy->Fill(asymmetry,t12diff);
 460         if(TMath::Abs(t12diff) > GetPairTimeCut()) continue;
 461       }
 462
 463       //Select good pair (good phi, pt cuts, aperture and invariant mass)
 464       if(GetNeutralMesonSelection()->SelectPair(mom1, mom2))
 465       {
 466         if(GetDebug()>1)
 467           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());
 468
 469         //Play with the MC stack if available
 470         if(IsDataMC()){
 471           //Check origin of the candidates
 472           Int_t  label1 = photon1->GetLabel();
 473           Int_t  label2 = photon2->GetLabel();
 474           if(label1>=0)tag1 = GetMCAnalysisUtils()->CheckOrigin(label1, GetReader(), photon1->GetInputFileIndex());
 475           if(label2>=0)tag2 = GetMCAnalysisUtils()->CheckOrigin(label2, GetReader(), photon2->GetInputFileIndex());
 476
 477           if(GetDebug() > 0) printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - Origin of: photon1 %d; photon2 %d \n",tag1, tag2);
 478           if(GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCPi0Decay) &&
 479              GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCPi0Decay)){
 480
 481             //Check if pi0 mother is the same
 482             if(GetReader()->ReadStack()){
 483               if(label1>=0){
 484                 TParticle * mother1 = GetMCStack()->Particle(label1);//photon in kine tree
 485                 label1 = mother1->GetFirstMother();
 486                 //mother1 = GetMCStack()->Particle(label1);//pi0
 487               }
 488               if(label2>=0){
 489                 TParticle * mother2 = GetMCStack()->Particle(label2);//photon in kine tree
 490                 label2 = mother2->GetFirstMother();
 491                 //mother2 = GetMCStack()->Particle(label2);//pi0
 492                }
 493             }
 494             else if(GetReader()->ReadAODMCParticles()){//&& (input > -1)){
 495               if(label1>=0){
 496                 AliAODMCParticle * mother1 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles(photon1->GetInputFileIndex()))->At(label1);//photon in kine tree
 497                 label1 = mother1->GetMother();
 498                 //mother1 = GetMCStack()->Particle(label1);//pi0
 499                }
 500               if(label2>=0){
 501                 AliAODMCParticle * mother2 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles(photon2->GetInputFileIndex()))->At(label2);//photon in kine tree
 502                 label2 = mother2->GetMother();
 503                 //mother2 = GetMCStack()->Particle(label2);//pi0
 504               }
 505             }
 506
 507             //printf("mother1 %d, mother2 %d\n",label1,label2);
 508             if(label1 == label2 && label1>=0)
 509               GetMCAnalysisUtils()->SetTagBit(tag,AliMCAnalysisUtils::kMCPi0);
 510           }
 511         }//Work with stack also
 512
 513         //Fill some histograms about shower shape
 514         if(clusters && GetReader()->GetDataType()!=AliCaloTrackReader::kMC){
 515           //Photon1
 516
 517           //printf("Signal Cl1: e %f, pt %f, disp %f, l1 %f, l0 %f, eta %f, phi %f \n",
 518           //       e,pt,disp,l1,l0,photon2->Eta(),photon2->Phi());
 519
 520           fhEDispPi0   ->Fill(e1, disp1);
 521           fhELambda0Pi0->Fill(e1, l01  );
 522           fhELambda1Pi0->Fill(e1, l11  );
 523
 524           if(TMath::Abs(eta1)<0.2){
 525             fhELambdaPi0EtaCen->Fill(e1, l01);
 526           }
 527           else if (TMath::Abs(eta1)>0.2 && TMath::Abs(eta1)<0.5){
 528             fhELambdaPi0EtaMid->Fill(e1, l01);
 529           }
 530           else {
 531             fhELambdaPi0EtaBor->Fill(e1, l01);
 532           }
 533
 534           //Photon2
 535           //printf("Signal Cl2: e %f, pt %f, disp %f, l1 %f, l0 %f, eta %f, phi %f \n",e
 536           //     ,pt,disp,l1,l0,photon2->Eta(),photon2->Phi());
 537
 538           fhEDispPi0   ->Fill(e2, disp2);
 539           fhELambda0Pi0->Fill(e2, l02  );
 540           fhELambda1Pi0->Fill(e2, l12  );
 541
 542           if(TMath::Abs(eta2)<0.2){
 543             fhELambdaPi0EtaCen->Fill(e2, l02);
 544           }
 545           else if (TMath::Abs(eta2)>0.2 && TMath::Abs(eta2)<0.5){
 546             fhELambdaPi0EtaMid->Fill(e2, l02);
 547           }
 548           else {
 549             fhELambdaPi0EtaBor->Fill(e2, l02);
 550           }
 551
 552           if(IsDataMC()) {
 553             //Photon1
 554             if( GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCPhoton) && !GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCConversion) ){
 555               fhELambdaPhNoConv->Fill(e1, l01);
 556             }//photon   no conversion
 557             else if  ( GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCElectron)){
 558               fhELambdaElectrons->Fill(e1, l01);
 559             }//electron
 560             else if  ( GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCConversion) && GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCConversion) ){
 561               fhELambdaConvPhotons->Fill(e1, l01);
 562             }//convesion photon
 563
 564             else if  ( GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCPi0)  ){
 565               fhELambdaPi0NoPh->Fill(e1, l01);
 566             }//pi0
 567             else {
 568               fhELambdaOtherParts->Fill(e1, l01);
 569             }//other particules
 570
 571             //Photon 2
 572             if( GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCPhoton) && !GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCConversion) ){
 573               fhELambdaPhNoConv->Fill(e2, l02);
 574             }//photon   no conversion
 575             else if  ( GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCElectron)){
 576               fhELambdaElectrons->Fill(e2, l02);
 577             }//electron
 578             else if  ( GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCConversion) && GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCConversion) ){
 579               fhELambdaConvPhotons->Fill(e2, l02);
 580             }//convesion photon
 581
 582             else if  ( GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCPi0)  ){
 583               fhELambdaPi0NoPh->Fill(e2, l02);
 584             }//pi0
 585             else {
 586               fhELambdaOtherParts->Fill(e2, l02);
 587             }//other particules
 588           }//is datamc
 589
 590         }//MC histograms
 591
 592         //Create AOD for analysis
 593         mom = mom1+mom2;
 594         AliAODPWG4Particle pi0 = AliAODPWG4Particle(mom);
 595         //pi0.SetLabel(calo->GetLabel());
 596         pi0.SetIdentifiedParticleType(AliCaloPID::kPi0);
 597         pi0.SetDetector(photon1->GetDetector());
 598         pi0.SetTag(tag);
 599         //Set the indeces of the original caloclusters
 600         pi0.SetCaloLabel(photon1->GetCaloLabel(0), photon2->GetCaloLabel(0));
 601         //pi0.SetInputFileIndex(input);
 602         AddAODParticle(pi0);
 603       }//pi0
 604       else{
 605         Float_t phi = (mom1+mom2).Phi();
 606         if(phi < 0) phi+=TMath::TwoPi();
 607
 608         //Fill some histograms about shower shape
 609         if(GetReader()->GetDataType()!=AliCaloTrackReader::kMC){
 610           //Photon1
 611
 612           //printf("Bkg Cl1: e %f, pt %f, disp %f, l1 %f, l0 %f, eta %f, phi %f \n",
 613           //       e,pt,disp,l1,l0,photon2->Eta(),photon2->Phi());
 614
 615           fhEDispBkg   ->Fill(e1, disp1);
 616           fhELambda0Bkg->Fill(e1, l01  );
 617           fhELambda1Bkg->Fill(e1, l11  );
 618
 619           //Photon2
 620
 621           //printf("Bkg Cl2: e %f, pt %f, disp %f, l1 %f, l0 %f, eta %f, phi %f \n",
 622           //e,pt,disp,l1,l0,photon2->Eta(),photon2->Phi());
 623
 624           fhEDispBkg   ->Fill(e2, disp2);
 625           fhELambda0Bkg->Fill(e2, l02  );
 626           fhELambda1Bkg->Fill(e2, l12  );
 627
 628         }
 629
 630       }//bkg pair
 631
 632     }//2n photon loop
 633
 634   }//1st photon loop
 635
 636   if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - End fill AODs \n");
 637
 638 }
 639
 640 //__________________________________________________________________
 641 void  AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS()
 642 {
 643   //Do analysis and fill aods
 644   //Search for the photon decay in calorimeters
 645   //Read photon list from AOD, produced in class AliAnaPhoton and AliGammaConversion
 646   //Check if 2 photons have the mass of the pi0.
 647
 648   TLorentzVector mom1;
 649   TLorentzVector mom2;
 650   TLorentzVector mom ;
 651   Int_t tag1 = 0;
 652   Int_t tag2 = 0;
 653   Int_t tag  = 0;
 654   Int_t evtIndex = 0;
 655   if(!GetInputAODBranch()){
 656     printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - No input calo photons in AOD branch with name < %s > , STOP\n",GetInputAODName().Data());
 657     abort();
 658   }
 659
 660   for(Int_t iphoton = 0; iphoton < GetInputAODBranch()->GetEntriesFast(); iphoton++){
 661     AliAODPWG4Particle * photon1 =  (AliAODPWG4Particle*) (GetInputAODBranch()->At(iphoton));
 662     mom1 = *(photon1->Momentum());
 663
 664     //Play with the MC stack if available
 665     fInputAODGammaConv = (TClonesArray *) GetReader()->GetOutputEvent()->FindListObject(fInputAODGammaConvName);
 666     if(!fInputAODGammaConv) {
 667       printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - No input gamma conversions in AOD branch with name < %s >, STOP\n",fInputAODGammaConvName.Data());
 668       abort();
 669     }
 670     for(Int_t jphoton = 0; jphoton < fInputAODGammaConv->GetEntriesFast(); jphoton++){
 671       AliAODPWG4Particle * photon2 =  (AliAODPWG4Particle*) (fInputAODGammaConv->At(jphoton));
 672       if(GetMixedEvent())
 673         evtIndex = GetMixedEvent()->EventIndexForCaloCluster(photon2->GetCaloLabel(0)) ;
 674       if(TMath::Abs(GetVertex(evtIndex)[2]) > GetZvertexCut()) continue ;  //vertex cut
 675
 676       mom2 = *(photon2->Momentum());
 677
 678       //Int_t input = -1;       //if -1 photons come from different files, not a pi0
 679       //if(photon1->GetInputFileIndex() == photon2->GetInputFileIndex()) input = photon1->GetInputFileIndex();
 680
 681       //Select good pair (good phi, pt cuts, aperture and invariant mass)
 682       if(GetNeutralMesonSelection()->SelectPair(mom1, mom2)){
 683         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());
 684
 685         if(IsDataMC()){
 686           Int_t label1 = photon1->GetLabel();
 687           Int_t label2 = photon2->GetLabel();
 688           if(label1>=0)tag1 = GetMCAnalysisUtils()->CheckOrigin(label1, GetReader(), photon1->GetInputFileIndex());
 689           if(label2>=0)tag2 = GetMCAnalysisUtils()->CheckOrigin(label2, GetReader(), photon2->GetInputFileIndex());
 690           if(GetDebug() > 0) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - Origin of: photon1 %d; photon2 %d \n",tag1, tag2);
 691           if(GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCPi0Decay) &&
 692              GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCPi0Decay)){
 693             //Check if pi0 mother is the same
 694
 695             if(GetReader()->ReadStack()){
 696               if(label1>=0){
 697                 TParticle * mother1 = GetMCStack()->Particle(label1);//photon in kine tree
 698                 label1 = mother1->GetFirstMother();
 699                 //mother1 = GetMCStack()->Particle(label1);//pi0
 700               }
 701               if(label2>=0){
 702                 TParticle * mother2 = GetMCStack()->Particle(label2);//photon in kine tree
 703                 label2 = mother2->GetFirstMother();
 704                 //mother2 = GetMCStack()->Particle(label2);//pi0
 705               }
 706             }
 707             else if(GetReader()->ReadAODMCParticles()&& label1>=0 && label2>=0){ //&& (input > -1)){
 708               if(label1>=0){
 709                 AliAODMCParticle * mother1 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles(photon1->GetInputFileIndex()))->At(label1);//photon in kine tree
 710                 label1 = mother1->GetMother();
 711                 //mother1 = GetMCStack()->Particle(label1);//pi0
 712               }
 713               if(label2>=0){
 714                 AliAODMCParticle * mother2 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles(photon2->GetInputFileIndex()))->At(label2);//photon in kine tree
 715                 label2 = mother2->GetMother();
 716                 //mother2 = GetMCStack()->Particle(label2);//pi0
 717               }
 718             }
 719
 720             //printf("mother1 %d, mother2 %d\n",label1,label2);
 721             if(label1 == label2 && label1>=0)
 722               GetMCAnalysisUtils()->SetTagBit(tag,AliMCAnalysisUtils::kMCPi0);
 723           }
 724         }//Work with stack also
 725
 726         //Create AOD for analysis
 727         mom = mom1+mom2;
 728         AliAODPWG4Particle pi0 = AliAODPWG4Particle(mom);
 729         //pi0.SetLabel(calo->GetLabel());
 730         pi0.SetIdentifiedParticleType(AliCaloPID::kPi0);
 731         pi0.SetDetector(photon1->GetDetector());
 732         pi0.SetTag(tag);
 733         //Set the indeces of the original tracks or caloclusters
 734         pi0.SetCaloLabel(photon1->GetCaloLabel(0), -1);
 735         pi0.SetTrackLabel(photon2->GetTrackLabel(0), photon2->GetTrackLabel(1));
 736         //pi0.SetInputFileIndex(input);
 737         AddAODParticle(pi0);
 738       }//pi0
 739     }//2n photon loop
 740
 741   }//1st photon loop
 742
 743   if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - End fill AODs \n");
 744
 745 }
 746
 747
 748 //__________________________________________________________________
 749 void  AliAnaPi0EbE::MakeShowerShapeIdentification()
 750 {
 751   //Search for pi0 in fCalorimeter with shower shape analysis
 752
 753   TObjArray * pl = 0x0;
 754   //Select the Calorimeter of the photon
 755   if(fCalorimeter == "PHOS")
 756     pl = GetPHOSClusters();
 757   else if (fCalorimeter == "EMCAL")
 758     pl = GetEMCALClusters();
 759
 760   if(!pl) {
 761     Info("MakeShowerShapeIdentification","TObjArray with %s clusters is NULL!\n",fCalorimeter.Data());
 762     return;
 763   }
 764
 765   //Get vertex for photon momentum calculation
 766   //Double_t vertex2[] = {0,0,0} ; //vertex from second aod input
 767   //if(GetReader()->GetDataType() != AliCaloTrackReader::kMC)
 768   //{
 769   //if(GetReader()->GetSecondInputAODTree()) GetReader()->GetSecondInputAODVertex(vertex2);
 770   //}
 771
 772
 773   TLorentzVector mom ;
 774   for(Int_t icalo = 0; icalo < pl->GetEntriesFast(); icalo++){
 775     AliVCluster * calo = (AliVCluster*) (pl->At(icalo));
 776
 777     Int_t evtIndex = 0 ;
 778     if (GetMixedEvent()) {
 779       evtIndex=GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ;
 780     }
 781     if(TMath::Abs(GetVertex(evtIndex)[2]) > GetZvertexCut()) continue ;  //vertex cut
 782
 783     //Cluster selection, not charged, with pi0 id and in fiducial cut
 784
 785     //Input from second AOD?
 786     //Int_t input = 0;
 787     //  if     (fCalorimeter == "EMCAL" && GetReader()->GetEMCALClustersNormalInputEntries() <= icalo) input = 1 ;
 788     //  else if(fCalorimeter == "PHOS"  && GetReader()->GetPHOSClustersNormalInputEntries()  <= icalo) input = 1;
 789
 790     //Get Momentum vector,
 791     //if     (input == 0)
 792     if(GetReader()->GetDataType() != AliCaloTrackReader::kMC){
 793       calo->GetMomentum(mom,GetVertex(evtIndex)) ;}//Assume that come from vertex in straight line
 794     else{
 795       Double_t vertex[]={0,0,0};
 796       calo->GetMomentum(mom,vertex) ;
 797     }
 798     //else if(input == 1) calo->GetMomentum(mom,vertex2);//Assume that come from vertex in straight line
 799
 800     //If too small or big pt, skip it
 801     if(mom.Pt() < GetMinPt() || mom.Pt() > GetMaxPt() ) continue ;
 802     //Check acceptance selection
 803     if(IsFiducialCutOn()){
 804       Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter) ;
 805       if(! in ) continue ;
 806     }
 807
 808     //Create AOD for analysis
 809     AliAODPWG4Particle aodpi0 = AliAODPWG4Particle(mom);
 810     aodpi0.SetLabel(calo->GetLabel());
 811     //Set the indeces of the original caloclusters
 812     aodpi0.SetCaloLabel(calo->GetID(),-1);
 813     aodpi0.SetDetector(fCalorimeter);
 814     if(GetDebug() > 1)
 815       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());
 816
 817     //Check Distance to Bad channel, set bit.
 818     Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
 819     if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;
 820     if(distBad < fMinDist) //In bad channel (PHOS cristal size 2.2x2.2 cm)
 821       continue ;
 822
 823     if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - FillAOD: Bad channel cut passed %4.2f\n",distBad);
 824
 825     if(distBad > fMinDist3) aodpi0.SetDistToBad(2) ;
 826     else if(distBad > fMinDist2) aodpi0.SetDistToBad(1) ;
 827     else aodpi0.SetDistToBad(0) ;
 828
 829     //Check PID
 830     //PID selection or bit setting
 831     if(GetReader()->GetDataType() == AliCaloTrackReader::kMC){
 832       //Get most probable PID, check PID weights (in MC this option is mandatory)
 833       aodpi0.SetIdentifiedParticleType(GetCaloPID()->GetIdentifiedParticleType(fCalorimeter,calo->GetPID(),mom.E()));//PID with weights
 834       if(GetDebug() > 1)
 835         printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - FillAOD: PDG of identified particle %d\n",aodpi0.GetIdentifiedParticleType());
 836       //If primary is not pi0, skip it.
 837       if(aodpi0.GetIdentifiedParticleType() != AliCaloPID::kPi0) continue ;
 838     }
 839     else if(IsCaloPIDOn()){
 840       //Skip matched clusters with tracks
 841       if(IsTrackMatched(calo)) continue ;
 842
 843       //Get most probable PID, 2 options check PID weights
 844       //or redo PID, recommended option for EMCal.
 845       if(!IsCaloPIDRecalculationOn())
 846         aodpi0.SetIdentifiedParticleType(GetCaloPID()->GetIdentifiedParticleType(fCalorimeter,calo->GetPID(),mom.E()));//PID with weights
 847       else
 848         aodpi0.SetIdentifiedParticleType(GetCaloPID()->GetIdentifiedParticleType(fCalorimeter,mom,calo));//PID recalculated
 849
 850       if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - PDG of identified particle %d\n",aodpi0.GetIdentifiedParticleType());
 851
 852       //If cluster does not pass pid, not pi0, skip it.
 853       if(aodpi0.GetIdentifiedParticleType() != AliCaloPID::kPi0) continue ;
 854
 855     }
 856     else{
 857       //Set PID bits for later selection (AliAnaPi0 for example)
 858       //GetPDG already called in SetPIDBits.
 859       GetCaloPID()->SetPIDBits(fCalorimeter,calo,&aodpi0, GetCaloUtils());
 860       if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - PID Bits set \n");
 861     }
 862
 863     if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - Pi0 selection cuts passed: pT %3.2f, pdg %d\n",aodpi0.Pt(), aodpi0.GetIdentifiedParticleType());
 864
 865     //Play with the MC stack if available
 866     //Check origin of the candidates
 867     if(IsDataMC()){
 868       if((GetReader()->GetDataType() == AliCaloTrackReader::kMC && fAnaType!=kSSCalo) ||
 869          GetReader()->GetDataType() != AliCaloTrackReader::kMC){
 870         //aodpi0.SetInputFileIndex(input);
 871         Int_t tag       =0;
 872         tag = GetMCAnalysisUtils()->CheckOrigin(calo->GetLabel(),GetReader(), aodpi0.GetInputFileIndex());
 873         //GetMCAnalysisUtils()->CheckMultipleOrigin(calo->GetLabels(),calo->GetNLabels(), GetReader(), aodpi0.GetInputFileIndex(), tag);
 874         aodpi0.SetTag(tag);
 875         if(GetDebug() > 0) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - Origin of candidate %d\n",aodpi0.GetTag());
 876       }
 877     }//Work with stack also
 878
 879     //Add AOD with pi0 object to aod branch
 880     AddAODParticle(aodpi0);
 881
 882   }//loop
 883
 884   if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - End fill AODs \n");
 885
 886 }
 887 //__________________________________________________________________
 888 void  AliAnaPi0EbE::MakeAnalysisFillHistograms()
 889 {
 890   //Do analysis and fill histograms
 891
 892   if(!GetOutputAODBranch()){
 893     printf("AliAnaPi0EbE::MakeAnalysisFillHistograms()  - No output pi0 in AOD branch with name < %s >,STOP \n",GetOutputAODName().Data());
 894     abort();
 895   }
 896   //Loop on stored AOD pi0
 897   Int_t naod = GetOutputAODBranch()->GetEntriesFast();
 898   if(GetDebug() > 0) printf("AliAnaPi0EbE::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod);
 899
 900   for(Int_t iaod = 0; iaod < naod ; iaod++){
 901
 902     AliAODPWG4Particle* pi0 =  (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
 903     Int_t pdg = pi0->GetIdentifiedParticleType();
 904
 905     if(IsCaloPIDOn() && pdg != AliCaloPID::kPi0) continue;
 906
 907     //Fill pi0 histograms
 908     Float_t ener  = pi0->E();
 909     Float_t pt    = pi0->Pt();
 910     Float_t phi   = pi0->Phi();
 911     if(phi < 0) phi+=TMath::TwoPi();
 912     Float_t eta = pi0->Eta();
 913
 914     fhPtPi0      ->Fill(pt);
 915     fhEPi0       ->Fill(ener);
 916     fhEEtaPhiPi0 ->Fill(ener,eta,phi);
 917
 918     if(IsDataMC()){
 919       if((GetReader()->GetDataType() == AliCaloTrackReader::kMC && fAnaType!=kSSCalo) ||
 920          GetReader()->GetDataType() != AliCaloTrackReader::kMC){
 921         if(GetMCAnalysisUtils()->CheckTagBit(pi0->GetTag(), AliMCAnalysisUtils::kMCPi0)){
 922           fhPtMCPi0  ->Fill(pt);
 923           fhPhiMCPi0 ->Fill(pt,phi);
 924           fhEtaMCPi0 ->Fill(pt,eta);
 925         }
 926         else{
 927           fhPtMCNoPi0  ->Fill(pt);
 928           fhPhiMCNoPi0 ->Fill(pt,phi);
 929           fhEtaMCNoPi0 ->Fill(pt,eta);
 930         }
 931       }
 932     }//Histograms with MC
 933
 934   }// aod loop
 935
 936 }
 937
 938
 939 //____________________________________________________________________________
 940 void AliAnaPi0EbE::Init()
 941 {
 942   //Init
 943   //Do some checks
 944   if(fCalorimeter == "PHOS" && !GetReader()->IsPHOSSwitchedOn() && NewOutputAOD()){
 945     printf("AliAnaPi0EbE::Init() - !!STOP: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!\n");
 946     abort();
 947   }
 948   else  if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn() && NewOutputAOD()){
 949     printf("AliAnaPi0EbE::Init() - !!STOP: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!\n");
 950     abort();
 951   }
 952
 953 }
 954
 955 //____________________________________________________________________________
 956 void AliAnaPi0EbE::InitParameters()
 957 {
 958   //Initialize the parameters of the analysis.
 959   AddToHistogramsName("AnaPi0EbE_");
 960
 961   fInputAODGammaConvName = "gammaconv" ;
 962   fAnaType = kIMCalo ;
 963   fCalorimeter = "EMCAL" ;
 964   fMinDist  = 2.;
 965   fMinDist2 = 4.;
 966   fMinDist3 = 5.;
 967
 968 }
 969
 970 //__________________________________________________________________
 971 void AliAnaPi0EbE::Print(const Option_t * opt) const
 972 {
 973   //Print some relevant parameters set for the analysis
 974   if(! opt)
 975     return;
 976
 977   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
 978   AliAnaPartCorrBaseClass::Print("");
 979   printf("Analysis Type = %d \n",  fAnaType) ;
 980   if(fAnaType == kSSCalo){
 981     printf("Calorimeter            =     %s\n", fCalorimeter.Data()) ;
 982     printf("Min Distance to Bad Channel   = %2.1f\n",fMinDist);
 983     printf("Min Distance to Bad Channel 2 = %2.1f\n",fMinDist2);
 984     printf("Min Distance to Bad Channel 3 = %2.1f\n",fMinDist3);
 985   }
 986   printf("    \n") ;
 987
 988 }