PWG4/PartCorrDep/AliAnaPhoton.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 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: AliAnaPhoton.cxx 28688 2008-09-11 15:04:07Z gconesab $ */
  16
  17 //_________________________________________________________________________
  18 //
  19 // Class for the photon identification.
  20 // Clusters from calorimeters are identified as photons
  21 // and kept in the AOD. Few histograms produced.
  22 // Produces input for other analysis classes like AliAnaPi0,
  23 // AliAnaParticleHadronCorrelation ...
  24 //
  25 // -- Author: Gustavo Conesa (LNF-INFN)
  26 //////////////////////////////////////////////////////////////////////////////
  27
  28
  29 // --- ROOT system ---
  30 #include <TH2F.h>
  31 #include <TH3D.h>
  32 #include <TClonesArray.h>
  33 #include <TObjString.h>
  34 #include "TParticle.h"
  35 #include "TDatabasePDG.h"
  36
  37 // --- Analysis system ---
  38 #include "AliAnaPhoton.h"
  39 #include "AliCaloTrackReader.h"
  40 #include "AliStack.h"
  41 #include "AliCaloPID.h"
  42 #include "AliMCAnalysisUtils.h"
  43 #include "AliFiducialCut.h"
  44 #include "AliVCluster.h"
  45 #include "AliAODMCParticle.h"
  46 #include "AliMixedEvent.h"
  47
  48
  49 ClassImp(AliAnaPhoton)
  50
  51 //____________________________
  52 AliAnaPhoton::AliAnaPhoton() :
  53     AliAnaPartCorrBaseClass(),    fCalorimeter(""),
  54     fMinDist(0.),                 fMinDist2(0.),                fMinDist3(0.),
  55     fRejectTrackMatch(0),         fTimeCutMin(-10000),          fTimeCutMax(10000),
  56     fNCellsCut(0),                fFillSSHistograms(kFALSE),
  57     fNOriginHistograms(8),        fNPrimaryHistograms(4),
  58
  59     // Histograms
  60     fhNCellsE(0),                 fhMaxCellDiffClusterE(0),     fhTimeE(0), // Control histograms
  61     fhEPhoton(0),                 fhPtPhoton(0),
  62     fhPhiPhoton(0),               fhEtaPhoton(0),
  63     fhEtaPhiPhoton(0),            fhEtaPhi05Photon(0),
  64
  65     // Shower shape histograms
  66     fhDispE(0),                   fhLam0E(0),                   fhLam1E(0),
  67     fhDispETRD(0),                fhLam0ETRD(0),                fhLam1ETRD(0),
  68
  69     fhNCellsLam0LowE(0),          fhNCellsLam1LowE(0),          fhNCellsDispLowE(0),
  70     fhNCellsLam0HighE(0),         fhNCellsLam1HighE(0),         fhNCellsDispHighE(0),
  71
  72     fhEtaLam0LowE(0),             fhPhiLam0LowE(0),
  73     fhEtaLam0HighE(0),            fhPhiLam0HighE(0),
  74     fhLam0DispLowE(0),            fhLam0DispHighE(0),
  75     fhLam1Lam0LowE(0),            fhLam1Lam0HighE(0),
  76     fhDispLam1LowE(0),            fhDispLam1HighE(0),
  77
  78     // MC histograms
  79     fhMCPhotonELambda0NoOverlap(0),     fhMCPhotonELambda0TwoOverlap(0),      fhMCPhotonELambda0NOverlap(0),
  80     //Embedding
  81     fhEmbeddedSignalFractionEnergy(0),
  82     fhEmbedPhotonELambda0FullSignal(0), fhEmbedPhotonELambda0MostlySignal(0),
  83     fhEmbedPhotonELambda0MostlyBkg(0),  fhEmbedPhotonELambda0FullBkg(0),
  84     fhEmbedPi0ELambda0FullSignal(0),    fhEmbedPi0ELambda0MostlySignal(0),
  85     fhEmbedPi0ELambda0MostlyBkg(0),     fhEmbedPi0ELambda0FullBkg(0)
  86 {
  87   //default ctor
  88
  89   for(Int_t i = 0; i < 14; i++){
  90     fhMCPt     [i] = 0;
  91     fhMCE      [i] = 0;
  92     fhMCPhi    [i] = 0;
  93     fhMCEta    [i] = 0;
  94     fhMCDeltaE [i] = 0;
  95     fhMCDeltaPt[i] = 0;
  96     fhMC2E     [i] = 0;
  97     fhMC2Pt    [i] = 0;
  98
  99   }
 100
 101   for(Int_t i = 0; i < 7; i++){
 102     fhPtPrimMC [i] = 0;
 103     fhEPrimMC  [i] = 0;
 104     fhPhiPrimMC[i] = 0;
 105     fhYPrimMC  [i] = 0;
 106
 107     fhPtPrimMCAcc [i] = 0;
 108     fhEPrimMCAcc  [i] = 0;
 109     fhPhiPrimMCAcc[i] = 0;
 110     fhYPrimMCAcc  [i] = 0;
 111   }
 112
 113   for(Int_t i = 0; i < 6; i++){
 114     fhMCELambda0    [i]                  = 0;
 115     fhMCELambda1    [i]                  = 0;
 116     fhMCEDispersion [i]                  = 0;
 117     fhMCNCellsE     [i]                  = 0;
 118     fhMCMaxCellDiffClusterE[i]           = 0;
 119     fhMCLambda0vsClusterMaxCellDiffE0[i] = 0;
 120     fhMCLambda0vsClusterMaxCellDiffE2[i] = 0;
 121     fhMCLambda0vsClusterMaxCellDiffE6[i] = 0;
 122     fhMCNCellsvsClusterMaxCellDiffE0 [i] = 0;
 123     fhMCNCellsvsClusterMaxCellDiffE2 [i] = 0;
 124     fhMCNCellsvsClusterMaxCellDiffE6 [i] = 0;
 125   }
 126
 127   //Initialize parameters
 128   InitParameters();
 129
 130 }
 131
 132 //__________________________________________________________________________
 133 Bool_t  AliAnaPhoton::ClusterSelected(AliVCluster* calo, TLorentzVector mom)
 134 {
 135   //Select clusters if they pass different cuts
 136   if(GetDebug() > 2)
 137     printf("AliAnaPhoton::ClusterSelected() Current Event %d; Before selection : E %2.2f, pT %2.2f, Ecl %2.2f, phi %2.2f, eta %2.2f\n",
 138            GetReader()->GetEventNumber(),
 139            mom.E(), mom.Pt(),calo->E(),mom.Phi()*TMath::RadToDeg(),mom.Eta());
 140
 141   //.......................................
 142   //If too small or big energy, skip it
 143   if(mom.E() < GetMinEnergy() || mom.E() > GetMaxEnergy() ) return kFALSE ;
 144   if(GetDebug() > 2) printf("\t Cluster %d Pass E Cut \n",calo->GetID());
 145
 146   //.......................................
 147   // TOF cut, BE CAREFUL WITH THIS CUT
 148   Double_t tof = calo->GetTOF()*1e9;
 149   if(tof < fTimeCutMin || tof > fTimeCutMax) return kFALSE;
 150   if(GetDebug() > 2)  printf("\t Cluster %d Pass Time Cut \n",calo->GetID());
 151
 152   //.......................................
 153   if(calo->GetNCells() <= fNCellsCut && GetReader()->GetDataType() != AliCaloTrackReader::kMC) return kFALSE;
 154   if(GetDebug() > 2) printf("\t Cluster %d Pass NCell Cut \n",calo->GetID());
 155
 156   //.......................................
 157   //Check acceptance selection
 158   if(IsFiducialCutOn()){
 159     Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter) ;
 160     if(! in ) return kFALSE ;
 161   }
 162   if(GetDebug() > 2) printf("Fiducial cut passed \n");
 163
 164   //.......................................
 165   //Skip matched clusters with tracks
 166   if(fRejectTrackMatch){
 167     if(IsTrackMatched(calo,GetReader()->GetInputEvent())) {
 168       if(GetDebug() > 2) printf("\t Reject track-matched clusters\n");
 169       return kFALSE ;
 170     }
 171     else
 172       if(GetDebug() > 2)  printf(" Track-matching cut passed \n");
 173   }// reject matched clusters
 174
 175   //.......................................
 176   //Check Distance to Bad channel, set bit.
 177   Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
 178   if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;
 179   if(distBad < fMinDist) {//In bad channel (PHOS cristal size 2.2x2.2 cm), EMCAL ( cell units )
 180     return kFALSE ;
 181   }
 182   else if(GetDebug() > 2) printf("\t Bad channel cut passed %4.2f > %2.2f \n",distBad, fMinDist);
 183
 184   if(GetDebug() > 0)
 185     printf("AliAnaPhoton::ClusterSelected() Current Event %d; After  selection : E %2.2f, pT %2.2f, Ecl %2.2f, phi %2.2f, eta %2.2f\n",
 186            GetReader()->GetEventNumber(),
 187            mom.E(), mom.Pt(),calo->E(),mom.Phi()*TMath::RadToDeg(),mom.Eta());
 188
 189   //All checks passed, cluster selected
 190   return kTRUE;
 191
 192 }
 193
 194 //_____________________________________________________________
 195 void AliAnaPhoton::FillAcceptanceHistograms(){
 196   //Fill acceptance histograms if MC data is available
 197
 198   if(GetReader()->ReadStack()){
 199     AliStack * stack = GetMCStack();
 200     if(stack){
 201       for(Int_t i=0 ; i<stack->GetNtrack(); i++){
 202         TParticle * prim = stack->Particle(i) ;
 203         Int_t pdg = prim->GetPdgCode();
 204         //printf("i %d, %s %d  %s %d \n",i, stack->Particle(i)->GetName(), stack->Particle(i)->GetPdgCode(),
 205         //                             prim->GetName(), prim->GetPdgCode());
 206
 207         if(pdg == 22){
 208
 209           // Get tag of this particle photon from fragmentation, decay, prompt ...
 210           Int_t tag = GetMCAnalysisUtils()->CheckOrigin(i,GetReader(), 0);
 211           if(!GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton)){
 212             //A conversion photon from a hadron, skip this kind of photon
 213             //            printf("AliAnaPhoton::FillAcceptanceHistograms() - not a photon, weird!: tag %d, conv %d, pi0 %d, hadron %d, electron %d, unk %d, muon %d,pion %d, proton %d, neutron %d, kaon %d, antiproton %d, antineutron %d\n",tag,
 214             //                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion),
 215             //                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0),
 216             //                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOther),
 217             //                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron),
 218             //                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCUnknown),
 219             //                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCMuon),
 220             //                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPion),
 221             //                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCProton),
 222             //                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiNeutron),
 223             //                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCKaon),
 224             //                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiProton),
 225             //                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiNeutron));
 226
 227             return;
 228           }
 229
 230           //Get photon kinematics
 231           if(prim->Energy() == TMath::Abs(prim->Pz()))  continue ; //Protection against floating point exception
 232
 233           Double_t photonY   = 0.5*TMath::Log((prim->Energy()-prim->Pz())/(prim->Energy()+prim->Pz())) ;
 234           Double_t photonE   = prim->Energy() ;
 235           Double_t photonPt  = prim->Pt() ;
 236           Double_t photonPhi = TMath::RadToDeg()*prim->Phi() ;
 237           if(photonPhi < 0) photonPhi+=TMath::TwoPi();
 238           Double_t photonEta = prim->Eta() ;
 239
 240
 241           //Check if photons hit the Calorimeter
 242           TLorentzVector lv;
 243           prim->Momentum(lv);
 244           Bool_t inacceptance = kFALSE;
 245           if     (fCalorimeter == "PHOS"){
 246             if(GetPHOSGeometry() && GetCaloUtils()->IsPHOSGeoMatrixSet()){
 247               Int_t mod ;
 248               Double_t x,z ;
 249               if(GetPHOSGeometry()->ImpactOnEmc(prim,mod,z,x))
 250                 inacceptance = kTRUE;
 251               if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
 252             }
 253             else{
 254               if(GetFiducialCut()->IsInFiducialCut(lv,fCalorimeter))
 255                 inacceptance = kTRUE ;
 256               if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
 257             }
 258           }
 259           else if(fCalorimeter == "EMCAL" && GetCaloUtils()->IsEMCALGeoMatrixSet()){
 260             if(GetEMCALGeometry()){
 261
 262               Int_t absID=0;
 263
 264               GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(prim->Eta(),prim->Phi(),absID);
 265
 266               if( absID >= 0)
 267                 inacceptance = kTRUE;
 268
 269               //                  if(GetEMCALGeometry()->Impact(phot1) && GetEMCALGeometry()->Impact(phot2))
 270               //                    inacceptance = kTRUE;
 271               if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
 272             }
 273             else{
 274               if(GetFiducialCut()->IsInFiducialCut(lv,fCalorimeter))
 275                 inacceptance = kTRUE ;
 276               if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
 277             }
 278           }       //In EMCAL
 279
 280           //Fill histograms
 281
 282           fhYPrimMC[mcPPhoton]->Fill(photonPt, photonY) ;
 283           if(TMath::Abs(photonY) < 1.0)
 284           {
 285             fhEPrimMC  [mcPPhoton]->Fill(photonE ) ;
 286             fhPtPrimMC [mcPPhoton]->Fill(photonPt) ;
 287             fhPhiPrimMC[mcPPhoton]->Fill(photonE , photonPhi) ;
 288             fhYPrimMC[mcPPhoton]  ->Fill(photonE , photonEta) ;
 289           }
 290           if(inacceptance){
 291             fhEPrimMCAcc[mcPPhoton]  ->Fill(photonE ) ;
 292             fhPtPrimMCAcc[mcPPhoton] ->Fill(photonPt) ;
 293             fhPhiPrimMCAcc[mcPPhoton]->Fill(photonE , photonPhi) ;
 294             fhYPrimMCAcc[mcPPhoton]  ->Fill(photonE , photonY) ;
 295           }//Accepted
 296
 297           //Origin of photon
 298           if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt) && fhEPrimMC[mcPPrompt])
 299           {
 300             fhYPrimMC[mcPPrompt]->Fill(photonPt, photonY) ;
 301             if(TMath::Abs(photonY) < 1.0){
 302               fhEPrimMC  [mcPPrompt]->Fill(photonE ) ;
 303               fhPtPrimMC [mcPPrompt]->Fill(photonPt) ;
 304               fhPhiPrimMC[mcPPrompt]->Fill(photonE , photonPhi) ;
 305               fhYPrimMC[mcPPrompt]  ->Fill(photonE , photonEta) ;
 306             }
 307             if(inacceptance){
 308               fhEPrimMCAcc[mcPPrompt]  ->Fill(photonE ) ;
 309               fhPtPrimMCAcc[mcPPrompt] ->Fill(photonPt) ;
 310               fhPhiPrimMCAcc[mcPPrompt]->Fill(photonE , photonPhi) ;
 311               fhYPrimMCAcc[mcPPrompt]  ->Fill(photonE , photonY) ;
 312             }//Accepted
 313           }
 314           else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation) && fhEPrimMC[mcPFragmentation])
 315           {
 316             fhYPrimMC[mcPFragmentation]->Fill(photonPt, photonY) ;
 317             if(TMath::Abs(photonY) < 1.0){
 318               fhEPrimMC  [mcPFragmentation]->Fill(photonE ) ;
 319               fhPtPrimMC [mcPFragmentation]->Fill(photonPt) ;
 320               fhPhiPrimMC[mcPFragmentation]->Fill(photonE , photonPhi) ;
 321               fhYPrimMC[mcPFragmentation]  ->Fill(photonE , photonEta) ;
 322             }
 323             if(inacceptance){
 324               fhEPrimMCAcc[mcPFragmentation]  ->Fill(photonE ) ;
 325               fhPtPrimMCAcc[mcPFragmentation] ->Fill(photonPt) ;
 326               fhPhiPrimMCAcc[mcPFragmentation]->Fill(photonE , photonPhi) ;
 327               fhYPrimMCAcc[mcPFragmentation]  ->Fill(photonE , photonY) ;
 328             }//Accepted
 329           }
 330           else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCISR) && fhEPrimMC[mcPISR])
 331           {
 332             fhYPrimMC[mcPISR]->Fill(photonPt, photonY) ;
 333             if(TMath::Abs(photonY) < 1.0){
 334               fhEPrimMC  [mcPISR]->Fill(photonE ) ;
 335               fhPtPrimMC [mcPISR]->Fill(photonPt) ;
 336               fhPhiPrimMC[mcPISR]->Fill(photonE , photonPhi) ;
 337               fhYPrimMC[mcPISR]->Fill(photonE , photonEta) ;
 338             }
 339             if(inacceptance){
 340               fhEPrimMCAcc[mcPISR]  ->Fill(photonE ) ;
 341               fhPtPrimMCAcc[mcPISR] ->Fill(photonPt) ;
 342               fhPhiPrimMCAcc[mcPISR]->Fill(photonE , photonPhi) ;
 343               fhYPrimMCAcc[mcPISR]  ->Fill(photonE , photonY) ;
 344             }//Accepted
 345           }
 346           else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay)&& fhEPrimMC[mcPPi0Decay])
 347           {
 348             fhYPrimMC[mcPPi0Decay]->Fill(photonPt, photonY) ;
 349             if(TMath::Abs(photonY) < 1.0){
 350               fhEPrimMC  [mcPPi0Decay]->Fill(photonE ) ;
 351               fhPtPrimMC [mcPPi0Decay]->Fill(photonPt) ;
 352               fhPhiPrimMC[mcPPi0Decay]->Fill(photonE , photonPhi) ;
 353               fhYPrimMC[mcPPi0Decay]  ->Fill(photonE , photonEta) ;
 354             }
 355             if(inacceptance){
 356               fhEPrimMCAcc[mcPPi0Decay]  ->Fill(photonE ) ;
 357               fhPtPrimMCAcc[mcPPi0Decay] ->Fill(photonPt) ;
 358               fhPhiPrimMCAcc[mcPPi0Decay]->Fill(photonE , photonPhi) ;
 359               fhYPrimMCAcc[mcPPi0Decay]  ->Fill(photonE , photonY) ;
 360             }//Accepted
 361           }
 362           else if( (GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay) ||
 363                     GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay)) && fhEPrimMC[mcPOtherDecay])
 364           {
 365             fhYPrimMC[mcPOtherDecay]->Fill(photonPt, photonY) ;
 366             if(TMath::Abs(photonY) < 1.0){
 367               fhEPrimMC  [mcPOtherDecay]->Fill(photonE ) ;
 368               fhPtPrimMC [mcPOtherDecay]->Fill(photonPt) ;
 369               fhPhiPrimMC[mcPOtherDecay]->Fill(photonE , photonPhi) ;
 370               fhYPrimMC[mcPOtherDecay]  ->Fill(photonE , photonEta) ;
 371             }
 372             if(inacceptance){
 373               fhEPrimMCAcc[mcPOtherDecay]  ->Fill(photonE ) ;
 374               fhPtPrimMCAcc[mcPOtherDecay] ->Fill(photonPt) ;
 375               fhPhiPrimMCAcc[mcPOtherDecay]->Fill(photonE , photonPhi) ;
 376               fhYPrimMCAcc[mcPOtherDecay]  ->Fill(photonE , photonY) ;
 377             }//Accepted
 378           }
 379           else if(fhEPrimMC[mcPOther])
 380           {
 381             fhYPrimMC[mcPOther]->Fill(photonPt, photonY) ;
 382             if(TMath::Abs(photonY) < 1.0){
 383               fhEPrimMC  [mcPOther]->Fill(photonE ) ;
 384               fhPtPrimMC [mcPOther]->Fill(photonPt) ;
 385               fhPhiPrimMC[mcPOther]->Fill(photonE , photonPhi) ;
 386               fhYPrimMC[mcPOther]  ->Fill(photonE , photonEta) ;
 387             }
 388             if(inacceptance){
 389               fhEPrimMCAcc[mcPOther]  ->Fill(photonE ) ;
 390               fhPtPrimMCAcc[mcPOther] ->Fill(photonPt) ;
 391               fhPhiPrimMCAcc[mcPOther]->Fill(photonE , photonPhi) ;
 392               fhYPrimMCAcc[mcPOther]  ->Fill(photonE , photonY) ;
 393             }//Accepted
 394           }//Other origin
 395         }// Primary photon
 396       }//loop on primaries
 397     }//stack exists and data is MC
 398   }//read stack
 399   else if(GetReader()->ReadAODMCParticles()){
 400     TClonesArray * mcparticles = GetReader()->GetAODMCParticles(0);
 401     if(mcparticles){
 402       Int_t nprim = mcparticles->GetEntriesFast();
 403
 404       for(Int_t i=0; i < nprim; i++)
 405       {
 406         AliAODMCParticle * prim = (AliAODMCParticle *) mcparticles->At(i);
 407
 408         Int_t pdg = prim->GetPdgCode();
 409
 410         if(pdg == 22){
 411
 412           // Get tag of this particle photon from fragmentation, decay, prompt ...
 413           Int_t tag = GetMCAnalysisUtils()->CheckOrigin(i,GetReader(), 0);
 414           if(!GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton)){
 415             //A conversion photon from a hadron, skip this kind of photon
 416 //            printf("AliAnaPhoton::FillAcceptanceHistograms() - not a photon, weird!: tag %d, conv %d, pi0 %d, hadron %d, electron %d, unk %d, muon %d,pion %d, proton %d, neutron %d, kaon %d, antiproton %d, antineutron %d\n",tag,
 417 //                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion),
 418 //                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0),
 419 //                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOther),
 420 //                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron),
 421 //                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCUnknown),
 422 //                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCMuon),
 423 //                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPion),
 424 //                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCProton),
 425 //                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiNeutron),
 426 //                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCKaon),
 427 //                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiProton),
 428 //                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiNeutron));
 429
 430             return;
 431           }
 432
 433           //Get photon kinematics
 434           if(prim->E() == TMath::Abs(prim->Pz()))  continue ; //Protection against floating point exception
 435
 436           Double_t photonY  = 0.5*TMath::Log((prim->E()-prim->Pz())/(prim->E()+prim->Pz())) ;
 437           Double_t photonE   = prim->E() ;
 438           Double_t photonPt  = prim->Pt() ;
 439           Double_t photonPhi = TMath::RadToDeg()*prim->Phi() ;
 440           if(photonPhi < 0) photonPhi+=TMath::TwoPi();
 441           Double_t photonEta = prim->Eta() ;
 442
 443           //Check if photons hit the Calorimeter
 444           TLorentzVector lv;
 445           lv.SetPxPyPzE(prim->Px(),prim->Py(),prim->Pz(),prim->E());
 446           Bool_t inacceptance = kFALSE;
 447           if     (fCalorimeter == "PHOS"){
 448             if(GetPHOSGeometry() && GetCaloUtils()->IsPHOSGeoMatrixSet()){
 449               Int_t mod ;
 450               Double_t x,z ;
 451               Double_t vtx[]={prim->Xv(),prim->Yv(),prim->Zv()};
 452               if(GetPHOSGeometry()->ImpactOnEmc(vtx, prim->Theta(),prim->Phi(),mod,z,x))
 453                 inacceptance = kTRUE;
 454               if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
 455             }
 456             else{
 457               if(GetFiducialCut()->IsInFiducialCut(lv,fCalorimeter))
 458                 inacceptance = kTRUE ;
 459               if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
 460             }
 461           }
 462           else if(fCalorimeter == "EMCAL" && GetCaloUtils()->IsEMCALGeoMatrixSet()){
 463             if(GetEMCALGeometry()){
 464
 465               Int_t absID=0;
 466
 467               GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(prim->Eta(),prim->Phi(),absID);
 468
 469               if( absID >= 0)
 470                 inacceptance = kTRUE;
 471
 472               if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
 473             }
 474             else{
 475               if(GetFiducialCut()->IsInFiducialCut(lv,fCalorimeter))
 476                 inacceptance = kTRUE ;
 477               if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
 478             }
 479           }       //In EMCAL
 480
 481           //Fill histograms
 482
 483           fhYPrimMC[mcPPhoton]->Fill(photonPt, photonY) ;
 484           if(TMath::Abs(photonY) < 1.0)
 485           {
 486             fhEPrimMC  [mcPPhoton]->Fill(photonE ) ;
 487             fhPtPrimMC [mcPPhoton]->Fill(photonPt) ;
 488             fhPhiPrimMC[mcPPhoton]->Fill(photonE , photonPhi) ;
 489             fhYPrimMC[mcPPhoton]->Fill(photonE , photonEta) ;
 490           }
 491           if(inacceptance){
 492             fhEPrimMCAcc[mcPPhoton]  ->Fill(photonE ) ;
 493             fhPtPrimMCAcc[mcPPhoton] ->Fill(photonPt) ;
 494             fhPhiPrimMCAcc[mcPPhoton]->Fill(photonE , photonPhi) ;
 495             fhYPrimMCAcc[mcPPhoton]  ->Fill(photonE , photonY) ;
 496           }//Accepted
 497
 498
 499           //Origin of photon
 500           if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt) && fhEPrimMC[mcPPrompt])
 501           {
 502             fhYPrimMC[mcPPrompt]->Fill(photonPt, photonY) ;
 503             if(TMath::Abs(photonY) < 1.0){
 504               fhEPrimMC  [mcPPrompt]->Fill(photonE ) ;
 505               fhPtPrimMC [mcPPrompt]->Fill(photonPt) ;
 506               fhPhiPrimMC[mcPPrompt]->Fill(photonE , photonPhi) ;
 507               fhYPrimMC[mcPPrompt]->Fill(photonE , photonEta) ;
 508             }
 509             if(inacceptance){
 510               fhEPrimMCAcc[mcPPrompt]  ->Fill(photonE ) ;
 511               fhPtPrimMCAcc[mcPPrompt] ->Fill(photonPt) ;
 512               fhPhiPrimMCAcc[mcPPrompt]->Fill(photonE , photonPhi) ;
 513               fhYPrimMCAcc[mcPPrompt]  ->Fill(photonE , photonY) ;
 514             }//Accepted
 515           }
 516           else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation) && fhEPrimMC[mcPFragmentation] )
 517           {
 518             fhYPrimMC[mcPFragmentation]->Fill(photonPt, photonY) ;
 519             if(TMath::Abs(photonY) < 1.0){
 520               fhEPrimMC  [mcPFragmentation]->Fill(photonE ) ;
 521               fhPtPrimMC [mcPFragmentation]->Fill(photonPt) ;
 522               fhPhiPrimMC[mcPFragmentation]->Fill(photonE , photonPhi) ;
 523               fhYPrimMC[mcPFragmentation]->Fill(photonE , photonEta) ;
 524             }
 525             if(inacceptance){
 526               fhEPrimMCAcc[mcPFragmentation]  ->Fill(photonE ) ;
 527               fhPtPrimMCAcc[mcPFragmentation] ->Fill(photonPt) ;
 528               fhPhiPrimMCAcc[mcPFragmentation]->Fill(photonE , photonPhi) ;
 529               fhYPrimMCAcc[mcPFragmentation]  ->Fill(photonE , photonY) ;
 530             }//Accepted
 531           }
 532           else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCISR) && fhEPrimMC[mcPISR])
 533           {
 534             fhYPrimMC[mcPISR]->Fill(photonPt, photonY) ;
 535             if(TMath::Abs(photonY) < 1.0){
 536               fhEPrimMC  [mcPISR]->Fill(photonE ) ;
 537               fhPtPrimMC [mcPISR]->Fill(photonPt) ;
 538               fhPhiPrimMC[mcPISR]->Fill(photonE , photonPhi) ;
 539               fhYPrimMC[mcPISR]->Fill(photonE , photonEta) ;
 540             }
 541             if(inacceptance){
 542               fhEPrimMCAcc[mcPISR]  ->Fill(photonE ) ;
 543               fhPtPrimMCAcc[mcPISR] ->Fill(photonPt) ;
 544               fhPhiPrimMCAcc[mcPISR]->Fill(photonE , photonPhi) ;
 545               fhYPrimMCAcc[mcPISR]  ->Fill(photonE , photonY) ;
 546             }//Accepted
 547           }
 548           else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay)&& fhEPrimMC[mcPPi0Decay])
 549           {
 550             fhYPrimMC[mcPPi0Decay]->Fill(photonPt, photonY) ;
 551             if(TMath::Abs(photonY) < 1.0){
 552               fhEPrimMC  [mcPPi0Decay]->Fill(photonE ) ;
 553               fhPtPrimMC [mcPPi0Decay]->Fill(photonPt) ;
 554               fhPhiPrimMC[mcPPi0Decay]->Fill(photonE , photonPhi) ;
 555               fhYPrimMC[mcPPi0Decay]->Fill(photonE , photonEta) ;
 556             }
 557             if(inacceptance){
 558               fhEPrimMCAcc[mcPPi0Decay]  ->Fill(photonE ) ;
 559               fhPtPrimMCAcc[mcPPi0Decay] ->Fill(photonPt) ;
 560               fhPhiPrimMCAcc[mcPPi0Decay]->Fill(photonE , photonPhi) ;
 561               fhYPrimMCAcc[mcPPi0Decay]  ->Fill(photonE , photonY) ;
 562             }//Accepted
 563           }
 564           else if((GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay) ||
 565                    GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay) ) && fhEPrimMC[mcPOtherDecay])
 566           {
 567             fhYPrimMC[mcPOtherDecay]->Fill(photonPt, photonY) ;
 568             if(TMath::Abs(photonY) < 1.0){
 569               fhEPrimMC  [mcPOtherDecay]->Fill(photonE ) ;
 570               fhPtPrimMC [mcPOtherDecay]->Fill(photonPt) ;
 571               fhPhiPrimMC[mcPOtherDecay]->Fill(photonE , photonPhi) ;
 572               fhYPrimMC[mcPOtherDecay]->Fill(photonE , photonEta) ;
 573             }
 574             if(inacceptance){
 575               fhEPrimMCAcc[mcPOtherDecay]  ->Fill(photonE ) ;
 576               fhPtPrimMCAcc[mcPOtherDecay] ->Fill(photonPt) ;
 577               fhPhiPrimMCAcc[mcPOtherDecay]->Fill(photonE , photonPhi) ;
 578               fhYPrimMCAcc[mcPOtherDecay]  ->Fill(photonE , photonY) ;
 579             }//Accepted
 580           }
 581           else if(fhEPrimMC[mcPOther])
 582           {
 583             fhYPrimMC[mcPOther]->Fill(photonPt, photonY) ;
 584             if(TMath::Abs(photonY) < 1.0){
 585               fhEPrimMC  [mcPOther]->Fill(photonE ) ;
 586               fhPtPrimMC [mcPOther]->Fill(photonPt) ;
 587               fhPhiPrimMC[mcPOther]->Fill(photonE , photonPhi) ;
 588               fhYPrimMC[mcPOther]->Fill(photonE , photonEta) ;
 589             }
 590             if(inacceptance){
 591               fhEPrimMCAcc[mcPOther]  ->Fill(photonE ) ;
 592               fhPtPrimMCAcc[mcPOther] ->Fill(photonPt) ;
 593               fhPhiPrimMCAcc[mcPOther]->Fill(photonE , photonPhi) ;
 594               fhYPrimMCAcc[mcPOther]  ->Fill(photonE , photonY) ;
 595             }//Accepted
 596           }//Other origin
 597         }// Primary photon
 598       }//loop on primaries
 599
 600     }//mc array exists and data is MC
 601   }     // read AOD MC
 602 }
 603
 604 //__________________________________________________________________
 605 void  AliAnaPhoton::FillShowerShapeHistograms(AliVCluster* cluster, const Int_t mcTag){
 606
 607   //Fill cluster Shower Shape histograms
 608
 609   if(!fFillSSHistograms || GetMixedEvent()) return;
 610
 611   Float_t energy  = cluster->E();
 612   Int_t   ncells  = cluster->GetNCells();
 613   Float_t lambda0 = cluster->GetM02();
 614   Float_t lambda1 = cluster->GetM20();
 615   Float_t disp    = cluster->GetDispersion()*cluster->GetDispersion();
 616
 617   TLorentzVector mom;
 618   if(GetReader()->GetDataType() != AliCaloTrackReader::kMC){
 619     cluster->GetMomentum(mom,GetVertex(0)) ;}//Assume that come from vertex in straight line
 620   else{
 621     Double_t vertex[]={0,0,0};
 622     cluster->GetMomentum(mom,vertex) ;
 623   }
 624
 625   Float_t eta = mom.Eta();
 626   Float_t phi = mom.Phi();
 627   if(phi < 0) phi+=TMath::TwoPi();
 628
 629   fhLam0E ->Fill(energy,lambda0);
 630   fhLam1E ->Fill(energy,lambda1);
 631   fhDispE ->Fill(energy,disp);
 632
 633   if(fCalorimeter == "EMCAL" && GetModuleNumber(cluster) > 5){
 634     fhLam0ETRD->Fill(energy,lambda0);
 635     fhLam1ETRD->Fill(energy,lambda1);
 636     fhDispETRD->Fill(energy,disp);
 637   }
 638
 639   if(energy < 2){
 640     fhNCellsLam0LowE ->Fill(ncells,lambda0);
 641     fhNCellsLam1LowE ->Fill(ncells,lambda1);
 642     fhNCellsDispLowE ->Fill(ncells,disp);
 643
 644     fhLam1Lam0LowE  ->Fill(lambda1,lambda0);
 645     fhLam0DispLowE  ->Fill(lambda0,disp);
 646     fhDispLam1LowE  ->Fill(disp,lambda1);
 647     fhEtaLam0LowE   ->Fill(eta,lambda0);
 648     fhPhiLam0LowE   ->Fill(phi,lambda0);
 649
 650   }
 651   else {
 652     fhNCellsLam0HighE ->Fill(ncells,lambda0);
 653     fhNCellsLam1HighE ->Fill(ncells,lambda1);
 654     fhNCellsDispHighE ->Fill(ncells,disp);
 655
 656     fhLam1Lam0HighE  ->Fill(lambda1,lambda0);
 657     fhLam0DispHighE  ->Fill(lambda0,disp);
 658     fhDispLam1HighE  ->Fill(disp,lambda1);
 659     fhEtaLam0HighE   ->Fill(eta, lambda0);
 660     fhPhiLam0HighE   ->Fill(phi, lambda0);
 661   }
 662
 663   if(IsDataMC()){
 664
 665     AliVCaloCells* cells = 0;
 666     if(fCalorimeter == "EMCAL") cells = GetEMCALCells();
 667     else                        cells = GetPHOSCells();
 668
 669     //Fill histograms to check shape of embedded clusters
 670     Float_t fraction = 0;
 671     if(GetReader()->IsEmbeddedClusterSelectionOn()){//Only working for EMCAL
 672
 673       Float_t clusterE = 0; // recalculate in case corrections applied.
 674       Float_t cellE    = 0;
 675       for(Int_t icell = 0; icell < cluster->GetNCells(); icell++){
 676         cellE    = cells->GetCellAmplitude(cluster->GetCellAbsId(icell));
 677         clusterE+=cellE;
 678         fraction+=cellE*cluster->GetCellAmplitudeFraction(icell);
 679       }
 680
 681       //Fraction of total energy due to the embedded signal
 682       fraction/=clusterE;
 683
 684       if(GetDebug() > 1 ) printf("AliAnaPhoton::FillShowerShapeHistogram() - Energy fraction of embedded signal %2.3f, Energy %2.3f\n",fraction, clusterE);
 685
 686       fhEmbeddedSignalFractionEnergy->Fill(clusterE,fraction);
 687
 688     }  // embedded fraction
 689
 690     // Get the fraction of the cluster energy that carries the cell with highest energy
 691     Int_t absID             =-1 ;
 692     Float_t maxCellFraction = 0.;
 693
 694     absID = GetCaloUtils()->GetMaxEnergyCell(cells, cluster,maxCellFraction);
 695
 696     // Check the origin and fill histograms
 697     if( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton) &&
 698        !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCConversion) &&
 699        !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0) &&
 700        !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEta)){
 701       fhMCELambda0[mcssPhoton]    ->Fill(energy, lambda0);
 702       fhMCELambda1[mcssPhoton]    ->Fill(energy, lambda1);
 703       fhMCEDispersion[mcssPhoton] ->Fill(energy, disp);
 704       fhMCNCellsE[mcssPhoton]     ->Fill(energy, ncells);
 705       fhMCMaxCellDiffClusterE[mcssPhoton]->Fill(energy,maxCellFraction);
 706
 707       if     (energy < 2.){
 708         fhMCLambda0vsClusterMaxCellDiffE0[mcssPhoton]->Fill(lambda0, maxCellFraction);
 709         fhMCNCellsvsClusterMaxCellDiffE0[mcssPhoton] ->Fill(ncells,  maxCellFraction);
 710       }
 711       else if(energy < 6.){
 712         fhMCLambda0vsClusterMaxCellDiffE2[mcssPhoton]->Fill(lambda0, maxCellFraction);
 713         fhMCNCellsvsClusterMaxCellDiffE2[mcssPhoton] ->Fill(ncells,  maxCellFraction);
 714       }
 715       else{
 716         fhMCLambda0vsClusterMaxCellDiffE6[mcssPhoton]->Fill(lambda0, maxCellFraction);
 717         fhMCNCellsvsClusterMaxCellDiffE6[mcssPhoton] ->Fill(ncells,  maxCellFraction);
 718       }
 719
 720       if(!GetReader()->IsEmbeddedClusterSelectionOn()){
 721         //Check particle overlaps in cluster
 722
 723         //Compare the primary depositing more energy with the rest, if no photon/electron as comon ancestor (conversions), count as other particle
 724         Int_t ancPDG = 0, ancStatus = -1;
 725         TLorentzVector momentum; TVector3 prodVertex;
 726         Int_t ancLabel = 0;
 727         Int_t noverlaps = 1;
 728         for (UInt_t ilab = 0; ilab < cluster->GetNLabels(); ilab++ ) {
 729           ancLabel = GetMCAnalysisUtils()->CheckCommonAncestor(cluster->GetLabels()[0],cluster->GetLabels()[ilab], GetReader(),ancPDG,ancStatus,momentum,prodVertex);
 730           if(ancPDG!=22 && TMath::Abs(ancPDG)!=11) noverlaps++;
 731         }
 732
 733         if(noverlaps == 1){
 734           fhMCPhotonELambda0NoOverlap  ->Fill(energy, lambda0);
 735         }
 736         else if(noverlaps == 2){
 737           fhMCPhotonELambda0TwoOverlap ->Fill(energy, lambda0);
 738         }
 739         else if(noverlaps > 2){
 740           fhMCPhotonELambda0NOverlap   ->Fill(energy, lambda0);
 741         }
 742         else {
 743           printf("AliAnaPhoton::FillShowerShapeHistogram() - n overlaps = %d!!", noverlaps);
 744         }
 745       }//No embedding
 746
 747       //Fill histograms to check shape of embedded clusters
 748       if(GetReader()->IsEmbeddedClusterSelectionOn()){
 749
 750         if     (fraction > 0.9)
 751         {
 752           fhEmbedPhotonELambda0FullSignal   ->Fill(energy, lambda0);
 753         }
 754         else if(fraction > 0.5)
 755         {
 756           fhEmbedPhotonELambda0MostlySignal ->Fill(energy, lambda0);
 757         }
 758         else if(fraction > 0.1)
 759         {
 760           fhEmbedPhotonELambda0MostlyBkg    ->Fill(energy, lambda0);
 761         }
 762         else
 763         {
 764           fhEmbedPhotonELambda0FullBkg      ->Fill(energy, lambda0);
 765         }
 766       } // embedded
 767
 768     }//photon   no conversion
 769     else if  ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCElectron)){
 770       fhMCELambda0[mcssElectron]    ->Fill(energy, lambda0);
 771       fhMCELambda1[mcssElectron]    ->Fill(energy, lambda1);
 772       fhMCEDispersion[mcssElectron] ->Fill(energy, disp);
 773       fhMCNCellsE[mcssElectron]     ->Fill(energy, ncells);
 774       fhMCMaxCellDiffClusterE[mcssElectron]->Fill(energy,maxCellFraction);
 775
 776       if     (energy < 2.){
 777         fhMCLambda0vsClusterMaxCellDiffE0[mcssElectron]->Fill(lambda0, maxCellFraction);
 778         fhMCNCellsvsClusterMaxCellDiffE0[mcssElectron] ->Fill(ncells,  maxCellFraction);
 779       }
 780       else if(energy < 6.){
 781         fhMCLambda0vsClusterMaxCellDiffE2[mcssElectron]->Fill(lambda0, maxCellFraction);
 782         fhMCNCellsvsClusterMaxCellDiffE2[mcssElectron] ->Fill(ncells,  maxCellFraction);
 783       }
 784       else{
 785         fhMCLambda0vsClusterMaxCellDiffE6[mcssElectron]->Fill(lambda0, maxCellFraction);
 786         fhMCNCellsvsClusterMaxCellDiffE6[mcssElectron] ->Fill(ncells,  maxCellFraction);
 787       }
 788     }//electron
 789     else if  ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton) &&
 790               GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCConversion) ){
 791       fhMCELambda0[mcssConversion]    ->Fill(energy, lambda0);
 792       fhMCELambda1[mcssConversion]    ->Fill(energy, lambda1);
 793       fhMCEDispersion[mcssConversion] ->Fill(energy, disp);
 794       fhMCNCellsE[mcssConversion]     ->Fill(energy, ncells);
 795       fhMCMaxCellDiffClusterE[mcssConversion]->Fill(energy,maxCellFraction);
 796
 797       if     (energy < 2.){
 798         fhMCLambda0vsClusterMaxCellDiffE0[mcssConversion]->Fill(lambda0, maxCellFraction);
 799         fhMCNCellsvsClusterMaxCellDiffE0[mcssConversion] ->Fill(ncells,  maxCellFraction);
 800       }
 801       else if(energy < 6.){
 802         fhMCLambda0vsClusterMaxCellDiffE2[mcssConversion]->Fill(lambda0, maxCellFraction);
 803         fhMCNCellsvsClusterMaxCellDiffE2[mcssConversion] ->Fill(ncells,  maxCellFraction);
 804       }
 805       else{
 806         fhMCLambda0vsClusterMaxCellDiffE6[mcssConversion]->Fill(lambda0, maxCellFraction);
 807         fhMCNCellsvsClusterMaxCellDiffE6[mcssConversion] ->Fill(ncells,  maxCellFraction);
 808       }
 809
 810     }//conversion photon
 811     else if  ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0)  ){
 812       fhMCELambda0[mcssPi0]    ->Fill(energy, lambda0);
 813       fhMCELambda1[mcssPi0]    ->Fill(energy, lambda1);
 814       fhMCEDispersion[mcssPi0] ->Fill(energy, disp);
 815       fhMCNCellsE[mcssPi0]     ->Fill(energy, ncells);
 816       fhMCMaxCellDiffClusterE[mcssPi0]->Fill(energy,maxCellFraction);
 817
 818       if     (energy < 2.){
 819         fhMCLambda0vsClusterMaxCellDiffE0[mcssPi0]->Fill(lambda0, maxCellFraction);
 820         fhMCNCellsvsClusterMaxCellDiffE0[mcssPi0] ->Fill(ncells,  maxCellFraction);
 821       }
 822       else if(energy < 6.){
 823         fhMCLambda0vsClusterMaxCellDiffE2[mcssPi0]->Fill(lambda0, maxCellFraction);
 824         fhMCNCellsvsClusterMaxCellDiffE2[mcssPi0] ->Fill(ncells,  maxCellFraction);
 825       }
 826       else{
 827         fhMCLambda0vsClusterMaxCellDiffE6[mcssPi0]->Fill(lambda0, maxCellFraction);
 828         fhMCNCellsvsClusterMaxCellDiffE6[mcssPi0] ->Fill(ncells,  maxCellFraction);
 829       }
 830
 831       //Fill histograms to check shape of embedded clusters
 832       if(GetReader()->IsEmbeddedClusterSelectionOn()){
 833
 834         if     (fraction > 0.9)
 835         {
 836           fhEmbedPi0ELambda0FullSignal   ->Fill(energy, lambda0);
 837         }
 838         else if(fraction > 0.5)
 839         {
 840           fhEmbedPi0ELambda0MostlySignal ->Fill(energy, lambda0);
 841         }
 842         else if(fraction > 0.1)
 843         {
 844           fhEmbedPi0ELambda0MostlyBkg    ->Fill(energy, lambda0);
 845         }
 846         else
 847         {
 848           fhEmbedPi0ELambda0FullBkg      ->Fill(energy, lambda0);
 849         }
 850       } // embedded
 851
 852     }//pi0
 853     else if  ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEta)  ){
 854       fhMCELambda0[mcssEta]    ->Fill(energy, lambda0);
 855       fhMCELambda1[mcssEta]    ->Fill(energy, lambda1);
 856       fhMCEDispersion[mcssEta] ->Fill(energy, disp);
 857       fhMCNCellsE[mcssEta]     ->Fill(energy, ncells);
 858       fhMCMaxCellDiffClusterE[mcssEta]->Fill(energy,maxCellFraction);
 859
 860       if     (energy < 2.){
 861         fhMCLambda0vsClusterMaxCellDiffE0[mcssEta]->Fill(lambda0, maxCellFraction);
 862         fhMCNCellsvsClusterMaxCellDiffE0[mcssEta] ->Fill(ncells,  maxCellFraction);
 863       }
 864       else if(energy < 6.){
 865         fhMCLambda0vsClusterMaxCellDiffE2[mcssEta]->Fill(lambda0, maxCellFraction);
 866         fhMCNCellsvsClusterMaxCellDiffE2[mcssEta] ->Fill(ncells,  maxCellFraction);
 867       }
 868       else{
 869         fhMCLambda0vsClusterMaxCellDiffE6[mcssEta]->Fill(lambda0, maxCellFraction);
 870         fhMCNCellsvsClusterMaxCellDiffE6[mcssEta] ->Fill(ncells,  maxCellFraction);
 871       }
 872
 873     }//eta
 874     else {
 875       fhMCELambda0[mcssOther]    ->Fill(energy, lambda0);
 876       fhMCELambda1[mcssOther]    ->Fill(energy, lambda1);
 877       fhMCEDispersion[mcssOther] ->Fill(energy, disp);
 878       fhMCNCellsE[mcssOther]     ->Fill(energy, ncells);
 879       fhMCMaxCellDiffClusterE[mcssOther]->Fill(energy,maxCellFraction);
 880
 881       if     (energy < 2.){
 882         fhMCLambda0vsClusterMaxCellDiffE0[mcssOther]->Fill(lambda0, maxCellFraction);
 883         fhMCNCellsvsClusterMaxCellDiffE0[mcssOther] ->Fill(ncells,  maxCellFraction);
 884       }
 885       else if(energy < 6.){
 886         fhMCLambda0vsClusterMaxCellDiffE2[mcssOther]->Fill(lambda0, maxCellFraction);
 887         fhMCNCellsvsClusterMaxCellDiffE2[mcssOther] ->Fill(ncells,  maxCellFraction);
 888       }
 889       else{
 890         fhMCLambda0vsClusterMaxCellDiffE6[mcssOther]->Fill(lambda0, maxCellFraction);
 891         fhMCNCellsvsClusterMaxCellDiffE6[mcssOther] ->Fill(ncells,  maxCellFraction);
 892       }
 893
 894     }//other particles
 895
 896   }//MC data
 897
 898 }
 899
 900 //________________________________________________________________________
 901 TObjString *  AliAnaPhoton::GetAnalysisCuts()
 902 {
 903   //Save parameters used for analysis
 904   TString parList ; //this will be list of parameters used for this analysis.
 905   const Int_t buffersize = 255;
 906   char onePar[buffersize] ;
 907
 908   snprintf(onePar,buffersize,"--- AliAnaPhoton ---\n") ;
 909   parList+=onePar ;
 910   snprintf(onePar,buffersize,"Calorimeter: %s\n",fCalorimeter.Data()) ;
 911   parList+=onePar ;
 912   snprintf(onePar,buffersize,"fMinDist =%2.2f (Minimal distance to bad channel to accept cluster) \n",fMinDist) ;
 913   parList+=onePar ;
 914   snprintf(onePar,buffersize,"fMinDist2=%2.2f (Cuts on Minimal distance to study acceptance evaluation) \n",fMinDist2) ;
 915   parList+=onePar ;
 916   snprintf(onePar,buffersize,"fMinDist3=%2.2f (One more cut on distance used for acceptance-efficiency study) \n",fMinDist3) ;
 917   parList+=onePar ;
 918   snprintf(onePar,buffersize,"fRejectTrackMatch: %d\n",fRejectTrackMatch) ;
 919   parList+=onePar ;
 920
 921   //Get parameters set in base class.
 922   parList += GetBaseParametersList() ;
 923
 924   //Get parameters set in PID class.
 925   parList += GetCaloPID()->GetPIDParametersList() ;
 926
 927   //Get parameters set in FiducialCut class (not available yet)
 928   //parlist += GetFidCut()->GetFidCutParametersList()
 929
 930   return new TObjString(parList) ;
 931 }
 932
 933 //________________________________________________________________________
 934 TList *  AliAnaPhoton::GetCreateOutputObjects()
 935 {
 936   // Create histograms to be saved in output file and
 937   // store them in outputContainer
 938   TList * outputContainer = new TList() ;
 939   outputContainer->SetName("PhotonHistos") ;
 940
 941   Int_t nptbins  = GetHistoPtBins();  Float_t ptmax  = GetHistoPtMax();  Float_t ptmin  = GetHistoPtMin();
 942   Int_t nphibins = GetHistoPhiBins(); Float_t phimax = GetHistoPhiMax(); Float_t phimin = GetHistoPhiMin();
 943   Int_t netabins = GetHistoEtaBins(); Float_t etamax = GetHistoEtaMax(); Float_t etamin = GetHistoEtaMin();
 944   Int_t ssbins   = GetHistoShowerShapeBins();  Float_t ssmax   = GetHistoShowerShapeMax();  Float_t ssmin   = GetHistoShowerShapeMin();
 945   Int_t nbins    = GetHistoNClusterCellBins(); Int_t   nmax    = GetHistoNClusterCellMax(); Int_t   nmin    = GetHistoNClusterCellMin();
 946   Int_t ntimebins= GetHistoTimeBins();         Float_t timemax = GetHistoTimeMax();         Float_t timemin = GetHistoTimeMin();
 947
 948   fhNCellsE  = new TH2F ("hNCellsE","# of cells in cluster vs E of clusters", nptbins,ptmin,ptmax, nbins,nmin,nmax);
 949   fhNCellsE->SetXTitle("E (GeV)");
 950   fhNCellsE->SetYTitle("# of cells in cluster");
 951   outputContainer->Add(fhNCellsE);
 952
 953   fhTimeE  = new TH2F ("hTimeE","time of cluster vs E of clusters", nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
 954   fhTimeE->SetXTitle("E (GeV)");
 955   fhTimeE->SetYTitle("time (ns)");
 956   outputContainer->Add(fhTimeE);
 957
 958   fhMaxCellDiffClusterE  = new TH2F ("hMaxCellDiffClusterE","energy vs difference of cluster energy - max cell energy / cluster energy, good clusters",
 959                                     nptbins,ptmin,ptmax, 500,0,1.);
 960   fhMaxCellDiffClusterE->SetXTitle("E_{cluster} (GeV) ");
 961   fhMaxCellDiffClusterE->SetYTitle("(E_{cluster} - E_{cell max})/ E_{cluster}");
 962   outputContainer->Add(fhMaxCellDiffClusterE);
 963
 964   fhEPhoton  = new TH1F("hEPhoton","Number of #gamma over calorimeter vs energy",nptbins,ptmin,ptmax);
 965   fhEPhoton->SetYTitle("N");
 966   fhEPhoton->SetXTitle("E_{#gamma}(GeV)");
 967   outputContainer->Add(fhEPhoton) ;
 968
 969   fhPtPhoton  = new TH1F("hPtPhoton","Number of #gamma over calorimeter vs p_{T}",nptbins,ptmin,ptmax);
 970   fhPtPhoton->SetYTitle("N");
 971   fhPtPhoton->SetXTitle("p_{T #gamma}(GeV/c)");
 972   outputContainer->Add(fhPtPhoton) ;
 973
 974   fhPhiPhoton  = new TH2F
 975     ("hPhiPhoton","#phi_{#gamma} vs p_{T}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 976   fhPhiPhoton->SetYTitle("#phi (rad)");
 977   fhPhiPhoton->SetXTitle("p_{T #gamma} (GeV/c)");
 978   outputContainer->Add(fhPhiPhoton) ;
 979
 980   fhEtaPhoton  = new TH2F
 981     ("hEtaPhoton","#eta_{#gamma} vs p_{T}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 982   fhEtaPhoton->SetYTitle("#eta");
 983   fhEtaPhoton->SetXTitle("p_{T #gamma} (GeV/c)");
 984   outputContainer->Add(fhEtaPhoton) ;
 985
 986   fhEtaPhiPhoton  = new TH2F
 987   ("hEtaPhiPhoton","#eta vs #phi",netabins,etamin,etamax,nphibins,phimin,phimax);
 988   fhEtaPhiPhoton->SetYTitle("#phi (rad)");
 989   fhEtaPhiPhoton->SetXTitle("#eta");
 990   outputContainer->Add(fhEtaPhiPhoton) ;
 991   if(GetMinPt() < 0.5){
 992     fhEtaPhi05Photon  = new TH2F
 993     ("hEtaPhi05Photon","#eta vs #phi, E > 0.5",netabins,etamin,etamax,nphibins,phimin,phimax);
 994     fhEtaPhi05Photon->SetYTitle("#phi (rad)");
 995     fhEtaPhi05Photon->SetXTitle("#eta");
 996     outputContainer->Add(fhEtaPhi05Photon) ;
 997   }
 998
 999   //Shower shape
1000   if(fFillSSHistograms){
1001
1002     fhLam0E  = new TH2F ("hLam0E","#lambda_{0}^{2} vs E", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1003     fhLam0E->SetYTitle("#lambda_{0}^{2}");
1004     fhLam0E->SetXTitle("E (GeV)");
1005     outputContainer->Add(fhLam0E);
1006
1007     fhLam1E  = new TH2F ("hLam1E","#lambda_{1}^{2} vs E", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1008     fhLam1E->SetYTitle("#lambda_{1}^{2}");
1009     fhLam1E->SetXTitle("E (GeV)");
1010     outputContainer->Add(fhLam1E);
1011
1012     fhDispE  = new TH2F ("hDispE"," dispersion^{2} vs E", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1013     fhDispE->SetYTitle("D^{2}");
1014     fhDispE->SetXTitle("E (GeV) ");
1015     outputContainer->Add(fhDispE);
1016
1017     if(fCalorimeter == "EMCAL"){
1018       fhLam0ETRD  = new TH2F ("hLam0ETRD","#lambda_{0}^{2} vs E, EMCAL SM covered by TRD", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1019       fhLam0ETRD->SetYTitle("#lambda_{0}^{2}");
1020       fhLam0ETRD->SetXTitle("E (GeV)");
1021       outputContainer->Add(fhLam0ETRD);
1022
1023       fhLam1ETRD  = new TH2F ("hLam1ETRD","#lambda_{1}^{2} vs E, EMCAL SM covered by TRD", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1024       fhLam1ETRD->SetYTitle("#lambda_{1}^{2}");
1025       fhLam1ETRD->SetXTitle("E (GeV)");
1026       outputContainer->Add(fhLam1ETRD);
1027
1028       fhDispETRD  = new TH2F ("hDispETRD"," dispersion^{2} vs E, EMCAL SM covered by TRD", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1029       fhDispETRD->SetYTitle("Dispersion^{2}");
1030       fhDispETRD->SetXTitle("E (GeV) ");
1031       outputContainer->Add(fhDispETRD);
1032     }
1033
1034     fhNCellsLam0LowE  = new TH2F ("hNCellsLam0LowE","N_{cells} in cluster vs #lambda_{0}^{2}, E < 2 GeV", nbins,nmin, nmax, ssbins,ssmin,ssmax);
1035     fhNCellsLam0LowE->SetXTitle("N_{cells}");
1036     fhNCellsLam0LowE->SetYTitle("#lambda_{0}^{2}");
1037     outputContainer->Add(fhNCellsLam0LowE);
1038
1039     fhNCellsLam0HighE  = new TH2F ("hNCellsLam0HighE","N_{cells} in cluster vs #lambda_{0}^{2}, E > 2 GeV", nbins,nmin, nmax, ssbins,ssmin,ssmax);
1040     fhNCellsLam0HighE->SetXTitle("N_{cells}");
1041     fhNCellsLam0HighE->SetYTitle("#lambda_{0}^{2}");
1042     outputContainer->Add(fhNCellsLam0HighE);
1043
1044     fhNCellsLam1LowE  = new TH2F ("hNCellsLam1LowE","N_{cells} in cluster vs #lambda_{1}^{2}, E < 2 GeV", nbins,nmin, nmax, ssbins,ssmin,ssmax);
1045     fhNCellsLam1LowE->SetXTitle("N_{cells}");
1046     fhNCellsLam1LowE->SetYTitle("#lambda_{0}^{2}");
1047     outputContainer->Add(fhNCellsLam1LowE);
1048
1049     fhNCellsLam1HighE  = new TH2F ("hNCellsLam1HighE","N_{cells} in cluster vs #lambda_{1}^{2}, E > 2 GeV", nbins,nmin, nmax, ssbins,ssmin,ssmax);
1050     fhNCellsLam1HighE->SetXTitle("N_{cells}");
1051     fhNCellsLam1HighE->SetYTitle("#lambda_{0}^{2}");
1052     outputContainer->Add(fhNCellsLam1HighE);
1053
1054     fhNCellsDispLowE  = new TH2F ("hNCellsDispLowE","N_{cells} in cluster vs dispersion^{2}, E < 2 GeV", nbins,nmin, nmax, ssbins,ssmin,ssmax);
1055     fhNCellsDispLowE->SetXTitle("N_{cells}");
1056     fhNCellsDispLowE->SetYTitle("D^{2}");
1057     outputContainer->Add(fhNCellsDispLowE);
1058
1059     fhNCellsDispHighE  = new TH2F ("hNCellsDispHighE","N_{cells} in cluster vs dispersion^{2}, E < 2 GeV", nbins,nmin, nmax, ssbins,ssmin,ssmax);
1060     fhNCellsDispHighE->SetXTitle("N_{cells}");
1061     fhNCellsDispHighE->SetYTitle("D^{2}");
1062     outputContainer->Add(fhNCellsDispHighE);
1063
1064     fhEtaLam0LowE  = new TH2F ("hEtaLam0LowE","#eta vs #lambda_{0}^{2}, E < 2 GeV", netabins,etamin,etamax, ssbins,ssmin,ssmax);
1065     fhEtaLam0LowE->SetYTitle("#lambda_{0}^{2}");
1066     fhEtaLam0LowE->SetXTitle("#eta");
1067     outputContainer->Add(fhEtaLam0LowE);
1068
1069     fhPhiLam0LowE  = new TH2F ("hPhiLam0LowE","#phi vs #lambda_{0}^{2}, E < 2 GeV", nphibins,phimin,phimax, ssbins,ssmin,ssmax);
1070     fhPhiLam0LowE->SetYTitle("#lambda_{0}^{2}");
1071     fhPhiLam0LowE->SetXTitle("#phi");
1072     outputContainer->Add(fhPhiLam0LowE);
1073
1074     fhEtaLam0HighE  = new TH2F ("hEtaLam0HighE","#eta vs #lambda_{0}^{2}, E > 2 GeV", netabins,etamin,etamax, ssbins,ssmin,ssmax);
1075     fhEtaLam0HighE->SetYTitle("#lambda_{0}^{2}");
1076     fhEtaLam0HighE->SetXTitle("#eta");
1077     outputContainer->Add(fhEtaLam0HighE);
1078
1079     fhPhiLam0HighE  = new TH2F ("hPhiLam0HighE","#phi vs #lambda_{0}^{2}, E > 2 GeV", nphibins,phimin,phimax, ssbins,ssmin,ssmax);
1080     fhPhiLam0HighE->SetYTitle("#lambda_{0}^{2}");
1081     fhPhiLam0HighE->SetXTitle("#phi");
1082     outputContainer->Add(fhPhiLam0HighE);
1083
1084     fhLam1Lam0LowE  = new TH2F ("hLam1Lam0LowE","#lambda_{0}^{2} vs #lambda_{1}^{2} in cluster of E < 2 GeV",  ssbins,ssmin,ssmax, ssbins,ssmin,ssmax);
1085     fhLam1Lam0LowE->SetYTitle("#lambda_{0}^{2}");
1086     fhLam1Lam0LowE->SetXTitle("#lambda_{1}^{2}");
1087     outputContainer->Add(fhLam1Lam0LowE);
1088
1089     fhLam1Lam0HighE  = new TH2F ("hLam1Lam0HighE","#lambda_{0}^{2} vs #lambda_{1}^{2} in cluster of E > 2 GeV",  ssbins,ssmin,ssmax, ssbins,ssmin,ssmax);
1090     fhLam1Lam0HighE->SetYTitle("#lambda_{0}^{2}");
1091     fhLam1Lam0HighE->SetXTitle("#lambda_{1}^{2}");
1092     outputContainer->Add(fhLam1Lam0HighE);
1093
1094     fhLam0DispLowE  = new TH2F ("hLam0DispLowE","#lambda_{0}^{2} vs dispersion^{2} in cluster of E < 2 GeV",  ssbins,ssmin,ssmax, ssbins,ssmin,ssmax);
1095     fhLam0DispLowE->SetXTitle("#lambda_{0}^{2}");
1096     fhLam0DispLowE->SetYTitle("D^{2}");
1097     outputContainer->Add(fhLam0DispLowE);
1098
1099     fhLam0DispHighE  = new TH2F ("hLam0DispHighE","#lambda_{0}^{2} vs dispersion^{2} in cluster of E > 2 GeV",  ssbins,ssmin,ssmax, ssbins,ssmin,ssmax);
1100     fhLam0DispHighE->SetXTitle("#lambda_{0}^{2}");
1101     fhLam0DispHighE->SetYTitle("D^{2}");
1102     outputContainer->Add(fhLam0DispHighE);
1103
1104     fhDispLam1LowE  = new TH2F ("hDispLam1LowE","Dispersion^{2} vs #lambda_{1}^{2} in cluster of E < 2 GeV",  ssbins,ssmin,ssmax, ssbins,ssmin,ssmax);
1105     fhDispLam1LowE->SetXTitle("D^{2}");
1106     fhDispLam1LowE->SetYTitle("#lambda_{1}^{2}");
1107     outputContainer->Add(fhDispLam1LowE);
1108
1109     fhDispLam1HighE  = new TH2F ("hDispLam1HighE","Dispersion^{2} vs #lambda_{1^{2}} in cluster of E > 2 GeV",  ssbins,ssmin,ssmax, ssbins,ssmin,ssmax);
1110     fhDispLam1HighE->SetXTitle("D^{2}");
1111     fhDispLam1HighE->SetYTitle("#lambda_{1}^{2}");
1112     outputContainer->Add(fhDispLam1HighE);
1113
1114   } // Shower shape
1115
1116
1117   if(IsDataMC()){
1118
1119     TString ptype[] = { "#gamma", "#gamma_{#pi decay}","#gamma_{other decay}", "#pi^{0}","#eta",
1120                         "e^{#pm}","#gamma->e^{#pm}","hadron?","Anti-N","Anti-P",
1121                         "#gamma_{prompt}","#gamma_{fragmentation}","#gamma_{ISR}","String"                                } ;
1122
1123     TString pname[] = { "Photon","PhotonPi0Decay","PhotonOtherDecay","Pi0","Eta","Electron",
1124                         "Conversion", "Hadron", "AntiNeutron","AntiProton",
1125                         "PhotonPrompt","PhotonFragmentation","PhotonISR","String" } ;
1126
1127     for(Int_t i = 0; i < fNOriginHistograms; i++){
1128
1129       fhMCE[i]  = new TH1F(Form("hE_MC%s",pname[i].Data()),
1130                                 Form("cluster from %s : E ",ptype[i].Data()),
1131                                 nptbins,ptmin,ptmax);
1132       fhMCE[i]->SetXTitle("E (GeV)");
1133       outputContainer->Add(fhMCE[i]) ;
1134
1135       fhMCPt[i]  = new TH1F(Form("hPt_MC%s",pname[i].Data()),
1136                            Form("cluster from %s : p_{T} ",ptype[i].Data()),
1137                            nptbins,ptmin,ptmax);
1138       fhMCPt[i]->SetXTitle("p_{T} (GeV/c)");
1139       outputContainer->Add(fhMCPt[i]) ;
1140
1141       fhMCEta[i]  = new TH2F(Form("hEta_MC%s",pname[i].Data()),
1142                            Form("cluster from %s : #eta ",ptype[i].Data()),
1143                            nptbins,ptmin,ptmax,netabins,etamin,etamax);
1144       fhMCEta[i]->SetYTitle("#eta");
1145       fhMCEta[i]->SetXTitle("E (GeV)");
1146       outputContainer->Add(fhMCEta[i]) ;
1147
1148       fhMCPhi[i]  = new TH2F(Form("hPhi_MC%s",pname[i].Data()),
1149                            Form("cluster from %s : #phi ",ptype[i].Data()),
1150                            nptbins,ptmin,ptmax,nphibins,phimin,phimax);
1151       fhMCPhi[i]->SetYTitle("#phi (rad)");
1152       fhMCPhi[i]->SetXTitle("E (GeV)");
1153       outputContainer->Add(fhMCPhi[i]) ;
1154
1155
1156       fhMCDeltaE[i]  = new TH2F (Form("hDeltaE_MC%s",pname[i].Data()),
1157                                  Form("MC - Reco E from %s",pname[i].Data()),
1158                                  nptbins,ptmin,ptmax, 200,-50,50);
1159       fhMCDeltaE[i]->SetXTitle("#Delta E (GeV)");
1160       outputContainer->Add(fhMCDeltaE[i]);
1161
1162       fhMCDeltaPt[i]  = new TH2F (Form("hDeltaPt_MC%s",pname[i].Data()),
1163                                   Form("MC - Reco p_{T} from %s",pname[i].Data()),
1164                                   nptbins,ptmin,ptmax, 200,-50,50);
1165       fhMCDeltaPt[i]->SetXTitle("#Delta p_{T} (GeV/c)");
1166       outputContainer->Add(fhMCDeltaPt[i]);
1167
1168       fhMC2E[i]  = new TH2F (Form("h2E_MC%s",pname[i].Data()),
1169                              Form("E distribution, reconstructed vs generated from %s",pname[i].Data()),
1170                              nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
1171       fhMC2E[i]->SetXTitle("E_{rec} (GeV)");
1172       fhMC2E[i]->SetYTitle("E_{gen} (GeV)");
1173       outputContainer->Add(fhMC2E[i]);
1174
1175       fhMC2Pt[i]  = new TH2F (Form("h2Pt_MC%s",pname[i].Data()),
1176                               Form("p_T distribution, reconstructed vs generated from %s",pname[i].Data()),
1177                               nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
1178       fhMC2Pt[i]->SetXTitle("p_{T,rec} (GeV/c)");
1179       fhMC2Pt[i]->SetYTitle("p_{T,gen} (GeV/c)");
1180       outputContainer->Add(fhMC2Pt[i]);
1181
1182
1183     }
1184
1185     TString pptype[] = { "#gamma", "#gamma_{#pi decay}","#gamma_{other decay}","hadron?",
1186                          "#gamma_{prompt}","#gamma_{fragmentation}","#gamma_{ISR}"} ;
1187
1188     TString ppname[] = { "Photon","PhotonPi0Decay","PhotonOtherDecay","Hadron",
1189                          "PhotonPrompt","PhotonFragmentation","PhotonISR"} ;
1190
1191     for(Int_t i = 0; i < fNPrimaryHistograms; i++){
1192       fhEPrimMC[i]  = new TH1F(Form("hEPrim_MC%s",ppname[i].Data()),
1193                            Form("primary photon %s : E ",pptype[i].Data()),
1194                            nptbins,ptmin,ptmax);
1195       fhEPrimMC[i]->SetXTitle("E (GeV)");
1196       outputContainer->Add(fhEPrimMC[i]) ;
1197
1198       fhPtPrimMC[i]  = new TH1F(Form("hPtPrim_MC%s",ppname[i].Data()),
1199                             Form("primary photon %s : p_{T} ",pptype[i].Data()),
1200                             nptbins,ptmin,ptmax);
1201       fhPtPrimMC[i]->SetXTitle("p_{T} (GeV/c)");
1202       outputContainer->Add(fhPtPrimMC[i]) ;
1203
1204       fhYPrimMC[i]  = new TH2F(Form("hYPrim_MC%s",ppname[i].Data()),
1205                              Form("primary photon %s : Rapidity ",pptype[i].Data()),
1206                              nptbins,ptmin,ptmax,800,-8,8);
1207       fhYPrimMC[i]->SetYTitle("Rapidity");
1208       fhYPrimMC[i]->SetXTitle("E (GeV)");
1209       outputContainer->Add(fhYPrimMC[i]) ;
1210
1211       fhPhiPrimMC[i]  = new TH2F(Form("hPhiPrim_MC%s",ppname[i].Data()),
1212                              Form("primary photon %s : #phi ",pptype[i].Data()),
1213                              nptbins,ptmin,ptmax,nphibins,phimin,phimax);
1214       fhPhiPrimMC[i]->SetYTitle("#phi (rad)");
1215       fhPhiPrimMC[i]->SetXTitle("E (GeV)");
1216       outputContainer->Add(fhPhiPrimMC[i]) ;
1217
1218
1219       fhEPrimMCAcc[i]  = new TH1F(Form("hEPrimAcc_MC%s",ppname[i].Data()),
1220                                Form("primary photon %s in acceptance: E ",pptype[i].Data()),
1221                                nptbins,ptmin,ptmax);
1222       fhEPrimMCAcc[i]->SetXTitle("E (GeV)");
1223       outputContainer->Add(fhEPrimMCAcc[i]) ;
1224
1225       fhPtPrimMCAcc[i]  = new TH1F(Form("hPtPrimAcc_MC%s",ppname[i].Data()),
1226                                 Form("primary photon %s in acceptance: p_{T} ",pptype[i].Data()),
1227                                 nptbins,ptmin,ptmax);
1228       fhPtPrimMCAcc[i]->SetXTitle("p_{T} (GeV/c)");
1229       outputContainer->Add(fhPtPrimMCAcc[i]) ;
1230
1231       fhYPrimMCAcc[i]  = new TH2F(Form("hYPrimAcc_MC%s",ppname[i].Data()),
1232                                  Form("primary photon %s in acceptance: Rapidity ",pptype[i].Data()),
1233                                  nptbins,ptmin,ptmax,100,-1,1);
1234       fhYPrimMCAcc[i]->SetYTitle("Rapidity");
1235       fhYPrimMCAcc[i]->SetXTitle("E (GeV)");
1236       outputContainer->Add(fhYPrimMCAcc[i]) ;
1237
1238       fhPhiPrimMCAcc[i]  = new TH2F(Form("hPhiPrimAcc_MC%s",ppname[i].Data()),
1239                                  Form("primary photon %s in acceptance: #phi ",pptype[i].Data()),
1240                                  nptbins,ptmin,ptmax,nphibins,phimin,phimax);
1241       fhPhiPrimMCAcc[i]->SetYTitle("#phi (rad)");
1242       fhPhiPrimMCAcc[i]->SetXTitle("E (GeV)");
1243       outputContainer->Add(fhPhiPrimMCAcc[i]) ;
1244
1245     }
1246
1247     if(fFillSSHistograms){
1248
1249       TString ptypess[] = { "#gamma","hadron?","#pi^{0}","#eta","#gamma->e^{#pm}","e^{#pm}"} ;
1250
1251       TString pnamess[] = { "Photon","Hadron","Pi0","Eta","Conversion","Electron"} ;
1252
1253       for(Int_t i = 0; i < 6; i++){
1254
1255         fhMCELambda0[i]  = new TH2F(Form("hELambda0_MC%s",pnamess[i].Data()),
1256                                     Form("cluster from %s : E vs #lambda_{0}^{2}",ptypess[i].Data()),
1257                                     nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1258         fhMCELambda0[i]->SetYTitle("#lambda_{0}^{2}");
1259         fhMCELambda0[i]->SetXTitle("E (GeV)");
1260         outputContainer->Add(fhMCELambda0[i]) ;
1261
1262         fhMCELambda1[i]  = new TH2F(Form("hELambda1_MC%s",pnamess[i].Data()),
1263                                     Form("cluster from %s : E vs #lambda_{1}^{2}",ptypess[i].Data()),
1264                                     nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1265         fhMCELambda1[i]->SetYTitle("#lambda_{1}^{2}");
1266         fhMCELambda1[i]->SetXTitle("E (GeV)");
1267         outputContainer->Add(fhMCELambda1[i]) ;
1268
1269         fhMCEDispersion[i]  = new TH2F(Form("hEDispersion_MC%s",pnamess[i].Data()),
1270                                        Form("cluster from %s : E vs dispersion^{2}",ptypess[i].Data()),
1271                                        nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1272         fhMCEDispersion[i]->SetYTitle("D^{2}");
1273         fhMCEDispersion[i]->SetXTitle("E (GeV)");
1274         outputContainer->Add(fhMCEDispersion[i]) ;
1275
1276         fhMCNCellsE[i]  = new TH2F (Form("hNCellsE_MC%s",pnamess[i].Data()),
1277                                Form("# of cells in cluster from %s vs E of clusters",ptypess[i].Data()),
1278                                     nptbins,ptmin,ptmax, nbins,nmin,nmax);
1279         fhMCNCellsE[i]->SetXTitle("E (GeV)");
1280         fhMCNCellsE[i]->SetYTitle("# of cells in cluster");
1281         outputContainer->Add(fhMCNCellsE[i]);
1282
1283         fhMCMaxCellDiffClusterE[i]  = new TH2F (Form("hMaxCellDiffClusterE_MC%s",pnamess[i].Data()),
1284                                              Form("energy vs difference of cluster energy from %s - max cell energy / cluster energy, good clusters",ptypess[i].Data()),
1285                                            nptbins,ptmin,ptmax, 500,0,1.);
1286         fhMCMaxCellDiffClusterE[i]->SetXTitle("E_{cluster} (GeV) ");
1287         fhMCMaxCellDiffClusterE[i]->SetYTitle("(E_{cluster} - E_{cell max})/ E_{cluster}");
1288         outputContainer->Add(fhMCMaxCellDiffClusterE[i]);
1289
1290         fhMCLambda0vsClusterMaxCellDiffE0[i]  = new TH2F(Form("hLambda0vsClusterMaxCellDiffE0_MC%s",pnamess[i].Data()),
1291                                     Form("cluster from %s : #lambda^{2}_{0} vs fraction of energy carried by max cell, E < 2 GeV",ptypess[i].Data()),
1292                                     ssbins,ssmin,ssmax,500,0,1.);
1293         fhMCLambda0vsClusterMaxCellDiffE0[i]->SetXTitle("#lambda_{0}^{2}");
1294         fhMCLambda0vsClusterMaxCellDiffE0[i]->SetYTitle("(E_{cluster} - E_{cell max})/ E_{cluster}");
1295         outputContainer->Add(fhMCLambda0vsClusterMaxCellDiffE0[i]) ;
1296
1297         fhMCLambda0vsClusterMaxCellDiffE2[i]  = new TH2F(Form("hLambda0vsClusterMaxCellDiffE2_MC%s",pnamess[i].Data()),
1298                                                          Form("cluster from %s : #lambda^{2}_{0} vs fraction of energy carried by max cell, 2< E < 6 GeV",ptypess[i].Data()),
1299                                                          ssbins,ssmin,ssmax,500,0,1.);
1300         fhMCLambda0vsClusterMaxCellDiffE2[i]->SetXTitle("#lambda_{0}^{2}");
1301         fhMCLambda0vsClusterMaxCellDiffE2[i]->SetYTitle("(E_{cluster} - E_{cell max})/ E_{cluster}");
1302         outputContainer->Add(fhMCLambda0vsClusterMaxCellDiffE2[i]) ;
1303
1304         fhMCLambda0vsClusterMaxCellDiffE6[i]  = new TH2F(Form("hLambda0vsClusterMaxCellDiffE6_MC%s",pnamess[i].Data()),
1305                                                          Form("cluster from %s : #lambda^{2}_{0} vs fraction of energy carried by max cell, E > 6 GeV",ptypess[i].Data()),
1306                                                          ssbins,ssmin,ssmax,500,0,1.);
1307         fhMCLambda0vsClusterMaxCellDiffE6[i]->SetXTitle("#lambda_{0}^{2}");
1308         fhMCLambda0vsClusterMaxCellDiffE6[i]->SetYTitle("(E_{cluster} - E_{cell max})/ E_{cluster}");
1309         outputContainer->Add(fhMCLambda0vsClusterMaxCellDiffE6[i]) ;
1310
1311         fhMCNCellsvsClusterMaxCellDiffE0[i]  = new TH2F(Form("hNCellsvsClusterMaxCellDiffE0_MC%s",pnamess[i].Data()),
1312                                                          Form("cluster from %s : N cells in cluster vs fraction of energy carried by max cell, E < 2 GeV",ptypess[i].Data()),
1313                                                          nbins/5,nmin,nmax/5,500,0,1.);
1314         fhMCNCellsvsClusterMaxCellDiffE0[i]->SetXTitle("N cells in cluster");
1315         fhMCNCellsvsClusterMaxCellDiffE0[i]->SetYTitle("(E_{cluster} - E_{cell max})/ E_{cluster}");
1316         outputContainer->Add(fhMCNCellsvsClusterMaxCellDiffE0[i]) ;
1317
1318         fhMCNCellsvsClusterMaxCellDiffE2[i]  = new TH2F(Form("hNCellsvsClusterMaxCellDiffE2_MC%s",pnamess[i].Data()),
1319                                                          Form("cluster from %s : N cells in cluster vs fraction of energy carried by max cell, 2< E < 6 GeV",ptypess[i].Data()),
1320                                                          nbins/5,nmin,nmax/5,500,0,1.);
1321         fhMCNCellsvsClusterMaxCellDiffE2[i]->SetXTitle("N cells in cluster");
1322         fhMCNCellsvsClusterMaxCellDiffE2[i]->SetYTitle("(E_{cluster} - E_{cell max})/ E_{cluster}");
1323         outputContainer->Add(fhMCNCellsvsClusterMaxCellDiffE2[i]) ;
1324
1325         fhMCNCellsvsClusterMaxCellDiffE6[i]  = new TH2F(Form("hNCellsvsClusterMaxCellDiffE6_MC%s",pnamess[i].Data()),
1326                                                          Form("cluster from %s : N cells in cluster vs fraction of energy carried by max cell, E > 6 GeV",ptypess[i].Data()),
1327                                                          nbins/5,nmin,nmax/5,500,0,1.);
1328         fhMCNCellsvsClusterMaxCellDiffE6[i]->SetXTitle("N cells in cluster");
1329         fhMCNCellsvsClusterMaxCellDiffE6[i]->SetYTitle("E (GeV)");
1330         outputContainer->Add(fhMCNCellsvsClusterMaxCellDiffE6[i]) ;
1331
1332        }// loop
1333
1334       if(!GetReader()->IsEmbeddedClusterSelectionOn())
1335       {
1336         fhMCPhotonELambda0NoOverlap  = new TH2F("hELambda0_MCPhoton_NoOverlap",
1337                                                 "cluster from Photon : E vs #lambda_{0}^{2}",
1338                                                 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1339         fhMCPhotonELambda0NoOverlap->SetYTitle("#lambda_{0}^{2}");
1340         fhMCPhotonELambda0NoOverlap->SetXTitle("E (GeV)");
1341         outputContainer->Add(fhMCPhotonELambda0NoOverlap) ;
1342
1343         fhMCPhotonELambda0TwoOverlap  = new TH2F("hELambda0_MCPhoton_TwoOverlap",
1344                                                  "cluster from Photon : E vs #lambda_{0}^{2}",
1345                                                  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1346         fhMCPhotonELambda0TwoOverlap->SetYTitle("#lambda_{0}^{2}");
1347         fhMCPhotonELambda0TwoOverlap->SetXTitle("E (GeV)");
1348         outputContainer->Add(fhMCPhotonELambda0TwoOverlap) ;
1349
1350         fhMCPhotonELambda0NOverlap  = new TH2F("hELambda0_MCPhoton_NOverlap",
1351                                                "cluster from Photon : E vs #lambda_{0}^{2}",
1352                                                nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1353         fhMCPhotonELambda0NOverlap->SetYTitle("#lambda_{0}^{2}");
1354         fhMCPhotonELambda0NOverlap->SetXTitle("E (GeV)");
1355         outputContainer->Add(fhMCPhotonELambda0NOverlap) ;
1356
1357       } //No embedding
1358
1359       //Fill histograms to check shape of embedded clusters
1360       if(GetReader()->IsEmbeddedClusterSelectionOn())
1361       {
1362
1363         fhEmbeddedSignalFractionEnergy  = new TH2F("hEmbeddedSignal_FractionEnergy",
1364                                                     "Energy Fraction of embedded signal versus cluster energy",
1365                                                     nptbins,ptmin,ptmax,100,0.,1.);
1366         fhEmbeddedSignalFractionEnergy->SetYTitle("Fraction");
1367         fhEmbeddedSignalFractionEnergy->SetXTitle("E (GeV)");
1368         outputContainer->Add(fhEmbeddedSignalFractionEnergy) ;
1369
1370         fhEmbedPhotonELambda0FullSignal  = new TH2F("hELambda0_EmbedPhoton_FullSignal",
1371                                                 "cluster from Photon embedded with more than 90% energy in cluster : E vs #lambda_{0}^{2}",
1372                                                 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1373         fhEmbedPhotonELambda0FullSignal->SetYTitle("#lambda_{0}^{2}");
1374         fhEmbedPhotonELambda0FullSignal->SetXTitle("E (GeV)");
1375         outputContainer->Add(fhEmbedPhotonELambda0FullSignal) ;
1376
1377         fhEmbedPhotonELambda0MostlySignal  = new TH2F("hELambda0_EmbedPhoton_MostlySignal",
1378                                           "cluster from Photon embedded with 50% to 90% energy in cluster : E vs #lambda_{0}^{2}",
1379                                           nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1380         fhEmbedPhotonELambda0MostlySignal->SetYTitle("#lambda_{0}^{2}");
1381         fhEmbedPhotonELambda0MostlySignal->SetXTitle("E (GeV)");
1382         outputContainer->Add(fhEmbedPhotonELambda0MostlySignal) ;
1383
1384         fhEmbedPhotonELambda0MostlyBkg  = new TH2F("hELambda0_EmbedPhoton_MostlyBkg",
1385                                           "cluster from Photon embedded with 10% to 50% energy in cluster : E vs #lambda_{0}^{2}",
1386                                           nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1387         fhEmbedPhotonELambda0MostlyBkg->SetYTitle("#lambda_{0}^{2}");
1388         fhEmbedPhotonELambda0MostlyBkg->SetXTitle("E (GeV)");
1389         outputContainer->Add(fhEmbedPhotonELambda0MostlyBkg) ;
1390
1391         fhEmbedPhotonELambda0FullBkg  = new TH2F("hELambda0_EmbedPhoton_FullBkg",
1392                                                    "cluster from Photonm embedded with 0% to 10% energy in cluster : E vs #lambda_{0}^{2}",
1393                                                    nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1394         fhEmbedPhotonELambda0FullBkg->SetYTitle("#lambda_{0}^{2}");
1395         fhEmbedPhotonELambda0FullBkg->SetXTitle("E (GeV)");
1396         outputContainer->Add(fhEmbedPhotonELambda0FullBkg) ;
1397
1398         fhEmbedPi0ELambda0FullSignal  = new TH2F("hELambda0_EmbedPi0_FullSignal",
1399                                                     "cluster from Pi0 embedded with more than 90% energy in cluster : E vs #lambda_{0}^{2}",
1400                                                     nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1401         fhEmbedPi0ELambda0FullSignal->SetYTitle("#lambda_{0}^{2}");
1402         fhEmbedPi0ELambda0FullSignal->SetXTitle("E (GeV)");
1403         outputContainer->Add(fhEmbedPi0ELambda0FullSignal) ;
1404
1405         fhEmbedPi0ELambda0MostlySignal  = new TH2F("hELambda0_EmbedPi0_MostlySignal",
1406                                                       "cluster from Pi0 embedded with 50% to 90% energy in cluster : E vs #lambda_{0}^{2}",
1407                                                       nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1408         fhEmbedPi0ELambda0MostlySignal->SetYTitle("#lambda_{0}^{2}");
1409         fhEmbedPi0ELambda0MostlySignal->SetXTitle("E (GeV)");
1410         outputContainer->Add(fhEmbedPi0ELambda0MostlySignal) ;
1411
1412         fhEmbedPi0ELambda0MostlyBkg  = new TH2F("hELambda0_EmbedPi0_MostlyBkg",
1413                                                    "cluster from Pi0 embedded with 10% to 50% energy in cluster : E vs #lambda_{0}^{2}",
1414                                                    nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1415         fhEmbedPi0ELambda0MostlyBkg->SetYTitle("#lambda_{0}^{2}");
1416         fhEmbedPi0ELambda0MostlyBkg->SetXTitle("E (GeV)");
1417         outputContainer->Add(fhEmbedPi0ELambda0MostlyBkg) ;
1418
1419         fhEmbedPi0ELambda0FullBkg  = new TH2F("hELambda0_EmbedPi0_FullBkg",
1420                                                  "cluster from Pi0 embedded with 0% to 10% energy in cluster : E vs #lambda_{0}^{2}",
1421                                                  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1422         fhEmbedPi0ELambda0FullBkg->SetYTitle("#lambda_{0}^{2}");
1423         fhEmbedPi0ELambda0FullBkg->SetXTitle("E (GeV)");
1424         outputContainer->Add(fhEmbedPi0ELambda0FullBkg) ;
1425
1426       }// embedded histograms
1427
1428
1429     }// Fill SS MC histograms
1430
1431   }//Histos with MC
1432
1433   //Store calo PID histograms
1434   if(fRejectTrackMatch){
1435     TList * caloPIDHistos = GetCaloPID()->GetCreateOutputObjects() ;
1436     for(Int_t i = 0; i < caloPIDHistos->GetEntries(); i++) {
1437       outputContainer->Add(caloPIDHistos->At(i)) ;
1438     }
1439     delete caloPIDHistos;
1440   }
1441
1442   return outputContainer ;
1443
1444 }
1445
1446 //____________________________________________________________________________
1447 void AliAnaPhoton::Init()
1448 {
1449
1450   //Init
1451   //Do some checks
1452   if(fCalorimeter == "PHOS" && !GetReader()->IsPHOSSwitchedOn() && NewOutputAOD()){
1453     printf("AliAnaPhoton::Init() - !!STOP: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!\n");
1454     abort();
1455   }
1456   else  if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn() && NewOutputAOD()){
1457     printf("AliAnaPhoton::Init() - !!STOP: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!\n");
1458     abort();
1459   }
1460
1461   if(GetReader()->GetDataType() == AliCaloTrackReader::kMC) GetCaloPID()->SwitchOnBayesian();
1462
1463 }
1464
1465 //____________________________________________________________________________
1466 void AliAnaPhoton::InitParameters()
1467 {
1468
1469   //Initialize the parameters of the analysis.
1470   AddToHistogramsName("AnaPhoton_");
1471
1472   fCalorimeter = "EMCAL" ;
1473   fMinDist     = 2.;
1474   fMinDist2    = 4.;
1475   fMinDist3    = 5.;
1476
1477   fTimeCutMin  = -1;
1478   fTimeCutMax  = 9999999;
1479   fNCellsCut   = 0;
1480
1481   fRejectTrackMatch       = kTRUE ;
1482
1483 }
1484
1485 //__________________________________________________________________
1486 void  AliAnaPhoton::MakeAnalysisFillAOD()
1487 {
1488   //Do photon analysis and fill aods
1489
1490   //Get the vertex
1491   Double_t v[3] = {0,0,0}; //vertex ;
1492   GetReader()->GetVertex(v);
1493
1494   //Select the Calorimeter of the photon
1495   TObjArray * pl = 0x0;
1496   if(fCalorimeter == "PHOS")
1497     pl = GetPHOSClusters();
1498   else if (fCalorimeter == "EMCAL")
1499     pl = GetEMCALClusters();
1500
1501   if(!pl) {
1502     Info("MakeAnalysisFillAOD","TObjArray with %s clusters is NULL!\n",fCalorimeter.Data());
1503     return;
1504   }
1505
1506   //Init arrays, variables, get number of clusters
1507   TLorentzVector mom, mom2 ;
1508   Int_t nCaloClusters = pl->GetEntriesFast();
1509
1510   if(GetDebug() > 0) printf("AliAnaPhoton::MakeAnalysisFillAOD() - input %s cluster entries %d\n", fCalorimeter.Data(), nCaloClusters);
1511
1512   //----------------------------------------------------
1513   // Fill AOD with PHOS/EMCAL AliAODPWG4Particle objects
1514   //----------------------------------------------------
1515   // Loop on clusters
1516   for(Int_t icalo = 0; icalo < nCaloClusters; icalo++){
1517
1518           AliVCluster * calo =  (AliVCluster*) (pl->At(icalo));
1519     //printf("calo %d, %f\n",icalo,calo->E());
1520
1521     //Get the index where the cluster comes, to retrieve the corresponding vertex
1522     Int_t evtIndex = 0 ;
1523     if (GetMixedEvent()) {
1524       evtIndex=GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ;
1525       //Get the vertex and check it is not too large in z
1526       if(TMath::Abs(GetVertex(evtIndex)[2])> GetZvertexCut()) continue;
1527     }
1528
1529     //Cluster selection, not charged, with photon id and in fiducial cut
1530     if(GetReader()->GetDataType() != AliCaloTrackReader::kMC){
1531       calo->GetMomentum(mom,GetVertex(evtIndex)) ;}//Assume that come from vertex in straight line
1532     else{
1533       Double_t vertex[]={0,0,0};
1534       calo->GetMomentum(mom,vertex) ;
1535     }
1536
1537     //--------------------------------------
1538     // Cluster selection
1539     //--------------------------------------
1540     if(!ClusterSelected(calo,mom)) continue;
1541
1542     //----------------------------
1543     //Create AOD for analysis
1544     //----------------------------
1545     AliAODPWG4Particle aodph = AliAODPWG4Particle(mom);
1546
1547     //...............................................
1548     //Set the indeces of the original caloclusters (MC, ID), and calorimeter
1549     Int_t label = calo->GetLabel();
1550     aodph.SetLabel(label);
1551     aodph.SetCaloLabel(calo->GetID(),-1);
1552     aodph.SetDetector(fCalorimeter);
1553     //printf("Index %d, Id %d, iaod %d\n",icalo, calo->GetID(),GetOutputAODBranch()->GetEntriesFast());
1554
1555     //...............................................
1556     //Set bad channel distance bit
1557     Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
1558     if     (distBad > fMinDist3) aodph.SetDistToBad(2) ;
1559     else if(distBad > fMinDist2) aodph.SetDistToBad(1) ;
1560     else                         aodph.SetDistToBad(0) ;
1561     //printf("DistBad %f Bit %d\n",distBad, aodph.DistToBad());
1562
1563     //--------------------------------------------------------------------------------------
1564     //Play with the MC stack if available
1565     //--------------------------------------------------------------------------------------
1566
1567     //Check origin of the candidates
1568     if(IsDataMC()){
1569       aodph.SetTag(GetMCAnalysisUtils()->CheckOrigin(calo->GetLabels(),calo->GetNLabels(),GetReader(), aodph.GetInputFileIndex()));
1570
1571       if(GetDebug() > 0)
1572         printf("AliAnaPhoton::MakeAnalysisFillAOD() - Origin of candidate, bit map %d\n",aodph.GetTag());
1573     }//Work with stack also
1574
1575     //--------------------------------------------------------------------------------------
1576     //Fill some shower shape histograms before PID is applied
1577     //--------------------------------------------------------------------------------------
1578
1579     FillShowerShapeHistograms(calo,aodph.GetTag());
1580
1581     //-------------------------------------
1582     //PID selection or bit setting
1583     //-------------------------------------
1584
1585     //...............................................
1586     // Data, PID check on
1587     if(IsCaloPIDOn()){
1588       // Get most probable PID, 2 options check bayesian PID weights or redo PID
1589       // By default, redo PID
1590
1591       aodph.SetIdentifiedParticleType(GetCaloPID()->GetIdentifiedParticleType(fCalorimeter,mom,calo));
1592
1593       if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PDG of identified particle %d\n",aodph.GetIdentifiedParticleType());
1594
1595       //If cluster does not pass pid, not photon, skip it.
1596       if(aodph.GetIdentifiedParticleType() != AliCaloPID::kPhoton) continue ;
1597
1598     }
1599     //...............................................
1600     // Data, PID check off
1601     else{
1602       //Set PID bits for later selection (AliAnaPi0 for example)
1603       //GetIdentifiedParticleType already called in SetPIDBits.
1604
1605       GetCaloPID()->SetPIDBits(fCalorimeter,calo,&aodph, GetCaloUtils(),GetReader()->GetInputEvent());
1606
1607       if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PID Bits set \n");
1608     }
1609
1610     if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - Photon selection cuts passed: pT %3.2f, pdg %d\n",aodph.Pt(), aodph.GetIdentifiedParticleType());
1611
1612     //Add AOD with photon object to aod branch
1613     AddAODParticle(aodph);
1614
1615   }//loop
1616
1617   if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD()  End fill AODs, with %d entries \n",GetOutputAODBranch()->GetEntriesFast());
1618
1619 }
1620
1621 //__________________________________________________________________
1622 void  AliAnaPhoton::MakeAnalysisFillHistograms()
1623 {
1624   //Fill histograms
1625
1626   //-------------------------------------------------------------------
1627   // Access MC information in stack if requested, check that it exists.
1628   AliStack         * stack       = 0x0;
1629   TParticle        * primary     = 0x0;
1630   TClonesArray     * mcparticles = 0x0;
1631   AliAODMCParticle * aodprimary  = 0x0;
1632
1633   if(IsDataMC()){
1634
1635     if(GetReader()->ReadStack()){
1636       stack =  GetMCStack() ;
1637       if(!stack) {
1638         printf("AliAnaPhoton::MakeAnalysisFillHistograms() - Stack not available, is the MC handler called? STOP\n");
1639         abort();
1640       }
1641
1642     }
1643     else if(GetReader()->ReadAODMCParticles()){
1644
1645       //Get the list of MC particles
1646       mcparticles = GetReader()->GetAODMCParticles(0);
1647       if(!mcparticles && GetDebug() > 0)        {
1648         printf("AliAnaPhoton::MakeAnalysisFillHistograms() -  Standard MCParticles not available!\n");
1649       }
1650     }
1651   }// is data and MC
1652
1653
1654   // Get vertex
1655   Double_t v[3] = {0,0,0}; //vertex ;
1656   GetReader()->GetVertex(v);
1657   //fhVertex->Fill(v[0],v[1],v[2]);
1658   if(TMath::Abs(v[2]) > GetZvertexCut()) return ; // done elsewhere for Single Event analysis, but there for mixed event
1659
1660   //----------------------------------
1661   //Loop on stored AOD photons
1662   Int_t naod = GetOutputAODBranch()->GetEntriesFast();
1663   if(GetDebug() > 0) printf("AliAnaPhoton::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod);
1664
1665   for(Int_t iaod = 0; iaod < naod ; iaod++){
1666     AliAODPWG4Particle* ph =  (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
1667     Int_t pdg = ph->GetIdentifiedParticleType();
1668
1669     if(GetDebug() > 3)
1670       printf("AliAnaPhoton::MakeAnalysisFillHistograms() - PDG %d, MC TAG %d, Calorimeter %s\n", ph->GetIdentifiedParticleType(),ph->GetTag(), (ph->GetDetector()).Data()) ;
1671
1672     //If PID used, fill histos with photons in Calorimeter fCalorimeter
1673     if(IsCaloPIDOn() && pdg != AliCaloPID::kPhoton) continue;
1674     if(ph->GetDetector() != fCalorimeter) continue;
1675
1676     if(GetDebug() > 2)
1677       printf("AliAnaPhoton::MakeAnalysisFillHistograms() - ID Photon: pt %f, phi %f, eta %f\n", ph->Pt(),ph->Phi(),ph->Eta()) ;
1678
1679     //................................
1680     //Fill photon histograms
1681     Float_t ptcluster  = ph->Pt();
1682     Float_t phicluster = ph->Phi();
1683     Float_t etacluster = ph->Eta();
1684     Float_t ecluster   = ph->E();
1685
1686     fhEPhoton   ->Fill(ecluster);
1687     fhPtPhoton  ->Fill(ptcluster);
1688     fhPhiPhoton ->Fill(ptcluster,phicluster);
1689     fhEtaPhoton ->Fill(ptcluster,etacluster);
1690     if     (ecluster > 0.5)   fhEtaPhiPhoton  ->Fill(etacluster, phicluster);
1691     else if(GetMinPt() < 0.5) fhEtaPhi05Photon->Fill(etacluster, phicluster);
1692
1693
1694     //Get original cluster, to recover some information
1695     Int_t absID             = 0;
1696     Float_t maxCellFraction = 0;
1697     AliVCaloCells* cells    = 0;
1698     TObjArray * clusters    = 0;
1699     if(fCalorimeter == "EMCAL"){
1700       cells    = GetEMCALCells();
1701       clusters = GetEMCALClusters();
1702     }
1703     else{
1704       cells    = GetPHOSCells();
1705       clusters = GetPHOSClusters();
1706     }
1707
1708     Int_t iclus = -1;
1709     AliVCluster *cluster = FindCluster(clusters,ph->GetCaloLabel(0),iclus);
1710     absID = GetCaloUtils()->GetMaxEnergyCell(cells, cluster,maxCellFraction);
1711
1712     // Control histograms
1713     fhMaxCellDiffClusterE->Fill(ph->E(),maxCellFraction);
1714     fhNCellsE            ->Fill(ph->E(),cluster->GetNCells());
1715     fhTimeE              ->Fill(ph->E(),cluster->GetTOF()*1.e9);
1716
1717     //.......................................
1718     //Play with the MC data if available
1719     if(IsDataMC()){
1720
1721       FillAcceptanceHistograms();
1722
1723       //....................................................................
1724       // Access MC information in stack if requested, check that it exists.
1725       Int_t label =ph->GetLabel();
1726       if(label < 0) {
1727         if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** bad label ***:  label %d \n", label);
1728         continue;
1729       }
1730
1731       Float_t eprim   = 0;
1732       Float_t ptprim  = 0;
1733       if(GetReader()->ReadStack()){
1734
1735         if(label >=  stack->GetNtrack()) {
1736           if(GetDebug() > 2)  printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** large label ***:  label %d, n tracks %d \n", label, stack->GetNtrack());
1737           continue ;
1738         }
1739
1740         primary = stack->Particle(label);
1741         if(!primary){
1742           printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** no primary ***:  label %d \n", label);
1743           continue;
1744         }
1745         eprim   = primary->Energy();
1746         ptprim  = primary->Pt();
1747
1748       }
1749       else if(GetReader()->ReadAODMCParticles()){
1750         //Check which is the input
1751         if(ph->GetInputFileIndex() == 0){
1752           if(!mcparticles) continue;
1753           if(label >=  mcparticles->GetEntriesFast()) {
1754             if(GetDebug() > 2)  printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** large label ***:  label %d, n tracks %d \n",
1755                                        label, mcparticles->GetEntriesFast());
1756             continue ;
1757           }
1758           //Get the particle
1759           aodprimary = (AliAODMCParticle*) mcparticles->At(label);
1760
1761         }
1762
1763         if(!aodprimary){
1764           printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** no primary ***:  label %d \n", label);
1765           continue;
1766         }
1767
1768         eprim   = aodprimary->E();
1769         ptprim  = aodprimary->Pt();
1770
1771       }
1772
1773       Int_t tag =ph->GetTag();
1774
1775       if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) && fhMCE[mcPhoton])
1776       {
1777         fhMCE  [mcPhoton] ->Fill(ecluster);
1778         fhMCPt [mcPhoton] ->Fill(ptcluster);
1779         fhMCPhi[mcPhoton] ->Fill(ecluster,phicluster);
1780         fhMCEta[mcPhoton] ->Fill(ecluster,etacluster);
1781
1782         fhMC2E[mcPhoton]     ->Fill(ecluster, eprim);
1783         fhMC2Pt[mcPhoton]    ->Fill(ptcluster, ptprim);
1784         fhMCDeltaE[mcPhoton] ->Fill(ecluster,eprim-ecluster);
1785         fhMCDeltaPt[mcPhoton]->Fill(ptcluster,ptprim-ptcluster);
1786
1787         if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion) && fhMCE[mcConversion])
1788         {
1789           fhMCE  [mcConversion] ->Fill(ecluster);
1790           fhMCPt [mcConversion] ->Fill(ptcluster);
1791           fhMCPhi[mcConversion] ->Fill(ecluster,phicluster);
1792           fhMCEta[mcConversion] ->Fill(ecluster,etacluster);
1793
1794           fhMC2E[mcConversion]     ->Fill(ecluster, eprim);
1795           fhMC2Pt[mcConversion]    ->Fill(ptcluster, ptprim);
1796           fhMCDeltaE[mcConversion] ->Fill(ecluster,eprim-ecluster);
1797           fhMCDeltaPt[mcConversion]->Fill(ptcluster,ptprim-ptcluster);
1798         }
1799
1800         if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt) && fhMCE[mcPrompt]){
1801           fhMCE  [mcPrompt] ->Fill(ecluster);
1802           fhMCPt [mcPrompt] ->Fill(ptcluster);
1803           fhMCPhi[mcPrompt] ->Fill(ecluster,phicluster);
1804           fhMCEta[mcPrompt] ->Fill(ecluster,etacluster);
1805
1806           fhMC2E[mcPrompt]     ->Fill(ecluster, eprim);
1807           fhMC2Pt[mcPrompt]    ->Fill(ptcluster, ptprim);
1808           fhMCDeltaE[mcPrompt] ->Fill(ecluster,eprim-ecluster);
1809           fhMCDeltaPt[mcPrompt]->Fill(ptcluster,ptprim-ptcluster);
1810
1811         }
1812         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation)&& fhMCE[mcFragmentation])
1813         {
1814           fhMCE  [mcFragmentation] ->Fill(ecluster);
1815           fhMCPt [mcFragmentation] ->Fill(ptcluster);
1816           fhMCPhi[mcFragmentation] ->Fill(ecluster,phicluster);
1817           fhMCEta[mcFragmentation] ->Fill(ecluster,etacluster);
1818
1819           fhMC2E[mcFragmentation]     ->Fill(ecluster, eprim);
1820           fhMC2Pt[mcFragmentation]    ->Fill(ptcluster, ptprim);
1821           fhMCDeltaE[mcFragmentation] ->Fill(ecluster,eprim-ecluster);
1822           fhMCDeltaPt[mcFragmentation]->Fill(ptcluster,ptprim-ptcluster);
1823
1824         }
1825         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCISR)&& fhMCE[mcISR])
1826         {
1827           fhMCE  [mcISR] ->Fill(ecluster);
1828           fhMCPt [mcISR] ->Fill(ptcluster);
1829           fhMCPhi[mcISR] ->Fill(ecluster,phicluster);
1830           fhMCEta[mcISR] ->Fill(ecluster,etacluster);
1831
1832           fhMC2E[mcISR]     ->Fill(ecluster, eprim);
1833           fhMC2Pt[mcISR]    ->Fill(ptcluster, ptprim);
1834           fhMCDeltaE[mcISR] ->Fill(ecluster,eprim-ecluster);
1835           fhMCDeltaPt[mcISR]->Fill(ptcluster,ptprim-ptcluster);
1836
1837         }
1838         else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay) &&
1839                 !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) && fhMCE[mcPi0Decay])
1840         {
1841           fhMCE  [mcPi0Decay] ->Fill(ecluster);
1842           fhMCPt [mcPi0Decay] ->Fill(ptcluster);
1843           fhMCPhi[mcPi0Decay] ->Fill(ecluster,phicluster);
1844           fhMCEta[mcPi0Decay] ->Fill(ecluster,etacluster);
1845
1846           fhMC2E[mcPi0Decay]     ->Fill(ecluster, eprim);
1847           fhMC2Pt[mcPi0Decay]    ->Fill(ptcluster, ptprim);
1848           fhMCDeltaE[mcPi0Decay] ->Fill(ecluster,eprim-ecluster);
1849           fhMCDeltaPt[mcPi0Decay]->Fill(ptcluster,ptprim-ptcluster);
1850
1851         }
1852         else if( (GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay) ||
1853                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay) ) && fhMCE[mcOtherDecay])
1854         {
1855           fhMCE  [mcOtherDecay] ->Fill(ecluster);
1856           fhMCPt [mcOtherDecay] ->Fill(ptcluster);
1857           fhMCPhi[mcOtherDecay] ->Fill(ecluster,phicluster);
1858           fhMCEta[mcOtherDecay] ->Fill(ecluster,etacluster);
1859
1860           fhMC2E[mcOtherDecay]     ->Fill(ecluster, eprim);
1861           fhMC2Pt[mcOtherDecay]    ->Fill(ptcluster, ptprim);
1862           fhMCDeltaE[mcOtherDecay] ->Fill(ecluster,eprim-ecluster);
1863           fhMCDeltaPt[mcOtherDecay]->Fill(ptcluster,ptprim-ptcluster);
1864
1865         }
1866         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) && fhMCE  [mcPi0])
1867         {
1868           fhMCE  [mcPi0] ->Fill(ecluster);
1869           fhMCPt [mcPi0] ->Fill(ptcluster);
1870           fhMCPhi[mcPi0] ->Fill(ecluster,phicluster);
1871           fhMCEta[mcPi0] ->Fill(ecluster,etacluster);
1872
1873           fhMC2E[mcPi0]     ->Fill(ecluster, eprim);
1874           fhMC2Pt[mcPi0]    ->Fill(ptcluster, ptprim);
1875           fhMCDeltaE[mcPi0] ->Fill(ecluster,eprim-ecluster);
1876           fhMCDeltaPt[mcPi0]->Fill(ptcluster,ptprim-ptcluster);
1877
1878         }
1879         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta) && fhMCE[mcEta])
1880         {
1881           fhMCE  [mcEta] ->Fill(ecluster);
1882           fhMCPt [mcEta] ->Fill(ptcluster);
1883           fhMCPhi[mcEta] ->Fill(ecluster,phicluster);
1884           fhMCEta[mcEta] ->Fill(ecluster,etacluster);
1885
1886           fhMC2E[mcEta]     ->Fill(ecluster, eprim);
1887           fhMC2Pt[mcEta]    ->Fill(ptcluster, ptprim);
1888           fhMCDeltaE[mcEta] ->Fill(ecluster,eprim-ecluster);
1889           fhMCDeltaPt[mcEta]->Fill(ptcluster,ptprim-ptcluster);
1890
1891         }
1892       }
1893       else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiNeutron) && fhMCE[mcAntiNeutron])
1894       {
1895         fhMCE  [mcAntiNeutron] ->Fill(ecluster);
1896         fhMCPt [mcAntiNeutron] ->Fill(ptcluster);
1897         fhMCPhi[mcAntiNeutron] ->Fill(ecluster,phicluster);
1898         fhMCEta[mcAntiNeutron] ->Fill(ecluster,etacluster);
1899
1900         fhMC2E[mcAntiNeutron]     ->Fill(ecluster, eprim);
1901         fhMC2Pt[mcAntiNeutron]    ->Fill(ptcluster, ptprim);
1902         fhMCDeltaE[mcAntiNeutron] ->Fill(ecluster,eprim-ecluster);
1903         fhMCDeltaPt[mcAntiNeutron]->Fill(ptcluster,ptprim-ptcluster);
1904
1905       }
1906       else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiProton) && fhMCE[mcAntiProton])
1907       {
1908         fhMCE  [mcAntiProton] ->Fill(ecluster);
1909         fhMCPt [mcAntiProton] ->Fill(ptcluster);
1910         fhMCPhi[mcAntiProton] ->Fill(ecluster,phicluster);
1911         fhMCEta[mcAntiProton] ->Fill(ecluster,etacluster);
1912
1913         fhMC2E[mcAntiProton]     ->Fill(ecluster, eprim);
1914         fhMC2Pt[mcAntiProton]    ->Fill(ptcluster, ptprim);
1915         fhMCDeltaE[mcAntiProton] ->Fill(ecluster,eprim-ecluster);
1916         fhMCDeltaPt[mcAntiProton]->Fill(ecluster,ptprim-ptcluster);
1917
1918       }
1919       else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron) && fhMCE[mcElectron])
1920       {
1921         fhMCE  [mcElectron] ->Fill(ecluster);
1922         fhMCPt [mcElectron] ->Fill(ptcluster);
1923         fhMCPhi[mcElectron] ->Fill(ecluster,phicluster);
1924         fhMCEta[mcElectron] ->Fill(ecluster,etacluster);
1925
1926         fhMC2E[mcElectron]     ->Fill(ecluster, eprim);
1927         fhMC2Pt[mcElectron]    ->Fill(ptcluster, ptprim);
1928         fhMCDeltaE[mcElectron] ->Fill(ecluster,eprim-ecluster);
1929         fhMCDeltaPt[mcElectron]->Fill(ecluster,ptprim-ptcluster);
1930       }
1931       else if( fhMCE[mcOther]){
1932         fhMCE  [mcOther] ->Fill(ecluster);
1933         fhMCPt [mcOther] ->Fill(ptcluster);
1934         fhMCPhi[mcOther] ->Fill(ecluster,phicluster);
1935         fhMCEta[mcOther] ->Fill(ecluster,etacluster);
1936
1937         fhMC2E[mcOther]     ->Fill(ecluster, eprim);
1938         fhMC2Pt[mcOther]    ->Fill(ptcluster, ptprim);
1939         fhMCDeltaE[mcOther] ->Fill(ecluster,eprim-ecluster);
1940         fhMCDeltaPt[mcOther]->Fill(ecluster,ptprim-ptcluster);
1941
1942         //               printf(" AliAnaPhoton::MakeAnalysisFillHistograms() - Label %d, pT %2.3f Unknown, bits set: ",
1943         //                                      ph->GetLabel(),ph->Pt());
1944         //                for(Int_t i = 0; i < 20; i++) {
1945         //                        if(GetMCAnalysisUtils()->CheckTagBit(tag,i)) printf(" %d, ",i);
1946         //                }
1947         //                printf("\n");
1948
1949       }
1950
1951     }//Histograms with MC
1952
1953   }// aod loop
1954
1955 }
1956
1957
1958 //__________________________________________________________________
1959 void AliAnaPhoton::Print(const Option_t * opt) const
1960 {
1961   //Print some relevant parameters set for the analysis
1962
1963   if(! opt)
1964     return;
1965
1966   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
1967   AliAnaPartCorrBaseClass::Print(" ");
1968
1969   printf("Calorimeter            =     %s\n", fCalorimeter.Data()) ;
1970   printf("Min Distance to Bad Channel   = %2.1f\n",fMinDist);
1971   printf("Min Distance to Bad Channel 2 = %2.1f\n",fMinDist2);
1972   printf("Min Distance to Bad Channel 3 = %2.1f\n",fMinDist3);
1973   printf("Reject clusters with a track matched = %d\n",fRejectTrackMatch);
1974   printf("Time Cut: %3.1f < TOF  < %3.1f\n", fTimeCutMin, fTimeCutMax);
1975   printf("Number of cells in cluster is        > %d \n", fNCellsCut);
1976   printf("    \n") ;
1977
1978 }