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