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