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
 693     if(cluster->IsEMCAL() && GetCaloUtils()->IsRecalculationOfClusterTrackMatchingOn()){
 694       dR = 2000., dZ = 2000.;
 695       GetCaloUtils()->GetEMCALRecoUtils()->GetMatchedResiduals(cluster->GetID(),dZ,dR);
 696     }
 697
 698     if(TMath::Abs(dZ) < 0.05 && TMath::Abs(dR) < 0.05)
 699     {
 700       fhLam0ETM ->Fill(energy,lambda0);
 701       fhLam1ETM ->Fill(energy,lambda1);
 702       fhDispETM ->Fill(energy,disp);
 703
 704       if(fCalorimeter == "EMCAL" && GetModuleNumber(cluster) > 5){
 705         fhLam0ETMTRD->Fill(energy,lambda0);
 706         fhLam1ETMTRD->Fill(energy,lambda1);
 707         fhDispETMTRD->Fill(energy,disp);
 708       }
 709     }
 710   }// if track-matching was of, check effect of matching residual cut
 711
 712   if(energy < 2){
 713     fhNCellsLam0LowE ->Fill(ncells,lambda0);
 714     fhNCellsLam1LowE ->Fill(ncells,lambda1);
 715     fhNCellsDispLowE ->Fill(ncells,disp);
 716
 717     fhLam1Lam0LowE  ->Fill(lambda1,lambda0);
 718     fhLam0DispLowE  ->Fill(lambda0,disp);
 719     fhDispLam1LowE  ->Fill(disp,lambda1);
 720     fhEtaLam0LowE   ->Fill(eta,lambda0);
 721     fhPhiLam0LowE   ->Fill(phi,lambda0);
 722
 723   }
 724   else {
 725     fhNCellsLam0HighE ->Fill(ncells,lambda0);
 726     fhNCellsLam1HighE ->Fill(ncells,lambda1);
 727     fhNCellsDispHighE ->Fill(ncells,disp);
 728
 729     fhLam1Lam0HighE  ->Fill(lambda1,lambda0);
 730     fhLam0DispHighE  ->Fill(lambda0,disp);
 731     fhDispLam1HighE  ->Fill(disp,lambda1);
 732     fhEtaLam0HighE   ->Fill(eta, lambda0);
 733     fhPhiLam0HighE   ->Fill(phi, lambda0);
 734   }
 735
 736   if(IsDataMC()){
 737
 738     AliVCaloCells* cells = 0;
 739     if(fCalorimeter == "EMCAL") cells = GetEMCALCells();
 740     else                        cells = GetPHOSCells();
 741
 742     //Fill histograms to check shape of embedded clusters
 743     Float_t fraction = 0;
 744     if(GetReader()->IsEmbeddedClusterSelectionOn()){//Only working for EMCAL
 745
 746       Float_t clusterE = 0; // recalculate in case corrections applied.
 747       Float_t cellE    = 0;
 748       for(Int_t icell = 0; icell < cluster->GetNCells(); icell++){
 749         cellE    = cells->GetCellAmplitude(cluster->GetCellAbsId(icell));
 750         clusterE+=cellE;
 751         fraction+=cellE*cluster->GetCellAmplitudeFraction(icell);
 752       }
 753
 754       //Fraction of total energy due to the embedded signal
 755       fraction/=clusterE;
 756
 757       if(GetDebug() > 1 )
 758         printf("AliAnaPhoton::FillShowerShapeHistogram() - Energy fraction of embedded signal %2.3f, Energy %2.3f\n",fraction, clusterE);
 759
 760       fhEmbeddedSignalFractionEnergy->Fill(clusterE,fraction);
 761
 762     }  // embedded fraction
 763
 764     // Get the fraction of the cluster energy that carries the cell with highest energy
 765     Int_t absID             =-1 ;
 766     Float_t maxCellFraction = 0.;
 767
 768     absID = GetCaloUtils()->GetMaxEnergyCell(cells, cluster,maxCellFraction);
 769
 770     // Check the origin and fill histograms
 771     if( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton) &&
 772        !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCConversion) &&
 773        !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0) &&
 774        !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEta)){
 775       fhMCELambda0[kmcssPhoton]    ->Fill(energy, lambda0);
 776       fhMCELambda1[kmcssPhoton]    ->Fill(energy, lambda1);
 777       fhMCEDispersion[kmcssPhoton] ->Fill(energy, disp);
 778       fhMCNCellsE[kmcssPhoton]     ->Fill(energy, ncells);
 779       fhMCMaxCellDiffClusterE[kmcssPhoton]->Fill(energy,maxCellFraction);
 780
 781       if     (energy < 2.){
 782         fhMCLambda0vsClusterMaxCellDiffE0[kmcssPhoton]->Fill(lambda0, maxCellFraction);
 783         fhMCNCellsvsClusterMaxCellDiffE0[kmcssPhoton] ->Fill(ncells,  maxCellFraction);
 784       }
 785       else if(energy < 6.){
 786         fhMCLambda0vsClusterMaxCellDiffE2[kmcssPhoton]->Fill(lambda0, maxCellFraction);
 787         fhMCNCellsvsClusterMaxCellDiffE2[kmcssPhoton] ->Fill(ncells,  maxCellFraction);
 788       }
 789       else{
 790         fhMCLambda0vsClusterMaxCellDiffE6[kmcssPhoton]->Fill(lambda0, maxCellFraction);
 791         fhMCNCellsvsClusterMaxCellDiffE6[kmcssPhoton] ->Fill(ncells,  maxCellFraction);
 792       }
 793
 794       if(!GetReader()->IsEmbeddedClusterSelectionOn()){
 795         //Check particle overlaps in cluster
 796
 797         // Compare the primary depositing more energy with the rest,
 798         // if no photon/electron as comon ancestor (conversions), count as other particle
 799         Int_t ancPDG = 0, ancStatus = -1;
 800         TLorentzVector momentum; TVector3 prodVertex;
 801         Int_t ancLabel = 0;
 802         Int_t noverlaps = 1;
 803         for (UInt_t ilab = 0; ilab < cluster->GetNLabels(); ilab++ ) {
 804           ancLabel = GetMCAnalysisUtils()->CheckCommonAncestor(cluster->GetLabels()[0],cluster->GetLabels()[ilab], GetReader(),ancPDG,ancStatus,momentum,prodVertex);
 805           if(ancPDG!=22 && TMath::Abs(ancPDG)!=11) noverlaps++;
 806         }
 807         //printf("N overlaps %d \n",noverlaps);
 808
 809         if(noverlaps == 1){
 810           fhMCPhotonELambda0NoOverlap  ->Fill(energy, lambda0);
 811         }
 812         else if(noverlaps == 2){
 813           fhMCPhotonELambda0TwoOverlap ->Fill(energy, lambda0);
 814         }
 815         else if(noverlaps > 2){
 816           fhMCPhotonELambda0NOverlap   ->Fill(energy, lambda0);
 817         }
 818         else {
 819           printf("AliAnaPhoton::FillShowerShapeHistogram() - n overlaps = %d!!", noverlaps);
 820         }
 821       }//No embedding
 822
 823       //Fill histograms to check shape of embedded clusters
 824       if(GetReader()->IsEmbeddedClusterSelectionOn()){
 825
 826         if     (fraction > 0.9)
 827         {
 828           fhEmbedPhotonELambda0FullSignal   ->Fill(energy, lambda0);
 829         }
 830         else if(fraction > 0.5)
 831         {
 832           fhEmbedPhotonELambda0MostlySignal ->Fill(energy, lambda0);
 833         }
 834         else if(fraction > 0.1)
 835         {
 836           fhEmbedPhotonELambda0MostlyBkg    ->Fill(energy, lambda0);
 837         }
 838         else
 839         {
 840           fhEmbedPhotonELambda0FullBkg      ->Fill(energy, lambda0);
 841         }
 842       } // embedded
 843
 844     }//photon   no conversion
 845     else if  ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCElectron)){
 846       fhMCELambda0[kmcssElectron]    ->Fill(energy, lambda0);
 847       fhMCELambda1[kmcssElectron]    ->Fill(energy, lambda1);
 848       fhMCEDispersion[kmcssElectron] ->Fill(energy, disp);
 849       fhMCNCellsE[kmcssElectron]     ->Fill(energy, ncells);
 850       fhMCMaxCellDiffClusterE[kmcssElectron]->Fill(energy,maxCellFraction);
 851
 852       if     (energy < 2.){
 853         fhMCLambda0vsClusterMaxCellDiffE0[kmcssElectron]->Fill(lambda0, maxCellFraction);
 854         fhMCNCellsvsClusterMaxCellDiffE0[kmcssElectron] ->Fill(ncells,  maxCellFraction);
 855       }
 856       else if(energy < 6.){
 857         fhMCLambda0vsClusterMaxCellDiffE2[kmcssElectron]->Fill(lambda0, maxCellFraction);
 858         fhMCNCellsvsClusterMaxCellDiffE2[kmcssElectron] ->Fill(ncells,  maxCellFraction);
 859       }
 860       else{
 861         fhMCLambda0vsClusterMaxCellDiffE6[kmcssElectron]->Fill(lambda0, maxCellFraction);
 862         fhMCNCellsvsClusterMaxCellDiffE6[kmcssElectron] ->Fill(ncells,  maxCellFraction);
 863       }
 864     }//electron
 865     else if  ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton) &&
 866               GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCConversion) ){
 867       fhMCELambda0[kmcssConversion]    ->Fill(energy, lambda0);
 868       fhMCELambda1[kmcssConversion]    ->Fill(energy, lambda1);
 869       fhMCEDispersion[kmcssConversion] ->Fill(energy, disp);
 870       fhMCNCellsE[kmcssConversion]     ->Fill(energy, ncells);
 871       fhMCMaxCellDiffClusterE[kmcssConversion]->Fill(energy,maxCellFraction);
 872
 873       if     (energy < 2.){
 874         fhMCLambda0vsClusterMaxCellDiffE0[kmcssConversion]->Fill(lambda0, maxCellFraction);
 875         fhMCNCellsvsClusterMaxCellDiffE0[kmcssConversion] ->Fill(ncells,  maxCellFraction);
 876       }
 877       else if(energy < 6.){
 878         fhMCLambda0vsClusterMaxCellDiffE2[kmcssConversion]->Fill(lambda0, maxCellFraction);
 879         fhMCNCellsvsClusterMaxCellDiffE2[kmcssConversion] ->Fill(ncells,  maxCellFraction);
 880       }
 881       else{
 882         fhMCLambda0vsClusterMaxCellDiffE6[kmcssConversion]->Fill(lambda0, maxCellFraction);
 883         fhMCNCellsvsClusterMaxCellDiffE6[kmcssConversion] ->Fill(ncells,  maxCellFraction);
 884       }
 885
 886     }//conversion photon
 887     else if  ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0)  ){
 888       fhMCELambda0[kmcssPi0]    ->Fill(energy, lambda0);
 889       fhMCELambda1[kmcssPi0]    ->Fill(energy, lambda1);
 890       fhMCEDispersion[kmcssPi0] ->Fill(energy, disp);
 891       fhMCNCellsE[kmcssPi0]     ->Fill(energy, ncells);
 892       fhMCMaxCellDiffClusterE[kmcssPi0]->Fill(energy,maxCellFraction);
 893
 894       if     (energy < 2.){
 895         fhMCLambda0vsClusterMaxCellDiffE0[kmcssPi0]->Fill(lambda0, maxCellFraction);
 896         fhMCNCellsvsClusterMaxCellDiffE0[kmcssPi0] ->Fill(ncells,  maxCellFraction);
 897       }
 898       else if(energy < 6.){
 899         fhMCLambda0vsClusterMaxCellDiffE2[kmcssPi0]->Fill(lambda0, maxCellFraction);
 900         fhMCNCellsvsClusterMaxCellDiffE2[kmcssPi0] ->Fill(ncells,  maxCellFraction);
 901       }
 902       else{
 903         fhMCLambda0vsClusterMaxCellDiffE6[kmcssPi0]->Fill(lambda0, maxCellFraction);
 904         fhMCNCellsvsClusterMaxCellDiffE6[kmcssPi0] ->Fill(ncells,  maxCellFraction);
 905       }
 906
 907       //Fill histograms to check shape of embedded clusters
 908       if(GetReader()->IsEmbeddedClusterSelectionOn()){
 909
 910         if     (fraction > 0.9)
 911         {
 912           fhEmbedPi0ELambda0FullSignal   ->Fill(energy, lambda0);
 913         }
 914         else if(fraction > 0.5)
 915         {
 916           fhEmbedPi0ELambda0MostlySignal ->Fill(energy, lambda0);
 917         }
 918         else if(fraction > 0.1)
 919         {
 920           fhEmbedPi0ELambda0MostlyBkg    ->Fill(energy, lambda0);
 921         }
 922         else
 923         {
 924           fhEmbedPi0ELambda0FullBkg      ->Fill(energy, lambda0);
 925         }
 926       } // embedded
 927
 928     }//pi0
 929     else if  ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEta)  ){
 930       fhMCELambda0[kmcssEta]    ->Fill(energy, lambda0);
 931       fhMCELambda1[kmcssEta]    ->Fill(energy, lambda1);
 932       fhMCEDispersion[kmcssEta] ->Fill(energy, disp);
 933       fhMCNCellsE[kmcssEta]     ->Fill(energy, ncells);
 934       fhMCMaxCellDiffClusterE[kmcssEta]->Fill(energy,maxCellFraction);
 935
 936       if     (energy < 2.){
 937         fhMCLambda0vsClusterMaxCellDiffE0[kmcssEta]->Fill(lambda0, maxCellFraction);
 938         fhMCNCellsvsClusterMaxCellDiffE0[kmcssEta] ->Fill(ncells,  maxCellFraction);
 939       }
 940       else if(energy < 6.){
 941         fhMCLambda0vsClusterMaxCellDiffE2[kmcssEta]->Fill(lambda0, maxCellFraction);
 942         fhMCNCellsvsClusterMaxCellDiffE2[kmcssEta] ->Fill(ncells,  maxCellFraction);
 943       }
 944       else{
 945         fhMCLambda0vsClusterMaxCellDiffE6[kmcssEta]->Fill(lambda0, maxCellFraction);
 946         fhMCNCellsvsClusterMaxCellDiffE6[kmcssEta] ->Fill(ncells,  maxCellFraction);
 947       }
 948
 949     }//eta
 950     else {
 951       fhMCELambda0[kmcssOther]    ->Fill(energy, lambda0);
 952       fhMCELambda1[kmcssOther]    ->Fill(energy, lambda1);
 953       fhMCEDispersion[kmcssOther] ->Fill(energy, disp);
 954       fhMCNCellsE[kmcssOther]     ->Fill(energy, ncells);
 955       fhMCMaxCellDiffClusterE[kmcssOther]->Fill(energy,maxCellFraction);
 956
 957       if     (energy < 2.){
 958         fhMCLambda0vsClusterMaxCellDiffE0[kmcssOther]->Fill(lambda0, maxCellFraction);
 959         fhMCNCellsvsClusterMaxCellDiffE0[kmcssOther] ->Fill(ncells,  maxCellFraction);
 960       }
 961       else if(energy < 6.){
 962         fhMCLambda0vsClusterMaxCellDiffE2[kmcssOther]->Fill(lambda0, maxCellFraction);
 963         fhMCNCellsvsClusterMaxCellDiffE2[kmcssOther] ->Fill(ncells,  maxCellFraction);
 964       }
 965       else{
 966         fhMCLambda0vsClusterMaxCellDiffE6[kmcssOther]->Fill(lambda0, maxCellFraction);
 967         fhMCNCellsvsClusterMaxCellDiffE6[kmcssOther] ->Fill(ncells,  maxCellFraction);
 968       }
 969
 970     }//other particles
 971
 972   }//MC data
 973
 974 }
 975
 976 //__________________________________________________________________________
 977 void AliAnaPhoton::FillTrackMatchingResidualHistograms(AliVCluster* cluster,
 978                                                        const Int_t cut)
 979 {
 980   // If selected, fill histograms with residuals of matched clusters, help to define track matching cut
 981   // Residual filled for different cuts 0 (No cut), after 1 PID cut
 982
 983   Float_t dZ = cluster->GetTrackDz();
 984   Float_t dR = cluster->GetTrackDx();
 985
 986   if(cluster->IsEMCAL() && GetCaloUtils()->IsRecalculationOfClusterTrackMatchingOn())
 987   {
 988     dR = 2000., dZ = 2000.;
 989     GetCaloUtils()->GetEMCALRecoUtils()->GetMatchedResiduals(cluster->GetID(),dZ,dR);
 990   }
 991
 992   if(fhTrackMatchedDEta[cut] && TMath::Abs(dR) < 999)
 993   {
 994     fhTrackMatchedDEta[cut]->Fill(cluster->E(),dZ);
 995     fhTrackMatchedDPhi[cut]->Fill(cluster->E(),dR);
 996
 997     if(cluster->E() > 0.5) fhTrackMatchedDEtaDPhi[cut]->Fill(dZ,dR);
 998
 999     Int_t nSMod = GetModuleNumber(cluster);
1000
1001     if(fCalorimeter=="EMCAL" &&  nSMod > 5)
1002     {
1003       fhTrackMatchedDEtaTRD[cut]->Fill(cluster->E(),dZ);
1004       fhTrackMatchedDPhiTRD[cut]->Fill(cluster->E(),dR);
1005     }
1006
1007     // Check dEdx and E/p of matched clusters
1008
1009     if(TMath::Abs(dZ) < 0.05 && TMath::Abs(dR) < 0.05)
1010     {
1011
1012       AliVTrack *track = GetCaloUtils()->GetMatchedTrack(cluster, GetReader()->GetInputEvent());
1013
1014       if(track)
1015       {
1016
1017         Float_t dEdx   = track->GetTPCsignal();
1018         Float_t eOverp = cluster->E()/track->P();
1019
1020         fhdEdx[cut]  ->Fill(cluster->E(), dEdx);
1021         fhEOverP[cut]->Fill(cluster->E(), eOverp);
1022
1023         if(fCalorimeter=="EMCAL" && nSMod > 5)
1024           fhEOverPTRD[cut]->Fill(cluster->E(), eOverp);
1025
1026
1027       }
1028       else
1029           printf("AliAnaPhoton::FillTrackMatchingResidualHistograms() - Residual OK but (dR, dZ)= (%2.4f,%2.4f) no track associated WHAT? \n", dR,dZ);
1030
1031
1032
1033       if(IsDataMC())
1034       {
1035
1036         Int_t tag = GetMCAnalysisUtils()->CheckOrigin(cluster->GetLabels(),cluster->GetNLabels(),GetReader(), 0);
1037
1038         if  ( !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion)  )
1039         {
1040           if       ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0)      ||
1041                      GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta)       ) fhTrackMatchedMCParticle[cut]->Fill(cluster->E(), 2.5 );
1042           else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton)    ) fhTrackMatchedMCParticle[cut]->Fill(cluster->E(), 0.5 );
1043           else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron)  ) fhTrackMatchedMCParticle[cut]->Fill(cluster->E(), 1.5 );
1044           else                                                                                 fhTrackMatchedMCParticle[cut]->Fill(cluster->E(), 3.5 );
1045
1046           // Check if several particles contributed to cluster and discard overlapped mesons
1047           if(!GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) ||
1048              !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta)){
1049             if(cluster->GetNLabels()==1)
1050             {
1051               fhTrackMatchedDEtaMCNoOverlap[cut]->Fill(cluster->E(),dZ);
1052               fhTrackMatchedDPhiMCNoOverlap[cut]->Fill(cluster->E(),dR);
1053             }
1054             else
1055             {
1056               fhTrackMatchedDEtaMCOverlap[cut]->Fill(cluster->E(),dZ);
1057               fhTrackMatchedDPhiMCOverlap[cut]->Fill(cluster->E(),dR);
1058             }
1059
1060           }// Check overlaps
1061
1062         }
1063         else
1064         {
1065           if       ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0)      ||
1066                      GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta)       ) fhTrackMatchedMCParticle[cut]->Fill(cluster->E(), 6.5 );
1067           else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton)    ) fhTrackMatchedMCParticle[cut]->Fill(cluster->E(), 4.5 );
1068           else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron)  ) fhTrackMatchedMCParticle[cut]->Fill(cluster->E(), 5.5 );
1069           else                                                                                 fhTrackMatchedMCParticle[cut]->Fill(cluster->E(), 7.5 );
1070
1071           // Check if several particles contributed to cluster
1072           if(!GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) ||
1073              !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta)){
1074
1075             fhTrackMatchedDEtaMCConversion[cut]->Fill(cluster->E(),dZ);
1076             fhTrackMatchedDPhiMCConversion[cut]->Fill(cluster->E(),dR);
1077
1078           }// Check overlaps
1079
1080         }
1081
1082       } // MC
1083
1084     } // residuals window
1085
1086   } // Small residual
1087
1088 }
1089
1090 //___________________________________________
1091 TObjString *  AliAnaPhoton::GetAnalysisCuts()
1092 {
1093   //Save parameters used for analysis
1094   TString parList ; //this will be list of parameters used for this analysis.
1095   const Int_t buffersize = 255;
1096   char onePar[buffersize] ;
1097
1098   snprintf(onePar,buffersize,"--- AliAnaPhoton ---\n") ;
1099   parList+=onePar ;
1100   snprintf(onePar,buffersize,"Calorimeter: %s\n",fCalorimeter.Data()) ;
1101   parList+=onePar ;
1102   snprintf(onePar,buffersize,"fMinDist =%2.2f (Minimal distance to bad channel to accept cluster) \n",fMinDist) ;
1103   parList+=onePar ;
1104   snprintf(onePar,buffersize,"fMinDist2=%2.2f (Cuts on Minimal distance to study acceptance evaluation) \n",fMinDist2) ;
1105   parList+=onePar ;
1106   snprintf(onePar,buffersize,"fMinDist3=%2.2f (One more cut on distance used for acceptance-efficiency study) \n",fMinDist3) ;
1107   parList+=onePar ;
1108   snprintf(onePar,buffersize,"fRejectTrackMatch: %d\n",fRejectTrackMatch) ;
1109   parList+=onePar ;
1110
1111   //Get parameters set in base class.
1112   parList += GetBaseParametersList() ;
1113
1114   //Get parameters set in PID class.
1115   parList += GetCaloPID()->GetPIDParametersList() ;
1116
1117   //Get parameters set in FiducialCut class (not available yet)
1118   //parlist += GetFidCut()->GetFidCutParametersList()
1119
1120   return new TObjString(parList) ;
1121 }
1122
1123 //________________________________________________________________________
1124 TList *  AliAnaPhoton::GetCreateOutputObjects()
1125 {
1126   // Create histograms to be saved in output file and
1127   // store them in outputContainer
1128   TList * outputContainer = new TList() ;
1129   outputContainer->SetName("PhotonHistos") ;
1130
1131   Int_t nptbins  = GetHistogramRanges()->GetHistoPtBins();  Float_t ptmax  = GetHistogramRanges()->GetHistoPtMax();  Float_t ptmin  = GetHistogramRanges()->GetHistoPtMin();
1132   Int_t nphibins = GetHistogramRanges()->GetHistoPhiBins(); Float_t phimax = GetHistogramRanges()->GetHistoPhiMax(); Float_t phimin = GetHistogramRanges()->GetHistoPhiMin();
1133   Int_t netabins = GetHistogramRanges()->GetHistoEtaBins(); Float_t etamax = GetHistogramRanges()->GetHistoEtaMax(); Float_t etamin = GetHistogramRanges()->GetHistoEtaMin();
1134   Int_t ssbins   = GetHistogramRanges()->GetHistoShowerShapeBins();  Float_t ssmax   = GetHistogramRanges()->GetHistoShowerShapeMax();  Float_t ssmin   = GetHistogramRanges()->GetHistoShowerShapeMin();
1135   Int_t nbins    = GetHistogramRanges()->GetHistoNClusterCellBins(); Int_t   nmax    = GetHistogramRanges()->GetHistoNClusterCellMax(); Int_t   nmin    = GetHistogramRanges()->GetHistoNClusterCellMin();
1136   Int_t ntimebins= GetHistogramRanges()->GetHistoTimeBins();         Float_t timemax = GetHistogramRanges()->GetHistoTimeMax();         Float_t timemin = GetHistogramRanges()->GetHistoTimeMin();
1137
1138   Int_t   nresetabins = GetHistogramRanges()->GetHistoTrackResidualEtaBins();
1139   Float_t resetamax   = GetHistogramRanges()->GetHistoTrackResidualEtaMax();
1140   Float_t resetamin   = GetHistogramRanges()->GetHistoTrackResidualEtaMin();
1141   Int_t   nresphibins = GetHistogramRanges()->GetHistoTrackResidualPhiBins();
1142   Float_t resphimax   = GetHistogramRanges()->GetHistoTrackResidualPhiMax();
1143   Float_t resphimin   = GetHistogramRanges()->GetHistoTrackResidualPhiMin();
1144
1145   Int_t   ndedxbins   = GetHistogramRanges()->GetHistodEdxBins();
1146   Float_t dedxmax     = GetHistogramRanges()->GetHistodEdxMax();
1147   Float_t dedxmin     = GetHistogramRanges()->GetHistodEdxMin();
1148   Int_t   nPoverEbins = GetHistogramRanges()->GetHistoPOverEBins();
1149   Float_t pOverEmax   = GetHistogramRanges()->GetHistoPOverEMax();
1150   Float_t pOverEmin   = GetHistogramRanges()->GetHistoPOverEMin();
1151
1152   TString cut[] = {"Open","Reader","E","Time","NCells","Fidutial","Matching","Bad","PID"};
1153   for (Int_t i = 0; i < 9 ;  i++)
1154   {
1155     fhClusterCuts[i] = new TH1F(Form("hCut_%d_%s", i, cut[i].Data()),
1156                                 Form("Number of clusters that pass cuts <= %d, %s", i, cut[i].Data()),
1157                                 nptbins,ptmin,ptmax);
1158     fhClusterCuts[i]->SetYTitle("dN/dE ");
1159     fhClusterCuts[i]->SetXTitle("E (GeV)");
1160     outputContainer->Add(fhClusterCuts[i]) ;
1161   }
1162
1163   fhNCellsE  = new TH2F ("hNCellsE","# of cells in cluster vs E of clusters", nptbins,ptmin,ptmax, nbins,nmin,nmax);
1164   fhNCellsE->SetXTitle("E (GeV)");
1165   fhNCellsE->SetYTitle("# of cells in cluster");
1166   outputContainer->Add(fhNCellsE);
1167
1168   fhCellsE  = new TH2F ("hCellsE","energy of cells in cluster vs E of clusters", nptbins,ptmin,ptmax, nptbins*2,ptmin,ptmax);
1169   fhCellsE->SetXTitle("E_{cluster} (GeV)");
1170   fhCellsE->SetYTitle("E_{cell} (GeV)");
1171   outputContainer->Add(fhCellsE);
1172
1173   fhTimeE  = new TH2F ("hTimeE","time of cluster vs E of clusters", nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
1174   fhTimeE->SetXTitle("E (GeV)");
1175   fhTimeE->SetYTitle("time (ns)");
1176   outputContainer->Add(fhTimeE);
1177
1178   fhMaxCellDiffClusterE  = new TH2F ("hMaxCellDiffClusterE","energy vs difference of cluster energy - max cell energy / cluster energy, good clusters",
1179                                     nptbins,ptmin,ptmax, 500,0,1.);
1180   fhMaxCellDiffClusterE->SetXTitle("E_{cluster} (GeV) ");
1181   fhMaxCellDiffClusterE->SetYTitle("(E_{cluster} - E_{cell max})/ E_{cluster}");
1182   outputContainer->Add(fhMaxCellDiffClusterE);
1183
1184   fhEPhoton  = new TH1F("hEPhoton","Number of #gamma over calorimeter vs energy",nptbins,ptmin,ptmax);
1185   fhEPhoton->SetYTitle("N");
1186   fhEPhoton->SetXTitle("E_{#gamma}(GeV)");
1187   outputContainer->Add(fhEPhoton) ;
1188
1189   fhPtPhoton  = new TH1F("hPtPhoton","Number of #gamma over calorimeter vs p_{T}",nptbins,ptmin,ptmax);
1190   fhPtPhoton->SetYTitle("N");
1191   fhPtPhoton->SetXTitle("p_{T #gamma}(GeV/c)");
1192   outputContainer->Add(fhPtPhoton) ;
1193
1194   fhPhiPhoton  = new TH2F
1195     ("hPhiPhoton","#phi_{#gamma} vs p_{T}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
1196   fhPhiPhoton->SetYTitle("#phi (rad)");
1197   fhPhiPhoton->SetXTitle("p_{T #gamma} (GeV/c)");
1198   outputContainer->Add(fhPhiPhoton) ;
1199
1200   fhEtaPhoton  = new TH2F
1201     ("hEtaPhoton","#eta_{#gamma} vs p_{T}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
1202   fhEtaPhoton->SetYTitle("#eta");
1203   fhEtaPhoton->SetXTitle("p_{T #gamma} (GeV/c)");
1204   outputContainer->Add(fhEtaPhoton) ;
1205
1206   fhEtaPhiPhoton  = new TH2F
1207   ("hEtaPhiPhoton","#eta vs #phi",netabins,etamin,etamax,nphibins,phimin,phimax);
1208   fhEtaPhiPhoton->SetYTitle("#phi (rad)");
1209   fhEtaPhiPhoton->SetXTitle("#eta");
1210   outputContainer->Add(fhEtaPhiPhoton) ;
1211   if(GetMinPt() < 0.5){
1212     fhEtaPhi05Photon  = new TH2F
1213     ("hEtaPhi05Photon","#eta vs #phi, E > 0.5",netabins,etamin,etamax,nphibins,phimin,phimax);
1214     fhEtaPhi05Photon->SetYTitle("#phi (rad)");
1215     fhEtaPhi05Photon->SetXTitle("#eta");
1216     outputContainer->Add(fhEtaPhi05Photon) ;
1217   }
1218
1219   //Shower shape
1220   if(fFillSSHistograms){
1221
1222     fhLam0E  = new TH2F ("hLam0E","#lambda_{0}^{2} vs E", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1223     fhLam0E->SetYTitle("#lambda_{0}^{2}");
1224     fhLam0E->SetXTitle("E (GeV)");
1225     outputContainer->Add(fhLam0E);
1226
1227     fhLam1E  = new TH2F ("hLam1E","#lambda_{1}^{2} vs E", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1228     fhLam1E->SetYTitle("#lambda_{1}^{2}");
1229     fhLam1E->SetXTitle("E (GeV)");
1230     outputContainer->Add(fhLam1E);
1231
1232     fhDispE  = new TH2F ("hDispE"," dispersion^{2} vs E", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1233     fhDispE->SetYTitle("D^{2}");
1234     fhDispE->SetXTitle("E (GeV) ");
1235     outputContainer->Add(fhDispE);
1236
1237     if(!fRejectTrackMatch)
1238     {
1239       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);
1240       fhLam0ETM->SetYTitle("#lambda_{0}^{2}");
1241       fhLam0ETM->SetXTitle("E (GeV)");
1242       outputContainer->Add(fhLam0ETM);
1243
1244       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);
1245       fhLam1ETM->SetYTitle("#lambda_{1}^{2}");
1246       fhLam1ETM->SetXTitle("E (GeV)");
1247       outputContainer->Add(fhLam1ETM);
1248
1249       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);
1250       fhDispETM->SetYTitle("D^{2}");
1251       fhDispETM->SetXTitle("E (GeV) ");
1252       outputContainer->Add(fhDispETM);
1253     }
1254
1255     if(fCalorimeter == "EMCAL")
1256     {
1257       fhLam0ETRD  = new TH2F ("hLam0ETRD","#lambda_{0}^{2} vs E, EMCAL SM covered by TRD", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1258       fhLam0ETRD->SetYTitle("#lambda_{0}^{2}");
1259       fhLam0ETRD->SetXTitle("E (GeV)");
1260       outputContainer->Add(fhLam0ETRD);
1261
1262       fhLam1ETRD  = new TH2F ("hLam1ETRD","#lambda_{1}^{2} vs E, EMCAL SM covered by TRD", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1263       fhLam1ETRD->SetYTitle("#lambda_{1}^{2}");
1264       fhLam1ETRD->SetXTitle("E (GeV)");
1265       outputContainer->Add(fhLam1ETRD);
1266
1267       fhDispETRD  = new TH2F ("hDispETRD"," dispersion^{2} vs E, EMCAL SM covered by TRD", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1268       fhDispETRD->SetYTitle("Dispersion^{2}");
1269       fhDispETRD->SetXTitle("E (GeV) ");
1270       outputContainer->Add(fhDispETRD);
1271
1272       if(!fRejectTrackMatch)
1273       {
1274         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);
1275         fhLam0ETMTRD->SetYTitle("#lambda_{0}^{2}");
1276         fhLam0ETMTRD->SetXTitle("E (GeV)");
1277         outputContainer->Add(fhLam0ETMTRD);
1278
1279         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);
1280         fhLam1ETMTRD->SetYTitle("#lambda_{1}^{2}");
1281         fhLam1ETMTRD->SetXTitle("E (GeV)");
1282         outputContainer->Add(fhLam1ETMTRD);
1283
1284         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);
1285         fhDispETMTRD->SetYTitle("Dispersion^{2}");
1286         fhDispETMTRD->SetXTitle("E (GeV) ");
1287         outputContainer->Add(fhDispETMTRD);
1288       }
1289     }
1290
1291     fhNCellsLam0LowE  = new TH2F ("hNCellsLam0LowE","N_{cells} in cluster vs #lambda_{0}^{2}, E < 2 GeV", nbins,nmin, nmax, ssbins,ssmin,ssmax);
1292     fhNCellsLam0LowE->SetXTitle("N_{cells}");
1293     fhNCellsLam0LowE->SetYTitle("#lambda_{0}^{2}");
1294     outputContainer->Add(fhNCellsLam0LowE);
1295
1296     fhNCellsLam0HighE  = new TH2F ("hNCellsLam0HighE","N_{cells} in cluster vs #lambda_{0}^{2}, E > 2 GeV", nbins,nmin, nmax, ssbins,ssmin,ssmax);
1297     fhNCellsLam0HighE->SetXTitle("N_{cells}");
1298     fhNCellsLam0HighE->SetYTitle("#lambda_{0}^{2}");
1299     outputContainer->Add(fhNCellsLam0HighE);
1300
1301     fhNCellsLam1LowE  = new TH2F ("hNCellsLam1LowE","N_{cells} in cluster vs #lambda_{1}^{2}, E < 2 GeV", nbins,nmin, nmax, ssbins,ssmin,ssmax);
1302     fhNCellsLam1LowE->SetXTitle("N_{cells}");
1303     fhNCellsLam1LowE->SetYTitle("#lambda_{0}^{2}");
1304     outputContainer->Add(fhNCellsLam1LowE);
1305
1306     fhNCellsLam1HighE  = new TH2F ("hNCellsLam1HighE","N_{cells} in cluster vs #lambda_{1}^{2}, E > 2 GeV", nbins,nmin, nmax, ssbins,ssmin,ssmax);
1307     fhNCellsLam1HighE->SetXTitle("N_{cells}");
1308     fhNCellsLam1HighE->SetYTitle("#lambda_{0}^{2}");
1309     outputContainer->Add(fhNCellsLam1HighE);
1310
1311     fhNCellsDispLowE  = new TH2F ("hNCellsDispLowE","N_{cells} in cluster vs dispersion^{2}, E < 2 GeV", nbins,nmin, nmax, ssbins,ssmin,ssmax);
1312     fhNCellsDispLowE->SetXTitle("N_{cells}");
1313     fhNCellsDispLowE->SetYTitle("D^{2}");
1314     outputContainer->Add(fhNCellsDispLowE);
1315
1316     fhNCellsDispHighE  = new TH2F ("hNCellsDispHighE","N_{cells} in cluster vs dispersion^{2}, E < 2 GeV", nbins,nmin, nmax, ssbins,ssmin,ssmax);
1317     fhNCellsDispHighE->SetXTitle("N_{cells}");
1318     fhNCellsDispHighE->SetYTitle("D^{2}");
1319     outputContainer->Add(fhNCellsDispHighE);
1320
1321     fhEtaLam0LowE  = new TH2F ("hEtaLam0LowE","#eta vs #lambda_{0}^{2}, E < 2 GeV", netabins,etamin,etamax, ssbins,ssmin,ssmax);
1322     fhEtaLam0LowE->SetYTitle("#lambda_{0}^{2}");
1323     fhEtaLam0LowE->SetXTitle("#eta");
1324     outputContainer->Add(fhEtaLam0LowE);
1325
1326     fhPhiLam0LowE  = new TH2F ("hPhiLam0LowE","#phi vs #lambda_{0}^{2}, E < 2 GeV", nphibins,phimin,phimax, ssbins,ssmin,ssmax);
1327     fhPhiLam0LowE->SetYTitle("#lambda_{0}^{2}");
1328     fhPhiLam0LowE->SetXTitle("#phi");
1329     outputContainer->Add(fhPhiLam0LowE);
1330
1331     fhEtaLam0HighE  = new TH2F ("hEtaLam0HighE","#eta vs #lambda_{0}^{2}, E > 2 GeV", netabins,etamin,etamax, ssbins,ssmin,ssmax);
1332     fhEtaLam0HighE->SetYTitle("#lambda_{0}^{2}");
1333     fhEtaLam0HighE->SetXTitle("#eta");
1334     outputContainer->Add(fhEtaLam0HighE);
1335
1336     fhPhiLam0HighE  = new TH2F ("hPhiLam0HighE","#phi vs #lambda_{0}^{2}, E > 2 GeV", nphibins,phimin,phimax, ssbins,ssmin,ssmax);
1337     fhPhiLam0HighE->SetYTitle("#lambda_{0}^{2}");
1338     fhPhiLam0HighE->SetXTitle("#phi");
1339     outputContainer->Add(fhPhiLam0HighE);
1340
1341     fhLam1Lam0LowE  = new TH2F ("hLam1Lam0LowE","#lambda_{0}^{2} vs #lambda_{1}^{2} in cluster of E < 2 GeV",  ssbins,ssmin,ssmax, ssbins,ssmin,ssmax);
1342     fhLam1Lam0LowE->SetYTitle("#lambda_{0}^{2}");
1343     fhLam1Lam0LowE->SetXTitle("#lambda_{1}^{2}");
1344     outputContainer->Add(fhLam1Lam0LowE);
1345
1346     fhLam1Lam0HighE  = new TH2F ("hLam1Lam0HighE","#lambda_{0}^{2} vs #lambda_{1}^{2} in cluster of E > 2 GeV",  ssbins,ssmin,ssmax, ssbins,ssmin,ssmax);
1347     fhLam1Lam0HighE->SetYTitle("#lambda_{0}^{2}");
1348     fhLam1Lam0HighE->SetXTitle("#lambda_{1}^{2}");
1349     outputContainer->Add(fhLam1Lam0HighE);
1350
1351     fhLam0DispLowE  = new TH2F ("hLam0DispLowE","#lambda_{0}^{2} vs dispersion^{2} in cluster of E < 2 GeV",  ssbins,ssmin,ssmax, ssbins,ssmin,ssmax);
1352     fhLam0DispLowE->SetXTitle("#lambda_{0}^{2}");
1353     fhLam0DispLowE->SetYTitle("D^{2}");
1354     outputContainer->Add(fhLam0DispLowE);
1355
1356     fhLam0DispHighE  = new TH2F ("hLam0DispHighE","#lambda_{0}^{2} vs dispersion^{2} in cluster of E > 2 GeV",  ssbins,ssmin,ssmax, ssbins,ssmin,ssmax);
1357     fhLam0DispHighE->SetXTitle("#lambda_{0}^{2}");
1358     fhLam0DispHighE->SetYTitle("D^{2}");
1359     outputContainer->Add(fhLam0DispHighE);
1360
1361     fhDispLam1LowE  = new TH2F ("hDispLam1LowE","Dispersion^{2} vs #lambda_{1}^{2} in cluster of E < 2 GeV",  ssbins,ssmin,ssmax, ssbins,ssmin,ssmax);
1362     fhDispLam1LowE->SetXTitle("D^{2}");
1363     fhDispLam1LowE->SetYTitle("#lambda_{1}^{2}");
1364     outputContainer->Add(fhDispLam1LowE);
1365
1366     fhDispLam1HighE  = new TH2F ("hDispLam1HighE","Dispersion^{2} vs #lambda_{1^{2}} in cluster of E > 2 GeV",  ssbins,ssmin,ssmax, ssbins,ssmin,ssmax);
1367     fhDispLam1HighE->SetXTitle("D^{2}");
1368     fhDispLam1HighE->SetYTitle("#lambda_{1}^{2}");
1369     outputContainer->Add(fhDispLam1HighE);
1370
1371   } // Shower shape
1372
1373   // Track Matching
1374
1375   if(fFillTMHisto)
1376   {
1377
1378     fhTrackMatchedDEta[0]  = new TH2F
1379     ("hTrackMatchedDEtaNoCut",
1380      "d#eta of cluster-track vs cluster energy, no photon cuts",
1381      nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
1382     fhTrackMatchedDEta[0]->SetYTitle("d#eta");
1383     fhTrackMatchedDEta[0]->SetXTitle("E_{cluster} (GeV)");
1384
1385     fhTrackMatchedDPhi[0]  = new TH2F
1386     ("hTrackMatchedDPhiNoCut",
1387      "d#phi of cluster-track vs cluster energy, no photon cuts",
1388      nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
1389     fhTrackMatchedDPhi[0]->SetYTitle("d#phi (rad)");
1390     fhTrackMatchedDPhi[0]->SetXTitle("E_{cluster} (GeV)");
1391
1392     fhTrackMatchedDEtaDPhi[0]  = new TH2F
1393     ("hTrackMatchedDEtaDPhiNoCut",
1394      "d#eta vs d#phi of cluster-track vs cluster energy, no photon cuts",
1395      nresetabins,resetamin,resetamax,nresphibins,resphimin,resphimax);
1396     fhTrackMatchedDEtaDPhi[0]->SetYTitle("d#phi (rad)");
1397     fhTrackMatchedDEtaDPhi[0]->SetXTitle("d#eta");
1398
1399     fhdEdx[0]  = new TH2F ("hdEdxNoCut","matched track <dE/dx> vs cluster E, no photon cuts ",
1400                            nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
1401     fhdEdx[0]->SetXTitle("E (GeV)");
1402     fhdEdx[0]->SetYTitle("<dE/dx>");
1403
1404     fhEOverP[0]  = new TH2F ("hEOverPNoCut","matched track E/p vs cluster E, no photon cuts ",
1405                              nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
1406     fhEOverP[0]->SetXTitle("E (GeV)");
1407     fhEOverP[0]->SetYTitle("E/p");
1408
1409     outputContainer->Add(fhTrackMatchedDEta[0]) ;
1410     outputContainer->Add(fhTrackMatchedDPhi[0]) ;
1411     outputContainer->Add(fhTrackMatchedDEtaDPhi[0]) ;
1412     outputContainer->Add(fhdEdx[0]);
1413     outputContainer->Add(fhEOverP[0]);
1414
1415     fhTrackMatchedDEta[1]  = new TH2F
1416     ("hTrackMatchedDEta",
1417      "d#eta of cluster-track vs cluster energy, no photon cuts",
1418      nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
1419     fhTrackMatchedDEta[1]->SetYTitle("d#eta");
1420     fhTrackMatchedDEta[1]->SetXTitle("E_{cluster} (GeV)");
1421
1422     fhTrackMatchedDPhi[1]  = new TH2F
1423     ("hTrackMatchedDPhi",
1424      "d#phi of cluster-track vs cluster energy, no photon cuts",
1425      nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
1426     fhTrackMatchedDPhi[1]->SetYTitle("d#phi (rad)");
1427     fhTrackMatchedDPhi[1]->SetXTitle("E_{cluster} (GeV)");
1428
1429     fhTrackMatchedDEtaDPhi[1]  = new TH2F
1430     ("hTrackMatchedDEtaDPhi",
1431      "d#eta vs d#phi of cluster-track vs cluster energy, no photon cuts",
1432      nresetabins,resetamin,resetamax,nresphibins,resphimin,resphimax);
1433     fhTrackMatchedDEtaDPhi[1]->SetYTitle("d#phi (rad)");
1434     fhTrackMatchedDEtaDPhi[1]->SetXTitle("d#eta");
1435
1436     fhdEdx[1]  = new TH2F ("hdEdx","matched track <dE/dx> vs cluster E ",
1437                            nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
1438     fhdEdx[1]->SetXTitle("E (GeV)");
1439     fhdEdx[1]->SetYTitle("<dE/dx>");
1440
1441     fhEOverP[1]  = new TH2F ("hEOverP","matched track E/p vs cluster E ",
1442                              nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
1443     fhEOverP[1]->SetXTitle("E (GeV)");
1444     fhEOverP[1]->SetYTitle("E/p");
1445
1446     outputContainer->Add(fhTrackMatchedDEta[1]) ;
1447     outputContainer->Add(fhTrackMatchedDPhi[1]) ;
1448     outputContainer->Add(fhTrackMatchedDEtaDPhi[1]) ;
1449     outputContainer->Add(fhdEdx[1]);
1450     outputContainer->Add(fhEOverP[1]);
1451
1452     if(fCalorimeter=="EMCAL")
1453     {
1454
1455       fhTrackMatchedDEtaTRD[0]  = new TH2F
1456       ("hTrackMatchedDEtaTRDNoCut",
1457        "d#eta of cluster-track vs cluster energy, SM behind TRD, no photon cuts",
1458        nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
1459       fhTrackMatchedDEtaTRD[0]->SetYTitle("d#eta");
1460       fhTrackMatchedDEtaTRD[0]->SetXTitle("E_{cluster} (GeV)");
1461
1462       fhTrackMatchedDPhiTRD[0]  = new TH2F
1463       ("hTrackMatchedDPhiTRDNoCut",
1464        "d#phi of cluster-track vs cluster energy, SM behing TRD, no photon cuts",
1465        nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
1466       fhTrackMatchedDPhiTRD[0]->SetYTitle("d#phi (rad)");
1467       fhTrackMatchedDPhiTRD[0]->SetXTitle("E_{cluster} (GeV)");
1468
1469       fhEOverPTRD[0]  = new TH2F ("hEOverPTRDNoCut","matched track E/p vs cluster E, behind TRD, no photon cuts ",
1470                                   nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
1471       fhEOverPTRD[0]->SetXTitle("E (GeV)");
1472       fhEOverPTRD[0]->SetYTitle("E/p");
1473
1474       outputContainer->Add(fhTrackMatchedDEtaTRD[0]) ;
1475       outputContainer->Add(fhTrackMatchedDPhiTRD[0]) ;
1476       outputContainer->Add(fhEOverPTRD[0]);
1477
1478       fhTrackMatchedDEtaTRD[1]  = new TH2F
1479       ("hTrackMatchedDEtaTRD",
1480        "d#eta of cluster-track vs cluster energy, SM behind TRD",
1481        nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
1482       fhTrackMatchedDEtaTRD[1]->SetYTitle("d#eta");
1483       fhTrackMatchedDEtaTRD[1]->SetXTitle("E_{cluster} (GeV)");
1484
1485       fhTrackMatchedDPhiTRD[1]  = new TH2F
1486       ("hTrackMatchedDPhiTRD",
1487        "d#phi of cluster-track vs cluster energy, SM behing TRD",
1488        nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
1489       fhTrackMatchedDPhiTRD[1]->SetYTitle("d#phi (rad)");
1490       fhTrackMatchedDPhiTRD[1]->SetXTitle("E_{cluster} (GeV)");
1491
1492       fhEOverPTRD[1]  = new TH2F ("hEOverPTRD","matched track E/p vs cluster E, behind TRD ",
1493                                   nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
1494       fhEOverPTRD[1]->SetXTitle("E (GeV)");
1495       fhEOverPTRD[1]->SetYTitle("E/p");
1496
1497       outputContainer->Add(fhTrackMatchedDEtaTRD[1]) ;
1498       outputContainer->Add(fhTrackMatchedDPhiTRD[1]) ;
1499       outputContainer->Add(fhEOverPTRD[1]);
1500
1501     }
1502
1503     if(IsDataMC())
1504     {
1505
1506       fhTrackMatchedDEtaMCNoOverlap[0]  = new TH2F
1507       ("hTrackMatchedDEtaMCNoOverlapNoCut",
1508        "d#eta of cluster-track vs cluster energy, no other MC particles overlap",
1509        nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
1510       fhTrackMatchedDEtaMCNoOverlap[0]->SetYTitle("d#eta");
1511       fhTrackMatchedDEtaMCNoOverlap[0]->SetXTitle("E_{cluster} (GeV)");
1512
1513       fhTrackMatchedDPhiMCNoOverlap[0]  = new TH2F
1514       ("hTrackMatchedDPhiMCNoOverlapNoCut",
1515        "d#phi of cluster-track vs cluster energy, no other MC particles overlap",
1516        nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
1517       fhTrackMatchedDPhiMCNoOverlap[0]->SetYTitle("d#phi (rad)");
1518       fhTrackMatchedDPhiMCNoOverlap[0]->SetXTitle("E_{cluster} (GeV)");
1519
1520       outputContainer->Add(fhTrackMatchedDEtaMCNoOverlap[0]) ;
1521       outputContainer->Add(fhTrackMatchedDPhiMCNoOverlap[0]) ;
1522
1523       fhTrackMatchedDEtaMCNoOverlap[1]  = new TH2F
1524       ("hTrackMatchedDEtaMCNoOverlap",
1525        "d#eta of cluster-track vs cluster energy, no other MC particles overlap",
1526        nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
1527       fhTrackMatchedDEtaMCNoOverlap[1]->SetYTitle("d#eta");
1528       fhTrackMatchedDEtaMCNoOverlap[1]->SetXTitle("E_{cluster} (GeV)");
1529
1530       fhTrackMatchedDPhiMCNoOverlap[1]  = new TH2F
1531       ("hTrackMatchedDPhiMCNoOverlap",
1532        "d#phi of cluster-track vs cluster energy, no other MC particles overlap",
1533        nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
1534       fhTrackMatchedDPhiMCNoOverlap[1]->SetYTitle("d#phi (rad)");
1535       fhTrackMatchedDPhiMCNoOverlap[1]->SetXTitle("E_{cluster} (GeV)");
1536
1537       outputContainer->Add(fhTrackMatchedDEtaMCNoOverlap[1]) ;
1538       outputContainer->Add(fhTrackMatchedDPhiMCNoOverlap[1]) ;
1539
1540       fhTrackMatchedDEtaMCOverlap[0]  = new TH2F
1541       ("hTrackMatchedDEtaMCOverlapNoCut",
1542        "d#eta of cluster-track vs cluster energy, several MC particles overlap",
1543        nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
1544       fhTrackMatchedDEtaMCOverlap[0]->SetYTitle("d#eta");
1545       fhTrackMatchedDEtaMCOverlap[0]->SetXTitle("E_{cluster} (GeV)");
1546
1547       fhTrackMatchedDPhiMCOverlap[0]  = new TH2F
1548       ("hTrackMatchedDPhiMCOverlapNoCut",
1549        "d#phi of cluster-track vs cluster energy, several MC particles overlap",
1550        nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
1551       fhTrackMatchedDPhiMCOverlap[0]->SetYTitle("d#phi (rad)");
1552       fhTrackMatchedDPhiMCOverlap[0]->SetXTitle("E_{cluster} (GeV)");
1553
1554       outputContainer->Add(fhTrackMatchedDEtaMCOverlap[0]) ;
1555       outputContainer->Add(fhTrackMatchedDPhiMCOverlap[0]) ;
1556
1557       fhTrackMatchedDEtaMCOverlap[1]  = new TH2F
1558       ("hTrackMatchedDEtaMCOverlap",
1559        "d#eta of cluster-track vs cluster energy, several MC particles overlap",
1560        nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
1561       fhTrackMatchedDEtaMCOverlap[1]->SetYTitle("d#eta");
1562       fhTrackMatchedDEtaMCOverlap[1]->SetXTitle("E_{cluster} (GeV)");
1563
1564       fhTrackMatchedDPhiMCOverlap[1]  = new TH2F
1565       ("hTrackMatchedDPhiMCOverlap",
1566        "d#phi of cluster-track vs cluster energy, several MC particles overlap",
1567        nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
1568       fhTrackMatchedDPhiMCOverlap[1]->SetYTitle("d#phi (rad)");
1569       fhTrackMatchedDPhiMCOverlap[1]->SetXTitle("E_{cluster} (GeV)");
1570
1571       outputContainer->Add(fhTrackMatchedDEtaMCOverlap[1]) ;
1572       outputContainer->Add(fhTrackMatchedDPhiMCOverlap[1]) ;
1573
1574       fhTrackMatchedDEtaMCConversion[0]  = new TH2F
1575       ("hTrackMatchedDEtaMCConversionNoCut",
1576        "d#eta of cluster-track vs cluster energy, no other MC particles overlap appart from conversions",
1577        nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
1578       fhTrackMatchedDEtaMCConversion[0]->SetYTitle("d#eta");
1579       fhTrackMatchedDEtaMCConversion[0]->SetXTitle("E_{cluster} (GeV)");
1580
1581       fhTrackMatchedDPhiMCConversion[0]  = new TH2F
1582       ("hTrackMatchedDPhiMCConversionNoCut",
1583        "d#phi of cluster-track vs cluster energy, no other MC particles overlap appart from conversions",
1584        nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
1585       fhTrackMatchedDPhiMCConversion[0]->SetYTitle("d#phi (rad)");
1586       fhTrackMatchedDPhiMCConversion[0]->SetXTitle("E_{cluster} (GeV)");
1587
1588       outputContainer->Add(fhTrackMatchedDEtaMCConversion[0]) ;
1589       outputContainer->Add(fhTrackMatchedDPhiMCConversion[0]) ;
1590
1591
1592       fhTrackMatchedDEtaMCConversion[1]  = new TH2F
1593       ("hTrackMatchedDEtaMCConversion",
1594        "d#eta of cluster-track vs cluster energy, no other MC particles overlap appart from conversions",
1595        nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
1596       fhTrackMatchedDEtaMCConversion[1]->SetYTitle("d#eta");
1597       fhTrackMatchedDEtaMCConversion[1]->SetXTitle("E_{cluster} (GeV)");
1598
1599       fhTrackMatchedDPhiMCConversion[1]  = new TH2F
1600       ("hTrackMatchedDPhiMCConversion",
1601        "d#phi of cluster-track vs cluster energy, no other MC particles overlap appart from conversions",
1602        nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
1603       fhTrackMatchedDPhiMCConversion[1]->SetYTitle("d#phi (rad)");
1604       fhTrackMatchedDPhiMCConversion[1]->SetXTitle("E_{cluster} (GeV)");
1605
1606       outputContainer->Add(fhTrackMatchedDEtaMCConversion[1]) ;
1607       outputContainer->Add(fhTrackMatchedDPhiMCConversion[1]) ;
1608
1609
1610       fhTrackMatchedMCParticle[0]  = new TH2F
1611       ("hTrackMatchedMCParticleNoCut",
1612        "Origin of particle vs energy",
1613        nptbins,ptmin,ptmax,8,0,8);
1614       fhTrackMatchedMCParticle[0]->SetXTitle("E (GeV)");
1615       //fhTrackMatchedMCParticle[0]->SetYTitle("Particle type");
1616
1617       fhTrackMatchedMCParticle[0]->GetYaxis()->SetBinLabel(1 ,"Photon");
1618       fhTrackMatchedMCParticle[0]->GetYaxis()->SetBinLabel(2 ,"Electron");
1619       fhTrackMatchedMCParticle[0]->GetYaxis()->SetBinLabel(3 ,"Meson Merged");
1620       fhTrackMatchedMCParticle[0]->GetYaxis()->SetBinLabel(4 ,"Rest");
1621       fhTrackMatchedMCParticle[0]->GetYaxis()->SetBinLabel(5 ,"Conv. Photon");
1622       fhTrackMatchedMCParticle[0]->GetYaxis()->SetBinLabel(6 ,"Conv. Electron");
1623       fhTrackMatchedMCParticle[0]->GetYaxis()->SetBinLabel(7 ,"Conv. Merged");
1624       fhTrackMatchedMCParticle[0]->GetYaxis()->SetBinLabel(8 ,"Conv. Rest");
1625
1626       fhTrackMatchedMCParticle[1]  = new TH2F
1627       ("hTrackMatchedMCParticle",
1628        "Origin of particle vs energy",
1629        nptbins,ptmin,ptmax,8,0,8);
1630       fhTrackMatchedMCParticle[1]->SetXTitle("E (GeV)");
1631       //fhTrackMatchedMCParticle[1]->SetYTitle("Particle type");
1632
1633       fhTrackMatchedMCParticle[1]->GetYaxis()->SetBinLabel(1 ,"Photon");
1634       fhTrackMatchedMCParticle[1]->GetYaxis()->SetBinLabel(2 ,"Electron");
1635       fhTrackMatchedMCParticle[1]->GetYaxis()->SetBinLabel(3 ,"Meson Merged");
1636       fhTrackMatchedMCParticle[1]->GetYaxis()->SetBinLabel(4 ,"Rest");
1637       fhTrackMatchedMCParticle[1]->GetYaxis()->SetBinLabel(5 ,"Conv. Photon");
1638       fhTrackMatchedMCParticle[1]->GetYaxis()->SetBinLabel(6 ,"Conv. Electron");
1639       fhTrackMatchedMCParticle[1]->GetYaxis()->SetBinLabel(7 ,"Conv. Merged");
1640       fhTrackMatchedMCParticle[1]->GetYaxis()->SetBinLabel(8 ,"Conv. Rest");
1641
1642       outputContainer->Add(fhTrackMatchedMCParticle[0]);
1643       outputContainer->Add(fhTrackMatchedMCParticle[1]);
1644
1645     }
1646   }
1647
1648
1649   if(IsDataMC()){
1650
1651     TString ptype[] = { "#gamma", "#gamma_{#pi decay}","#gamma_{other decay}", "#pi^{0}","#eta",
1652                         "e^{#pm}","#gamma->e^{#pm}","hadron?","Anti-N","Anti-P",
1653                         "#gamma_{prompt}","#gamma_{fragmentation}","#gamma_{ISR}","String"                                } ;
1654
1655     TString pname[] = { "Photon","PhotonPi0Decay","PhotonOtherDecay","Pi0","Eta","Electron",
1656                         "Conversion", "Hadron", "AntiNeutron","AntiProton",
1657                         "PhotonPrompt","PhotonFragmentation","PhotonISR","String" } ;
1658
1659     for(Int_t i = 0; i < fNOriginHistograms; i++){
1660
1661       fhMCE[i]  = new TH1F(Form("hE_MC%s",pname[i].Data()),
1662                                 Form("cluster from %s : E ",ptype[i].Data()),
1663                                 nptbins,ptmin,ptmax);
1664       fhMCE[i]->SetXTitle("E (GeV)");
1665       outputContainer->Add(fhMCE[i]) ;
1666
1667       fhMCPt[i]  = new TH1F(Form("hPt_MC%s",pname[i].Data()),
1668                            Form("cluster from %s : p_{T} ",ptype[i].Data()),
1669                            nptbins,ptmin,ptmax);
1670       fhMCPt[i]->SetXTitle("p_{T} (GeV/c)");
1671       outputContainer->Add(fhMCPt[i]) ;
1672
1673       fhMCEta[i]  = new TH2F(Form("hEta_MC%s",pname[i].Data()),
1674                            Form("cluster from %s : #eta ",ptype[i].Data()),
1675                            nptbins,ptmin,ptmax,netabins,etamin,etamax);
1676       fhMCEta[i]->SetYTitle("#eta");
1677       fhMCEta[i]->SetXTitle("E (GeV)");
1678       outputContainer->Add(fhMCEta[i]) ;
1679
1680       fhMCPhi[i]  = new TH2F(Form("hPhi_MC%s",pname[i].Data()),
1681                            Form("cluster from %s : #phi ",ptype[i].Data()),
1682                            nptbins,ptmin,ptmax,nphibins,phimin,phimax);
1683       fhMCPhi[i]->SetYTitle("#phi (rad)");
1684       fhMCPhi[i]->SetXTitle("E (GeV)");
1685       outputContainer->Add(fhMCPhi[i]) ;
1686
1687
1688       fhMCDeltaE[i]  = new TH2F (Form("hDeltaE_MC%s",pname[i].Data()),
1689                                  Form("MC - Reco E from %s",pname[i].Data()),
1690                                  nptbins,ptmin,ptmax, 200,-50,50);
1691       fhMCDeltaE[i]->SetXTitle("#Delta E (GeV)");
1692       outputContainer->Add(fhMCDeltaE[i]);
1693
1694       fhMCDeltaPt[i]  = new TH2F (Form("hDeltaPt_MC%s",pname[i].Data()),
1695                                   Form("MC - Reco p_{T} from %s",pname[i].Data()),
1696                                   nptbins,ptmin,ptmax, 200,-50,50);
1697       fhMCDeltaPt[i]->SetXTitle("#Delta p_{T} (GeV/c)");
1698       outputContainer->Add(fhMCDeltaPt[i]);
1699
1700       fhMC2E[i]  = new TH2F (Form("h2E_MC%s",pname[i].Data()),
1701                              Form("E distribution, reconstructed vs generated from %s",pname[i].Data()),
1702                              nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
1703       fhMC2E[i]->SetXTitle("E_{rec} (GeV)");
1704       fhMC2E[i]->SetYTitle("E_{gen} (GeV)");
1705       outputContainer->Add(fhMC2E[i]);
1706
1707       fhMC2Pt[i]  = new TH2F (Form("h2Pt_MC%s",pname[i].Data()),
1708                               Form("p_T distribution, reconstructed vs generated from %s",pname[i].Data()),
1709                               nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
1710       fhMC2Pt[i]->SetXTitle("p_{T,rec} (GeV/c)");
1711       fhMC2Pt[i]->SetYTitle("p_{T,gen} (GeV/c)");
1712       outputContainer->Add(fhMC2Pt[i]);
1713
1714
1715     }
1716
1717     TString pptype[] = { "#gamma", "#gamma_{#pi decay}","#gamma_{other decay}","hadron?",
1718                          "#gamma_{prompt}","#gamma_{fragmentation}","#gamma_{ISR}"} ;
1719
1720     TString ppname[] = { "Photon","PhotonPi0Decay","PhotonOtherDecay","Hadron",
1721                          "PhotonPrompt","PhotonFragmentation","PhotonISR"} ;
1722
1723     for(Int_t i = 0; i < fNPrimaryHistograms; i++){
1724       fhEPrimMC[i]  = new TH1F(Form("hEPrim_MC%s",ppname[i].Data()),
1725                            Form("primary photon %s : E ",pptype[i].Data()),
1726                            nptbins,ptmin,ptmax);
1727       fhEPrimMC[i]->SetXTitle("E (GeV)");
1728       outputContainer->Add(fhEPrimMC[i]) ;
1729
1730       fhPtPrimMC[i]  = new TH1F(Form("hPtPrim_MC%s",ppname[i].Data()),
1731                             Form("primary photon %s : p_{T} ",pptype[i].Data()),
1732                             nptbins,ptmin,ptmax);
1733       fhPtPrimMC[i]->SetXTitle("p_{T} (GeV/c)");
1734       outputContainer->Add(fhPtPrimMC[i]) ;
1735
1736       fhYPrimMC[i]  = new TH2F(Form("hYPrim_MC%s",ppname[i].Data()),
1737                              Form("primary photon %s : Rapidity ",pptype[i].Data()),
1738                              nptbins,ptmin,ptmax,800,-8,8);
1739       fhYPrimMC[i]->SetYTitle("Rapidity");
1740       fhYPrimMC[i]->SetXTitle("E (GeV)");
1741       outputContainer->Add(fhYPrimMC[i]) ;
1742
1743       fhPhiPrimMC[i]  = new TH2F(Form("hPhiPrim_MC%s",ppname[i].Data()),
1744                              Form("primary photon %s : #phi ",pptype[i].Data()),
1745                              nptbins,ptmin,ptmax,nphibins,phimin,phimax);
1746       fhPhiPrimMC[i]->SetYTitle("#phi (rad)");
1747       fhPhiPrimMC[i]->SetXTitle("E (GeV)");
1748       outputContainer->Add(fhPhiPrimMC[i]) ;
1749
1750
1751       fhEPrimMCAcc[i]  = new TH1F(Form("hEPrimAcc_MC%s",ppname[i].Data()),
1752                                Form("primary photon %s in acceptance: E ",pptype[i].Data()),
1753                                nptbins,ptmin,ptmax);
1754       fhEPrimMCAcc[i]->SetXTitle("E (GeV)");
1755       outputContainer->Add(fhEPrimMCAcc[i]) ;
1756
1757       fhPtPrimMCAcc[i]  = new TH1F(Form("hPtPrimAcc_MC%s",ppname[i].Data()),
1758                                 Form("primary photon %s in acceptance: p_{T} ",pptype[i].Data()),
1759                                 nptbins,ptmin,ptmax);
1760       fhPtPrimMCAcc[i]->SetXTitle("p_{T} (GeV/c)");
1761       outputContainer->Add(fhPtPrimMCAcc[i]) ;
1762
1763       fhYPrimMCAcc[i]  = new TH2F(Form("hYPrimAcc_MC%s",ppname[i].Data()),
1764                                  Form("primary photon %s in acceptance: Rapidity ",pptype[i].Data()),
1765                                  nptbins,ptmin,ptmax,100,-1,1);
1766       fhYPrimMCAcc[i]->SetYTitle("Rapidity");
1767       fhYPrimMCAcc[i]->SetXTitle("E (GeV)");
1768       outputContainer->Add(fhYPrimMCAcc[i]) ;
1769
1770       fhPhiPrimMCAcc[i]  = new TH2F(Form("hPhiPrimAcc_MC%s",ppname[i].Data()),
1771                                  Form("primary photon %s in acceptance: #phi ",pptype[i].Data()),
1772                                  nptbins,ptmin,ptmax,nphibins,phimin,phimax);
1773       fhPhiPrimMCAcc[i]->SetYTitle("#phi (rad)");
1774       fhPhiPrimMCAcc[i]->SetXTitle("E (GeV)");
1775       outputContainer->Add(fhPhiPrimMCAcc[i]) ;
1776
1777     }
1778
1779     if(fFillSSHistograms){
1780
1781       TString ptypess[] = { "#gamma","hadron?","#pi^{0}","#eta","#gamma->e^{#pm}","e^{#pm}"} ;
1782
1783       TString pnamess[] = { "Photon","Hadron","Pi0","Eta","Conversion","Electron"} ;
1784
1785       for(Int_t i = 0; i < 6; i++){
1786
1787         fhMCELambda0[i]  = new TH2F(Form("hELambda0_MC%s",pnamess[i].Data()),
1788                                     Form("cluster from %s : E vs #lambda_{0}^{2}",ptypess[i].Data()),
1789                                     nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1790         fhMCELambda0[i]->SetYTitle("#lambda_{0}^{2}");
1791         fhMCELambda0[i]->SetXTitle("E (GeV)");
1792         outputContainer->Add(fhMCELambda0[i]) ;
1793
1794         fhMCELambda1[i]  = new TH2F(Form("hELambda1_MC%s",pnamess[i].Data()),
1795                                     Form("cluster from %s : E vs #lambda_{1}^{2}",ptypess[i].Data()),
1796                                     nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1797         fhMCELambda1[i]->SetYTitle("#lambda_{1}^{2}");
1798         fhMCELambda1[i]->SetXTitle("E (GeV)");
1799         outputContainer->Add(fhMCELambda1[i]) ;
1800
1801         fhMCEDispersion[i]  = new TH2F(Form("hEDispersion_MC%s",pnamess[i].Data()),
1802                                        Form("cluster from %s : E vs dispersion^{2}",ptypess[i].Data()),
1803                                        nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1804         fhMCEDispersion[i]->SetYTitle("D^{2}");
1805         fhMCEDispersion[i]->SetXTitle("E (GeV)");
1806         outputContainer->Add(fhMCEDispersion[i]) ;
1807
1808         fhMCNCellsE[i]  = new TH2F (Form("hNCellsE_MC%s",pnamess[i].Data()),
1809                                Form("# of cells in cluster from %s vs E of clusters",ptypess[i].Data()),
1810                                     nptbins,ptmin,ptmax, nbins,nmin,nmax);
1811         fhMCNCellsE[i]->SetXTitle("E (GeV)");
1812         fhMCNCellsE[i]->SetYTitle("# of cells in cluster");
1813         outputContainer->Add(fhMCNCellsE[i]);
1814
1815         fhMCMaxCellDiffClusterE[i]  = new TH2F (Form("hMaxCellDiffClusterE_MC%s",pnamess[i].Data()),
1816                                              Form("energy vs difference of cluster energy from %s - max cell energy / cluster energy, good clusters",ptypess[i].Data()),
1817                                            nptbins,ptmin,ptmax, 500,0,1.);
1818         fhMCMaxCellDiffClusterE[i]->SetXTitle("E_{cluster} (GeV) ");
1819         fhMCMaxCellDiffClusterE[i]->SetYTitle("(E_{cluster} - E_{cell max})/ E_{cluster}");
1820         outputContainer->Add(fhMCMaxCellDiffClusterE[i]);
1821
1822         fhMCLambda0vsClusterMaxCellDiffE0[i]  = new TH2F(Form("hLambda0vsClusterMaxCellDiffE0_MC%s",pnamess[i].Data()),
1823                                     Form("cluster from %s : #lambda^{2}_{0} vs fraction of energy carried by max cell, E < 2 GeV",ptypess[i].Data()),
1824                                     ssbins,ssmin,ssmax,500,0,1.);
1825         fhMCLambda0vsClusterMaxCellDiffE0[i]->SetXTitle("#lambda_{0}^{2}");
1826         fhMCLambda0vsClusterMaxCellDiffE0[i]->SetYTitle("(E_{cluster} - E_{cell max})/ E_{cluster}");
1827         outputContainer->Add(fhMCLambda0vsClusterMaxCellDiffE0[i]) ;
1828
1829         fhMCLambda0vsClusterMaxCellDiffE2[i]  = new TH2F(Form("hLambda0vsClusterMaxCellDiffE2_MC%s",pnamess[i].Data()),
1830                                                          Form("cluster from %s : #lambda^{2}_{0} vs fraction of energy carried by max cell, 2< E < 6 GeV",ptypess[i].Data()),
1831                                                          ssbins,ssmin,ssmax,500,0,1.);
1832         fhMCLambda0vsClusterMaxCellDiffE2[i]->SetXTitle("#lambda_{0}^{2}");
1833         fhMCLambda0vsClusterMaxCellDiffE2[i]->SetYTitle("(E_{cluster} - E_{cell max})/ E_{cluster}");
1834         outputContainer->Add(fhMCLambda0vsClusterMaxCellDiffE2[i]) ;
1835
1836         fhMCLambda0vsClusterMaxCellDiffE6[i]  = new TH2F(Form("hLambda0vsClusterMaxCellDiffE6_MC%s",pnamess[i].Data()),
1837                                                          Form("cluster from %s : #lambda^{2}_{0} vs fraction of energy carried by max cell, E > 6 GeV",ptypess[i].Data()),
1838                                                          ssbins,ssmin,ssmax,500,0,1.);
1839         fhMCLambda0vsClusterMaxCellDiffE6[i]->SetXTitle("#lambda_{0}^{2}");
1840         fhMCLambda0vsClusterMaxCellDiffE6[i]->SetYTitle("(E_{cluster} - E_{cell max})/ E_{cluster}");
1841         outputContainer->Add(fhMCLambda0vsClusterMaxCellDiffE6[i]) ;
1842
1843         fhMCNCellsvsClusterMaxCellDiffE0[i]  = new TH2F(Form("hNCellsvsClusterMaxCellDiffE0_MC%s",pnamess[i].Data()),
1844                                                          Form("cluster from %s : N cells in cluster vs fraction of energy carried by max cell, E < 2 GeV",ptypess[i].Data()),
1845                                                          nbins/5,nmin,nmax/5,500,0,1.);
1846         fhMCNCellsvsClusterMaxCellDiffE0[i]->SetXTitle("N cells in cluster");
1847         fhMCNCellsvsClusterMaxCellDiffE0[i]->SetYTitle("(E_{cluster} - E_{cell max})/ E_{cluster}");
1848         outputContainer->Add(fhMCNCellsvsClusterMaxCellDiffE0[i]) ;
1849
1850         fhMCNCellsvsClusterMaxCellDiffE2[i]  = new TH2F(Form("hNCellsvsClusterMaxCellDiffE2_MC%s",pnamess[i].Data()),
1851                                                          Form("cluster from %s : N cells in cluster vs fraction of energy carried by max cell, 2< E < 6 GeV",ptypess[i].Data()),
1852                                                          nbins/5,nmin,nmax/5,500,0,1.);
1853         fhMCNCellsvsClusterMaxCellDiffE2[i]->SetXTitle("N cells in cluster");
1854         fhMCNCellsvsClusterMaxCellDiffE2[i]->SetYTitle("(E_{cluster} - E_{cell max})/ E_{cluster}");
1855         outputContainer->Add(fhMCNCellsvsClusterMaxCellDiffE2[i]) ;
1856
1857         fhMCNCellsvsClusterMaxCellDiffE6[i]  = new TH2F(Form("hNCellsvsClusterMaxCellDiffE6_MC%s",pnamess[i].Data()),
1858                                                          Form("cluster from %s : N cells in cluster vs fraction of energy carried by max cell, E > 6 GeV",ptypess[i].Data()),
1859                                                          nbins/5,nmin,nmax/5,500,0,1.);
1860         fhMCNCellsvsClusterMaxCellDiffE6[i]->SetXTitle("N cells in cluster");
1861         fhMCNCellsvsClusterMaxCellDiffE6[i]->SetYTitle("E (GeV)");
1862         outputContainer->Add(fhMCNCellsvsClusterMaxCellDiffE6[i]) ;
1863
1864        }// loop
1865
1866       if(!GetReader()->IsEmbeddedClusterSelectionOn())
1867       {
1868         fhMCPhotonELambda0NoOverlap  = new TH2F("hELambda0_MCPhoton_NoOverlap",
1869                                                 "cluster from Photon : E vs #lambda_{0}^{2}",
1870                                                 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1871         fhMCPhotonELambda0NoOverlap->SetYTitle("#lambda_{0}^{2}");
1872         fhMCPhotonELambda0NoOverlap->SetXTitle("E (GeV)");
1873         outputContainer->Add(fhMCPhotonELambda0NoOverlap) ;
1874
1875         fhMCPhotonELambda0TwoOverlap  = new TH2F("hELambda0_MCPhoton_TwoOverlap",
1876                                                  "cluster from Photon : E vs #lambda_{0}^{2}",
1877                                                  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1878         fhMCPhotonELambda0TwoOverlap->SetYTitle("#lambda_{0}^{2}");
1879         fhMCPhotonELambda0TwoOverlap->SetXTitle("E (GeV)");
1880         outputContainer->Add(fhMCPhotonELambda0TwoOverlap) ;
1881
1882         fhMCPhotonELambda0NOverlap  = new TH2F("hELambda0_MCPhoton_NOverlap",
1883                                                "cluster from Photon : E vs #lambda_{0}^{2}",
1884                                                nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1885         fhMCPhotonELambda0NOverlap->SetYTitle("#lambda_{0}^{2}");
1886         fhMCPhotonELambda0NOverlap->SetXTitle("E (GeV)");
1887         outputContainer->Add(fhMCPhotonELambda0NOverlap) ;
1888
1889       } //No embedding
1890
1891       //Fill histograms to check shape of embedded clusters
1892       if(GetReader()->IsEmbeddedClusterSelectionOn())
1893       {
1894
1895         fhEmbeddedSignalFractionEnergy  = new TH2F("hEmbeddedSignal_FractionEnergy",
1896                                                     "Energy Fraction of embedded signal versus cluster energy",
1897                                                     nptbins,ptmin,ptmax,100,0.,1.);
1898         fhEmbeddedSignalFractionEnergy->SetYTitle("Fraction");
1899         fhEmbeddedSignalFractionEnergy->SetXTitle("E (GeV)");
1900         outputContainer->Add(fhEmbeddedSignalFractionEnergy) ;
1901
1902         fhEmbedPhotonELambda0FullSignal  = new TH2F("hELambda0_EmbedPhoton_FullSignal",
1903                                                 "cluster from Photon embedded with more than 90% energy in cluster : E vs #lambda_{0}^{2}",
1904                                                 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1905         fhEmbedPhotonELambda0FullSignal->SetYTitle("#lambda_{0}^{2}");
1906         fhEmbedPhotonELambda0FullSignal->SetXTitle("E (GeV)");
1907         outputContainer->Add(fhEmbedPhotonELambda0FullSignal) ;
1908
1909         fhEmbedPhotonELambda0MostlySignal  = new TH2F("hELambda0_EmbedPhoton_MostlySignal",
1910                                           "cluster from Photon embedded with 50% to 90% energy in cluster : E vs #lambda_{0}^{2}",
1911                                           nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1912         fhEmbedPhotonELambda0MostlySignal->SetYTitle("#lambda_{0}^{2}");
1913         fhEmbedPhotonELambda0MostlySignal->SetXTitle("E (GeV)");
1914         outputContainer->Add(fhEmbedPhotonELambda0MostlySignal) ;
1915
1916         fhEmbedPhotonELambda0MostlyBkg  = new TH2F("hELambda0_EmbedPhoton_MostlyBkg",
1917                                           "cluster from Photon embedded with 10% to 50% energy in cluster : E vs #lambda_{0}^{2}",
1918                                           nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1919         fhEmbedPhotonELambda0MostlyBkg->SetYTitle("#lambda_{0}^{2}");
1920         fhEmbedPhotonELambda0MostlyBkg->SetXTitle("E (GeV)");
1921         outputContainer->Add(fhEmbedPhotonELambda0MostlyBkg) ;
1922
1923         fhEmbedPhotonELambda0FullBkg  = new TH2F("hELambda0_EmbedPhoton_FullBkg",
1924                                                    "cluster from Photonm embedded with 0% to 10% energy in cluster : E vs #lambda_{0}^{2}",
1925                                                    nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1926         fhEmbedPhotonELambda0FullBkg->SetYTitle("#lambda_{0}^{2}");
1927         fhEmbedPhotonELambda0FullBkg->SetXTitle("E (GeV)");
1928         outputContainer->Add(fhEmbedPhotonELambda0FullBkg) ;
1929
1930         fhEmbedPi0ELambda0FullSignal  = new TH2F("hELambda0_EmbedPi0_FullSignal",
1931                                                     "cluster from Pi0 embedded with more than 90% energy in cluster : E vs #lambda_{0}^{2}",
1932                                                     nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1933         fhEmbedPi0ELambda0FullSignal->SetYTitle("#lambda_{0}^{2}");
1934         fhEmbedPi0ELambda0FullSignal->SetXTitle("E (GeV)");
1935         outputContainer->Add(fhEmbedPi0ELambda0FullSignal) ;
1936
1937         fhEmbedPi0ELambda0MostlySignal  = new TH2F("hELambda0_EmbedPi0_MostlySignal",
1938                                                       "cluster from Pi0 embedded with 50% to 90% energy in cluster : E vs #lambda_{0}^{2}",
1939                                                       nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1940         fhEmbedPi0ELambda0MostlySignal->SetYTitle("#lambda_{0}^{2}");
1941         fhEmbedPi0ELambda0MostlySignal->SetXTitle("E (GeV)");
1942         outputContainer->Add(fhEmbedPi0ELambda0MostlySignal) ;
1943
1944         fhEmbedPi0ELambda0MostlyBkg  = new TH2F("hELambda0_EmbedPi0_MostlyBkg",
1945                                                    "cluster from Pi0 embedded with 10% to 50% energy in cluster : E vs #lambda_{0}^{2}",
1946                                                    nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1947         fhEmbedPi0ELambda0MostlyBkg->SetYTitle("#lambda_{0}^{2}");
1948         fhEmbedPi0ELambda0MostlyBkg->SetXTitle("E (GeV)");
1949         outputContainer->Add(fhEmbedPi0ELambda0MostlyBkg) ;
1950
1951         fhEmbedPi0ELambda0FullBkg  = new TH2F("hELambda0_EmbedPi0_FullBkg",
1952                                                  "cluster from Pi0 embedded with 0% to 10% energy in cluster : E vs #lambda_{0}^{2}",
1953                                                  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1954         fhEmbedPi0ELambda0FullBkg->SetYTitle("#lambda_{0}^{2}");
1955         fhEmbedPi0ELambda0FullBkg->SetXTitle("E (GeV)");
1956         outputContainer->Add(fhEmbedPi0ELambda0FullBkg) ;
1957
1958       }// embedded histograms
1959
1960
1961     }// Fill SS MC histograms
1962
1963   }//Histos with MC
1964
1965   return outputContainer ;
1966
1967 }
1968
1969 //____________________________________________________________________________
1970 void AliAnaPhoton::Init()
1971 {
1972
1973   //Init
1974   //Do some checks
1975   if(fCalorimeter == "PHOS" && !GetReader()->IsPHOSSwitchedOn() && NewOutputAOD()){
1976     printf("AliAnaPhoton::Init() - !!STOP: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!\n");
1977     abort();
1978   }
1979   else  if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn() && NewOutputAOD()){
1980     printf("AliAnaPhoton::Init() - !!STOP: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!\n");
1981     abort();
1982   }
1983
1984   if(GetReader()->GetDataType() == AliCaloTrackReader::kMC) GetCaloPID()->SwitchOnBayesian();
1985
1986 }
1987
1988 //____________________________________________________________________________
1989 void AliAnaPhoton::InitParameters()
1990 {
1991
1992   //Initialize the parameters of the analysis.
1993   AddToHistogramsName("AnaPhoton_");
1994
1995   fCalorimeter = "EMCAL" ;
1996   fMinDist     = 2.;
1997   fMinDist2    = 4.;
1998   fMinDist3    = 5.;
1999
2000   fTimeCutMin  =-1000000;
2001   fTimeCutMax  = 1000000;
2002   fNCellsCut   = 0;
2003
2004   fRejectTrackMatch       = kTRUE ;
2005
2006 }
2007
2008 //__________________________________________________________________
2009 void  AliAnaPhoton::MakeAnalysisFillAOD()
2010 {
2011   //Do photon analysis and fill aods
2012
2013   //Get the vertex
2014   Double_t v[3] = {0,0,0}; //vertex ;
2015   GetReader()->GetVertex(v);
2016
2017   //Select the Calorimeter of the photon
2018   TObjArray * pl = 0x0;
2019   if(fCalorimeter == "PHOS")
2020     pl = GetPHOSClusters();
2021   else if (fCalorimeter == "EMCAL")
2022     pl = GetEMCALClusters();
2023
2024   if(!pl) {
2025     Info("MakeAnalysisFillAOD","TObjArray with %s clusters is NULL!\n",fCalorimeter.Data());
2026     return;
2027   }
2028
2029   // Loop on raw clusters before filtering in the reader and fill control histogram
2030   if((GetReader()->GetEMCALClusterListName()=="" && fCalorimeter=="EMCAL") || fCalorimeter=="PHOS"){
2031     for(Int_t iclus = 0; iclus < GetReader()->GetInputEvent()->GetNumberOfCaloClusters(); iclus++ ){
2032       AliVCluster * clus = GetReader()->GetInputEvent()->GetCaloCluster(iclus);
2033       if     (fCalorimeter == "PHOS"  && clus->IsPHOS()  && clus->E() > GetReader()->GetPHOSPtMin() ) fhClusterCuts[0]->Fill(clus->E());
2034       else if(fCalorimeter == "EMCAL" && clus->IsEMCAL() && clus->E() > GetReader()->GetEMCALPtMin()) fhClusterCuts[0]->Fill(clus->E());
2035     }
2036   }
2037   else { // reclusterized
2038     TClonesArray * clusterList = 0;
2039     if(GetReader()->GetOutputEvent())
2040       clusterList = dynamic_cast<TClonesArray*> (GetReader()->GetOutputEvent()->FindListObject(GetReader()->GetEMCALClusterListName()));
2041     if(clusterList){
2042       Int_t nclusters = clusterList->GetEntriesFast();
2043       for (Int_t iclus =  0; iclus <  nclusters; iclus++) {
2044         AliVCluster * clus = dynamic_cast<AliVCluster*> (clusterList->At(iclus));
2045         if(clus)fhClusterCuts[0]->Fill(clus->E());
2046       }
2047     }
2048   }
2049
2050   //Init arrays, variables, get number of clusters
2051   TLorentzVector mom, mom2 ;
2052   Int_t nCaloClusters = pl->GetEntriesFast();
2053
2054   if(GetDebug() > 0) printf("AliAnaPhoton::MakeAnalysisFillAOD() - input %s cluster entries %d\n", fCalorimeter.Data(), nCaloClusters);
2055
2056   //----------------------------------------------------
2057   // Fill AOD with PHOS/EMCAL AliAODPWG4Particle objects
2058   //----------------------------------------------------
2059   // Loop on clusters
2060   for(Int_t icalo = 0; icalo < nCaloClusters; icalo++){
2061
2062           AliVCluster * calo =  (AliVCluster*) (pl->At(icalo));
2063     //printf("calo %d, %f\n",icalo,calo->E());
2064
2065     //Get the index where the cluster comes, to retrieve the corresponding vertex
2066     Int_t evtIndex = 0 ;
2067     if (GetMixedEvent()) {
2068       evtIndex=GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ;
2069       //Get the vertex and check it is not too large in z
2070       if(TMath::Abs(GetVertex(evtIndex)[2])> GetZvertexCut()) continue;
2071     }
2072
2073     //Cluster selection, not charged, with photon id and in fiducial cut
2074     if(GetReader()->GetDataType() != AliCaloTrackReader::kMC){
2075       calo->GetMomentum(mom,GetVertex(evtIndex)) ;}//Assume that come from vertex in straight line
2076     else{
2077       Double_t vertex[]={0,0,0};
2078       calo->GetMomentum(mom,vertex) ;
2079     }
2080
2081     //--------------------------------------
2082     // Cluster selection
2083     //--------------------------------------
2084     if(!ClusterSelected(calo,mom)) continue;
2085
2086     //----------------------------
2087     //Create AOD for analysis
2088     //----------------------------
2089     AliAODPWG4Particle aodph = AliAODPWG4Particle(mom);
2090
2091     //...............................................
2092     //Set the indeces of the original caloclusters (MC, ID), and calorimeter
2093     Int_t label = calo->GetLabel();
2094     aodph.SetLabel(label);
2095     aodph.SetCaloLabel(calo->GetID(),-1);
2096     aodph.SetDetector(fCalorimeter);
2097     //printf("Index %d, Id %d, iaod %d\n",icalo, calo->GetID(),GetOutputAODBranch()->GetEntriesFast());
2098
2099     //...............................................
2100     //Set bad channel distance bit
2101     Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
2102     if     (distBad > fMinDist3) aodph.SetDistToBad(2) ;
2103     else if(distBad > fMinDist2) aodph.SetDistToBad(1) ;
2104     else                         aodph.SetDistToBad(0) ;
2105     //printf("DistBad %f Bit %d\n",distBad, aodph.DistToBad());
2106
2107     //--------------------------------------------------------------------------------------
2108     // Play with the MC stack if available
2109     //--------------------------------------------------------------------------------------
2110
2111     //Check origin of the candidates
2112     Int_t tag = -1;
2113
2114     if(IsDataMC()){
2115
2116       tag = GetMCAnalysisUtils()->CheckOrigin(calo->GetLabels(),calo->GetNLabels(),GetReader(), aodph.GetInputFileIndex());
2117       aodph.SetTag(tag);
2118
2119       if(GetDebug() > 0)
2120         printf("AliAnaPhoton::MakeAnalysisFillAOD() - Origin of candidate, bit map %d\n",aodph.GetTag());
2121     }//Work with stack also
2122
2123
2124     //--------------------------------------------------------------------------------------
2125     //Fill some shower shape histograms before PID is applied
2126     //--------------------------------------------------------------------------------------
2127
2128     FillShowerShapeHistograms(calo,tag);
2129
2130     //-------------------------------------
2131     //PID selection or bit setting
2132     //-------------------------------------
2133
2134     //...............................................
2135     // Data, PID check on
2136     if(IsCaloPIDOn()){
2137       // Get most probable PID, 2 options check bayesian PID weights or redo PID
2138       // By default, redo PID
2139
2140       aodph.SetIdentifiedParticleType(GetCaloPID()->GetIdentifiedParticleType(fCalorimeter,mom,calo));
2141
2142       if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PDG of identified particle %d\n",aodph.GetIdentifiedParticleType());
2143
2144       //If cluster does not pass pid, not photon, skip it.
2145       if(aodph.GetIdentifiedParticleType() != AliCaloPID::kPhoton) continue ;
2146
2147     }
2148
2149     //...............................................
2150     // Data, PID check off
2151     else{
2152       //Set PID bits for later selection (AliAnaPi0 for example)
2153       //GetIdentifiedParticleType already called in SetPIDBits.
2154
2155       GetCaloPID()->SetPIDBits(fCalorimeter,calo,&aodph, GetCaloUtils(),GetReader()->GetInputEvent());
2156
2157       if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PID Bits set \n");
2158     }
2159
2160     if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - Photon selection cuts passed: pT %3.2f, pdg %d\n",aodph.Pt(), aodph.GetIdentifiedParticleType());
2161
2162     fhClusterCuts[8]->Fill(calo->E());
2163
2164     // Matching after cuts
2165     if(fFillTMHisto) FillTrackMatchingResidualHistograms(calo,1);
2166
2167     // Add number of local maxima to AOD, method name in AOD to be FIXED
2168     AliVCaloCells* cells    = 0;
2169     if(fCalorimeter == "EMCAL") cells    = GetEMCALCells();
2170     else                        cells    = GetPHOSCells();
2171     const Int_t   nc = calo->GetNCells();
2172     Int_t   *absIdList = new Int_t  [nc];
2173     Float_t *maxEList  = new Float_t[nc];
2174
2175     aodph.SetFiducialArea(GetCaloUtils()->GetNumberOfLocalMaxima(calo, cells, absIdList, maxEList));
2176
2177     delete [] absIdList;
2178     delete [] maxEList;
2179
2180     //Add AOD with photon object to aod branch
2181     AddAODParticle(aodph);
2182
2183   }//loop
2184
2185   if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD()  End fill AODs, with %d entries \n",GetOutputAODBranch()->GetEntriesFast());
2186
2187 }
2188
2189 //__________________________________________________________________
2190 void  AliAnaPhoton::MakeAnalysisFillHistograms()
2191 {
2192   //Fill histograms
2193
2194   //-------------------------------------------------------------------
2195   // Access MC information in stack if requested, check that it exists.
2196   AliStack         * stack       = 0x0;
2197   TParticle        * primary     = 0x0;
2198   TClonesArray     * mcparticles = 0x0;
2199   AliAODMCParticle * aodprimary  = 0x0;
2200
2201   if(IsDataMC()){
2202
2203     if(GetReader()->ReadStack()){
2204       stack =  GetMCStack() ;
2205       if(!stack) {
2206         printf("AliAnaPhoton::MakeAnalysisFillHistograms() - Stack not available, is the MC handler called? STOP\n");
2207         abort();
2208       }
2209
2210     }
2211     else if(GetReader()->ReadAODMCParticles()){
2212
2213       //Get the list of MC particles
2214       mcparticles = GetReader()->GetAODMCParticles(0);
2215       if(!mcparticles && GetDebug() > 0)        {
2216         printf("AliAnaPhoton::MakeAnalysisFillHistograms() -  Standard MCParticles not available!\n");
2217       }
2218     }
2219   }// is data and MC
2220
2221
2222   // Get vertex
2223   Double_t v[3] = {0,0,0}; //vertex ;
2224   GetReader()->GetVertex(v);
2225   //fhVertex->Fill(v[0],v[1],v[2]);
2226   if(TMath::Abs(v[2]) > GetZvertexCut()) return ; // done elsewhere for Single Event analysis, but there for mixed event
2227
2228   //----------------------------------
2229   //Loop on stored AOD photons
2230   Int_t naod = GetOutputAODBranch()->GetEntriesFast();
2231   if(GetDebug() > 0) printf("AliAnaPhoton::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod);
2232
2233   for(Int_t iaod = 0; iaod < naod ; iaod++){
2234     AliAODPWG4Particle* ph =  (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
2235     Int_t pdg = ph->GetIdentifiedParticleType();
2236
2237     if(GetDebug() > 3)
2238       printf("AliAnaPhoton::MakeAnalysisFillHistograms() - PDG %d, MC TAG %d, Calorimeter %s\n", ph->GetIdentifiedParticleType(),ph->GetTag(), (ph->GetDetector()).Data()) ;
2239
2240     //If PID used, fill histos with photons in Calorimeter fCalorimeter
2241     if(IsCaloPIDOn() && pdg != AliCaloPID::kPhoton) continue;
2242     if(ph->GetDetector() != fCalorimeter) continue;
2243
2244     if(GetDebug() > 2)
2245       printf("AliAnaPhoton::MakeAnalysisFillHistograms() - ID Photon: pt %f, phi %f, eta %f\n", ph->Pt(),ph->Phi(),ph->Eta()) ;
2246
2247     //................................
2248     //Fill photon histograms
2249     Float_t ptcluster  = ph->Pt();
2250     Float_t phicluster = ph->Phi();
2251     Float_t etacluster = ph->Eta();
2252     Float_t ecluster   = ph->E();
2253
2254     fhEPhoton   ->Fill(ecluster);
2255     fhPtPhoton  ->Fill(ptcluster);
2256     fhPhiPhoton ->Fill(ptcluster,phicluster);
2257     fhEtaPhoton ->Fill(ptcluster,etacluster);
2258     if     (ecluster > 0.5)   fhEtaPhiPhoton  ->Fill(etacluster, phicluster);
2259     else if(GetMinPt() < 0.5) fhEtaPhi05Photon->Fill(etacluster, phicluster);
2260
2261
2262     //Get original cluster, to recover some information
2263     Int_t absID             = 0;
2264     Float_t maxCellFraction = 0;
2265     AliVCaloCells* cells    = 0;
2266     TObjArray * clusters    = 0;
2267     if(fCalorimeter == "EMCAL"){
2268       cells    = GetEMCALCells();
2269       clusters = GetEMCALClusters();
2270     }
2271     else{
2272       cells    = GetPHOSCells();
2273       clusters = GetPHOSClusters();
2274     }
2275
2276     Int_t iclus = -1;
2277     AliVCluster *cluster = FindCluster(clusters,ph->GetCaloLabel(0),iclus);
2278     if(cluster)
2279     {
2280       absID = GetCaloUtils()->GetMaxEnergyCell(cells, cluster,maxCellFraction);
2281
2282       // Control histograms
2283       fhMaxCellDiffClusterE->Fill(ph->E(),maxCellFraction);
2284       fhNCellsE            ->Fill(ph->E(),cluster->GetNCells());
2285       fhTimeE              ->Fill(ph->E(),cluster->GetTOF()*1.e9);
2286       if(cells)
2287       {
2288       for(Int_t icell = 0; icell <  cluster->GetNCells(); icell++)
2289         fhCellsE->Fill(ph->E(),cells->GetCellAmplitude(cluster->GetCellsAbsId()[icell]));
2290       }
2291     }
2292
2293     //.......................................
2294     //Play with the MC data if available
2295     if(IsDataMC()){
2296
2297       FillAcceptanceHistograms();
2298
2299       //....................................................................
2300       // Access MC information in stack if requested, check that it exists.
2301       Int_t label =ph->GetLabel();
2302       if(label < 0) {
2303         if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** bad label ***:  label %d \n", label);
2304         continue;
2305       }
2306
2307       Float_t eprim   = 0;
2308       Float_t ptprim  = 0;
2309       if(GetReader()->ReadStack()){
2310
2311         if(label >=  stack->GetNtrack()) {
2312           if(GetDebug() > 2)  printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** large label ***:  label %d, n tracks %d \n", label, stack->GetNtrack());
2313           continue ;
2314         }
2315
2316         primary = stack->Particle(label);
2317         if(!primary){
2318           printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** no primary ***:  label %d \n", label);
2319           continue;
2320         }
2321         eprim   = primary->Energy();
2322         ptprim  = primary->Pt();
2323
2324       }
2325       else if(GetReader()->ReadAODMCParticles()){
2326         //Check which is the input
2327         if(ph->GetInputFileIndex() == 0){
2328           if(!mcparticles) continue;
2329           if(label >=  mcparticles->GetEntriesFast()) {
2330             if(GetDebug() > 2)  printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** large label ***:  label %d, n tracks %d \n",
2331                                        label, mcparticles->GetEntriesFast());
2332             continue ;
2333           }
2334           //Get the particle
2335           aodprimary = (AliAODMCParticle*) mcparticles->At(label);
2336
2337         }
2338
2339         if(!aodprimary){
2340           printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** no primary ***:  label %d \n", label);
2341           continue;
2342         }
2343
2344         eprim   = aodprimary->E();
2345         ptprim  = aodprimary->Pt();
2346
2347       }
2348
2349       Int_t tag =ph->GetTag();
2350
2351       if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) && fhMCE[kmcPhoton])
2352       {
2353         fhMCE  [kmcPhoton] ->Fill(ecluster);
2354         fhMCPt [kmcPhoton] ->Fill(ptcluster);
2355         fhMCPhi[kmcPhoton] ->Fill(ecluster,phicluster);
2356         fhMCEta[kmcPhoton] ->Fill(ecluster,etacluster);
2357
2358         fhMC2E[kmcPhoton]     ->Fill(ecluster, eprim);
2359         fhMC2Pt[kmcPhoton]    ->Fill(ptcluster, ptprim);
2360         fhMCDeltaE[kmcPhoton] ->Fill(ecluster,eprim-ecluster);
2361         fhMCDeltaPt[kmcPhoton]->Fill(ptcluster,ptprim-ptcluster);
2362
2363         if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion) && fhMCE[kmcConversion])
2364         {
2365           fhMCE  [kmcConversion] ->Fill(ecluster);
2366           fhMCPt [kmcConversion] ->Fill(ptcluster);
2367           fhMCPhi[kmcConversion] ->Fill(ecluster,phicluster);
2368           fhMCEta[kmcConversion] ->Fill(ecluster,etacluster);
2369
2370           fhMC2E[kmcConversion]     ->Fill(ecluster, eprim);
2371           fhMC2Pt[kmcConversion]    ->Fill(ptcluster, ptprim);
2372           fhMCDeltaE[kmcConversion] ->Fill(ecluster,eprim-ecluster);
2373           fhMCDeltaPt[kmcConversion]->Fill(ptcluster,ptprim-ptcluster);
2374         }
2375
2376         if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt) && fhMCE[kmcPrompt]){
2377           fhMCE  [kmcPrompt] ->Fill(ecluster);
2378           fhMCPt [kmcPrompt] ->Fill(ptcluster);
2379           fhMCPhi[kmcPrompt] ->Fill(ecluster,phicluster);
2380           fhMCEta[kmcPrompt] ->Fill(ecluster,etacluster);
2381
2382           fhMC2E[kmcPrompt]     ->Fill(ecluster, eprim);
2383           fhMC2Pt[kmcPrompt]    ->Fill(ptcluster, ptprim);
2384           fhMCDeltaE[kmcPrompt] ->Fill(ecluster,eprim-ecluster);
2385           fhMCDeltaPt[kmcPrompt]->Fill(ptcluster,ptprim-ptcluster);
2386
2387         }
2388         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation)&& fhMCE[kmcFragmentation])
2389         {
2390           fhMCE  [kmcFragmentation] ->Fill(ecluster);
2391           fhMCPt [kmcFragmentation] ->Fill(ptcluster);
2392           fhMCPhi[kmcFragmentation] ->Fill(ecluster,phicluster);
2393           fhMCEta[kmcFragmentation] ->Fill(ecluster,etacluster);
2394
2395           fhMC2E[kmcFragmentation]     ->Fill(ecluster, eprim);
2396           fhMC2Pt[kmcFragmentation]    ->Fill(ptcluster, ptprim);
2397           fhMCDeltaE[kmcFragmentation] ->Fill(ecluster,eprim-ecluster);
2398           fhMCDeltaPt[kmcFragmentation]->Fill(ptcluster,ptprim-ptcluster);
2399
2400         }
2401         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCISR)&& fhMCE[kmcISR])
2402         {
2403           fhMCE  [kmcISR] ->Fill(ecluster);
2404           fhMCPt [kmcISR] ->Fill(ptcluster);
2405           fhMCPhi[kmcISR] ->Fill(ecluster,phicluster);
2406           fhMCEta[kmcISR] ->Fill(ecluster,etacluster);
2407
2408           fhMC2E[kmcISR]     ->Fill(ecluster, eprim);
2409           fhMC2Pt[kmcISR]    ->Fill(ptcluster, ptprim);
2410           fhMCDeltaE[kmcISR] ->Fill(ecluster,eprim-ecluster);
2411           fhMCDeltaPt[kmcISR]->Fill(ptcluster,ptprim-ptcluster);
2412
2413         }
2414         else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay) &&
2415                 !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) && fhMCE[kmcPi0Decay])
2416         {
2417           fhMCE  [kmcPi0Decay] ->Fill(ecluster);
2418           fhMCPt [kmcPi0Decay] ->Fill(ptcluster);
2419           fhMCPhi[kmcPi0Decay] ->Fill(ecluster,phicluster);
2420           fhMCEta[kmcPi0Decay] ->Fill(ecluster,etacluster);
2421
2422           fhMC2E[kmcPi0Decay]     ->Fill(ecluster, eprim);
2423           fhMC2Pt[kmcPi0Decay]    ->Fill(ptcluster, ptprim);
2424           fhMCDeltaE[kmcPi0Decay] ->Fill(ecluster,eprim-ecluster);
2425           fhMCDeltaPt[kmcPi0Decay]->Fill(ptcluster,ptprim-ptcluster);
2426
2427         }
2428         else if( (GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay) ||
2429                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay) ) && fhMCE[kmcOtherDecay])
2430         {
2431           fhMCE  [kmcOtherDecay] ->Fill(ecluster);
2432           fhMCPt [kmcOtherDecay] ->Fill(ptcluster);
2433           fhMCPhi[kmcOtherDecay] ->Fill(ecluster,phicluster);
2434           fhMCEta[kmcOtherDecay] ->Fill(ecluster,etacluster);
2435
2436           fhMC2E[kmcOtherDecay]     ->Fill(ecluster, eprim);
2437           fhMC2Pt[kmcOtherDecay]    ->Fill(ptcluster, ptprim);
2438           fhMCDeltaE[kmcOtherDecay] ->Fill(ecluster,eprim-ecluster);
2439           fhMCDeltaPt[kmcOtherDecay]->Fill(ptcluster,ptprim-ptcluster);
2440
2441         }
2442         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) && fhMCE  [kmcPi0])
2443         {
2444           fhMCE  [kmcPi0] ->Fill(ecluster);
2445           fhMCPt [kmcPi0] ->Fill(ptcluster);
2446           fhMCPhi[kmcPi0] ->Fill(ecluster,phicluster);
2447           fhMCEta[kmcPi0] ->Fill(ecluster,etacluster);
2448
2449           fhMC2E[kmcPi0]     ->Fill(ecluster, eprim);
2450           fhMC2Pt[kmcPi0]    ->Fill(ptcluster, ptprim);
2451           fhMCDeltaE[kmcPi0] ->Fill(ecluster,eprim-ecluster);
2452           fhMCDeltaPt[kmcPi0]->Fill(ptcluster,ptprim-ptcluster);
2453
2454         }
2455         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta) && fhMCE[kmcEta])
2456         {
2457           fhMCE  [kmcEta] ->Fill(ecluster);
2458           fhMCPt [kmcEta] ->Fill(ptcluster);
2459           fhMCPhi[kmcEta] ->Fill(ecluster,phicluster);
2460           fhMCEta[kmcEta] ->Fill(ecluster,etacluster);
2461
2462           fhMC2E[kmcEta]     ->Fill(ecluster, eprim);
2463           fhMC2Pt[kmcEta]    ->Fill(ptcluster, ptprim);
2464           fhMCDeltaE[kmcEta] ->Fill(ecluster,eprim-ecluster);
2465           fhMCDeltaPt[kmcEta]->Fill(ptcluster,ptprim-ptcluster);
2466
2467         }
2468       }
2469       else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiNeutron) && fhMCE[kmcAntiNeutron])
2470       {
2471         fhMCE  [kmcAntiNeutron] ->Fill(ecluster);
2472         fhMCPt [kmcAntiNeutron] ->Fill(ptcluster);
2473         fhMCPhi[kmcAntiNeutron] ->Fill(ecluster,phicluster);
2474         fhMCEta[kmcAntiNeutron] ->Fill(ecluster,etacluster);
2475
2476         fhMC2E[kmcAntiNeutron]     ->Fill(ecluster, eprim);
2477         fhMC2Pt[kmcAntiNeutron]    ->Fill(ptcluster, ptprim);
2478         fhMCDeltaE[kmcAntiNeutron] ->Fill(ecluster,eprim-ecluster);
2479         fhMCDeltaPt[kmcAntiNeutron]->Fill(ptcluster,ptprim-ptcluster);
2480
2481       }
2482       else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiProton) && fhMCE[kmcAntiProton])
2483       {
2484         fhMCE  [kmcAntiProton] ->Fill(ecluster);
2485         fhMCPt [kmcAntiProton] ->Fill(ptcluster);
2486         fhMCPhi[kmcAntiProton] ->Fill(ecluster,phicluster);
2487         fhMCEta[kmcAntiProton] ->Fill(ecluster,etacluster);
2488
2489         fhMC2E[kmcAntiProton]     ->Fill(ecluster, eprim);
2490         fhMC2Pt[kmcAntiProton]    ->Fill(ptcluster, ptprim);
2491         fhMCDeltaE[kmcAntiProton] ->Fill(ecluster,eprim-ecluster);
2492         fhMCDeltaPt[kmcAntiProton]->Fill(ecluster,ptprim-ptcluster);
2493
2494       }
2495       else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron) && fhMCE[kmcElectron])
2496       {
2497         fhMCE  [kmcElectron] ->Fill(ecluster);
2498         fhMCPt [kmcElectron] ->Fill(ptcluster);
2499         fhMCPhi[kmcElectron] ->Fill(ecluster,phicluster);
2500         fhMCEta[kmcElectron] ->Fill(ecluster,etacluster);
2501
2502         fhMC2E[kmcElectron]     ->Fill(ecluster, eprim);
2503         fhMC2Pt[kmcElectron]    ->Fill(ptcluster, ptprim);
2504         fhMCDeltaE[kmcElectron] ->Fill(ecluster,eprim-ecluster);
2505         fhMCDeltaPt[kmcElectron]->Fill(ecluster,ptprim-ptcluster);
2506       }
2507       else if( fhMCE[kmcOther]){
2508         fhMCE  [kmcOther] ->Fill(ecluster);
2509         fhMCPt [kmcOther] ->Fill(ptcluster);
2510         fhMCPhi[kmcOther] ->Fill(ecluster,phicluster);
2511         fhMCEta[kmcOther] ->Fill(ecluster,etacluster);
2512
2513         fhMC2E[kmcOther]     ->Fill(ecluster, eprim);
2514         fhMC2Pt[kmcOther]    ->Fill(ptcluster, ptprim);
2515         fhMCDeltaE[kmcOther] ->Fill(ecluster,eprim-ecluster);
2516         fhMCDeltaPt[kmcOther]->Fill(ecluster,ptprim-ptcluster);
2517
2518         //               printf(" AliAnaPhoton::MakeAnalysisFillHistograms() - Label %d, pT %2.3f Unknown, bits set: ",
2519         //                                      ph->GetLabel(),ph->Pt());
2520         //                for(Int_t i = 0; i < 20; i++) {
2521         //                        if(GetMCAnalysisUtils()->CheckTagBit(tag,i)) printf(" %d, ",i);
2522         //                }
2523         //                printf("\n");
2524
2525       }
2526
2527     }//Histograms with MC
2528
2529   }// aod loop
2530
2531 }
2532
2533
2534 //__________________________________________________________________
2535 void AliAnaPhoton::Print(const Option_t * opt) const
2536 {
2537   //Print some relevant parameters set for the analysis
2538
2539   if(! opt)
2540     return;
2541
2542   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
2543   AliAnaCaloTrackCorrBaseClass::Print(" ");
2544
2545   printf("Calorimeter            =     %s\n", fCalorimeter.Data()) ;
2546   printf("Min Distance to Bad Channel   = %2.1f\n",fMinDist);
2547   printf("Min Distance to Bad Channel 2 = %2.1f\n",fMinDist2);
2548   printf("Min Distance to Bad Channel 3 = %2.1f\n",fMinDist3);
2549   printf("Reject clusters with a track matched = %d\n",fRejectTrackMatch);
2550   printf("Time Cut: %3.1f < TOF  < %3.1f\n", fTimeCutMin, fTimeCutMax);
2551   printf("Number of cells in cluster is        > %d \n", fNCellsCut);
2552   printf("    \n") ;
2553
2554 }