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