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