PWGGA/CaloTrackCorrelations/AliAnaPhoton.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes hereby granted      *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 //_________________________________________________________________________
  17 //
  18 // Class for the photon identification.
  19 // Clusters from calorimeters are identified as photons
  20 // and kept in the AOD. Few histograms produced.
  21 // Produces input for other analysis classes like AliAnaPi0,
  22 // AliAnaParticleHadronCorrelation ...
  23 //
  24 // -- Author: Gustavo Conesa (LNF-INFN)
  25 //////////////////////////////////////////////////////////////////////////////
  26
  27
  28 // --- ROOT system ---
  29 #include <TH2F.h>
  30 #include <TH3D.h>
  31 #include <TClonesArray.h>
  32 #include <TObjString.h>
  33 #include "TParticle.h"
  34 #include "TDatabasePDG.h"
  35
  36 // --- Analysis system ---
  37 #include "AliAnaPhoton.h"
  38 #include "AliCaloTrackReader.h"
  39 #include "AliStack.h"
  40 #include "AliCaloPID.h"
  41 #include "AliMCAnalysisUtils.h"
  42 #include "AliFiducialCut.h"
  43 #include "AliVCluster.h"
  44 #include "AliAODMCParticle.h"
  45 #include "AliMixedEvent.h"
  46 #include "AliAODEvent.h"
  47 #include "AliESDEvent.h"
  48
  49 // --- Detectors ---
  50 #include "AliPHOSGeoUtils.h"
  51 #include "AliEMCALGeometry.h"
  52
  53 ClassImp(AliAnaPhoton)
  54
  55 //____________________________
  56 AliAnaPhoton::AliAnaPhoton() :
  57 AliAnaCaloTrackCorrBaseClass(), fCalorimeter(""),
  58 fMinDist(0.),                 fMinDist2(0.),                fMinDist3(0.),
  59 fRejectTrackMatch(0),         fFillTMHisto(kFALSE),
  60 fTimeCutMin(-10000),          fTimeCutMax(10000),
  61 fNCellsCut(0),
  62 fNLMCutMin(-1),               fNLMCutMax(10),
  63 fFillSSHistograms(kFALSE),    fFillOnlySimpleSSHisto(1),
  64 fNOriginHistograms(8),        fNPrimaryHistograms(4),
  65 fFillPileUpHistograms(0),     fFillEMCALBCHistograms(0),
  66 // Histograms
  67 fhNCellsE(0),                 fhCellsE(0),   // Control histograms
  68 fhMaxCellDiffClusterE(0),     fhTimePt(0),   // Control histograms
  69 fhEtaPhi(0),                  fhEtaPhiEMCALBC0(0),
  70 fhEtaPhiEMCALBC1(0),          fhEtaPhiEMCALBCN(0),
  71 fhTimeTriggerEMCALBCCluster(0),
  72 fhTimeTriggerEMCALBCUMCluster(0),
  73 fhEtaPhiTriggerEMCALBCClusterOverTh(0),
  74 fhEtaPhiTriggerEMCALBCUMClusterOverTh(0),
  75 fhEtaPhiTriggerEMCALBCClusterBelowTh1(0),
  76 fhEtaPhiTriggerEMCALBCUMClusterBelowTh1(0),
  77 fhEtaPhiTriggerEMCALBCClusterBelowTh2(0),
  78 fhEtaPhiTriggerEMCALBCUMClusterBelowTh2(0),
  79 fhEtaPhiTriggerEMCALBCExotic(0),             fhTimeTriggerEMCALBCExotic(0),
  80 fhEtaPhiTriggerEMCALBCUMExotic(0),           fhTimeTriggerEMCALBCUMExotic(0),
  81 fhEtaPhiTriggerEMCALBCBad(0),                fhTimeTriggerEMCALBCBad(0),
  82 fhEtaPhiTriggerEMCALBCUMBad(0),              fhTimeTriggerEMCALBCUMBad(0),
  83 fhEtaPhiTriggerEMCALBCBadExotic(0),          fhTimeTriggerEMCALBCBadExotic(0),
  84 fhEtaPhiTriggerEMCALBCUMBadExotic(0),        fhTimeTriggerEMCALBCUMBadExotic(0),
  85 fhEtaPhiTriggerEMCALBCExoticCluster(0),      fhTimeTriggerEMCALBCExoticCluster(0),
  86 fhEtaPhiTriggerEMCALBCUMExoticCluster(0),    fhTimeTriggerEMCALBCUMExoticCluster(0),
  87 fhEtaPhiTriggerEMCALBCBadCluster(0),         fhTimeTriggerEMCALBCBadCluster(0),
  88 fhEtaPhiTriggerEMCALBCUMBadCluster(0),       fhTimeTriggerEMCALBCUMBadCluster(0),
  89 fhEtaPhiTriggerEMCALBCBadExoticCluster(0),   fhTimeTriggerEMCALBCBadExoticCluster(0),
  90 fhEtaPhiTriggerEMCALBCUMBadExoticCluster(0), fhTimeTriggerEMCALBCUMBadExoticCluster(0),
  91 fhTimeTriggerEMCALBCBadMaxCell(0),           fhTimeTriggerEMCALBCUMBadMaxCell(0),
  92 fhTimeTriggerEMCALBCBadMaxCellExotic(0),     fhTimeTriggerEMCALBCUMBadMaxCellExotic(0),
  93 fhEtaPhiTriggerEMCALBCUMReMatchOpenTimeCluster (0), fhTimeTriggerEMCALBCUMReMatchOpenTimeCluster(0),
  94 fhEtaPhiTriggerEMCALBCUMReMatchCheckNeighCluster(0),fhTimeTriggerEMCALBCUMReMatchCheckNeighCluster(0),
  95 fhEtaPhiTriggerEMCALBCUMReMatchBothCluster(0),      fhTimeTriggerEMCALBCUMReMatchBothCluster(0),
  96 fhTimeTriggerEMCALBC0UMReMatchOpenTime(0),
  97 fhTimeTriggerEMCALBC0UMReMatchCheckNeigh(0),
  98 fhTimeTriggerEMCALBC0UMReMatchBoth(0),
  99
 100 fhEtaPhiNoTrigger(0),                        fhTimeNoTrigger(0),
 101
 102 fhEPhoton(0),                 fhPtPhoton(0),
 103 fhPhiPhoton(0),               fhEtaPhoton(0),
 104 fhEtaPhiPhoton(0),            fhEtaPhi05Photon(0),
 105 fhEtaPhiPhotonEMCALBC0(0),    fhEtaPhiPhotonEMCALBC1(0),   fhEtaPhiPhotonEMCALBCN(0),
 106 fhTimePhotonTriggerEMCALBC0UMReMatchOpenTime(0),
 107 fhTimePhotonTriggerEMCALBC0UMReMatchCheckNeigh(0),
 108 fhTimePhotonTriggerEMCALBC0UMReMatchBoth(0),
 109
 110 fhPtCentralityPhoton(0),      fhPtEventPlanePhoton(0),
 111
 112 // Shower shape histograms
 113 fhNLocMax(0),
 114 fhDispE(0),                   fhLam0E(0),                   fhLam1E(0),
 115 fhDispETRD(0),                fhLam0ETRD(0),                fhLam1ETRD(0),
 116 fhDispETM(0),                 fhLam0ETM(0),                 fhLam1ETM(0),
 117 fhDispETMTRD(0),              fhLam0ETMTRD(0),              fhLam1ETMTRD(0),
 118
 119 fhNCellsLam0LowE(0),          fhNCellsLam1LowE(0),          fhNCellsDispLowE(0),
 120 fhNCellsLam0HighE(0),         fhNCellsLam1HighE(0),         fhNCellsDispHighE(0),
 121
 122 fhEtaLam0LowE(0),             fhPhiLam0LowE(0),
 123 fhEtaLam0HighE(0),            fhPhiLam0HighE(0),
 124 fhLam0DispLowE(0),            fhLam0DispHighE(0),
 125 fhLam1Lam0LowE(0),            fhLam1Lam0HighE(0),
 126 fhDispLam1LowE(0),            fhDispLam1HighE(0),
 127 fhDispEtaE(0),                fhDispPhiE(0),
 128 fhSumEtaE(0),                 fhSumPhiE(0),                 fhSumEtaPhiE(0),
 129 fhDispEtaPhiDiffE(0),         fhSphericityE(0),
 130 fhDispSumEtaDiffE(0),         fhDispSumPhiDiffE(0),
 131
 132 // MC histograms
 133 fhMCPhotonELambda0NoOverlap(0),       fhMCPhotonELambda0TwoOverlap(0),      fhMCPhotonELambda0NOverlap(0),
 134 // Embedding
 135 fhEmbeddedSignalFractionEnergy(0),
 136 fhEmbedPhotonELambda0FullSignal(0),   fhEmbedPhotonELambda0MostlySignal(0),
 137 fhEmbedPhotonELambda0MostlyBkg(0),    fhEmbedPhotonELambda0FullBkg(0),
 138 fhEmbedPi0ELambda0FullSignal(0),      fhEmbedPi0ELambda0MostlySignal(0),
 139 fhEmbedPi0ELambda0MostlyBkg(0),       fhEmbedPi0ELambda0FullBkg(0),
 140 // PileUp
 141 fhTimePtNoCut(0),                     fhTimePtSPD(0),
 142 fhTimePtPhotonNoCut(0),               fhTimePtPhotonSPD(0),
 143 fhTimeNPileUpVertSPD(0),              fhTimeNPileUpVertTrack(0),
 144 fhTimeNPileUpVertContributors(0),
 145 fhTimePileUpMainVertexZDistance(0),   fhTimePileUpMainVertexZDiamond(0),
 146 fhClusterMultSPDPileUp(),             fhClusterMultNoPileUp(),
 147 fhEtaPhiBC0(0),  fhEtaPhiBCPlus(0),   fhEtaPhiBCMinus(0),
 148 fhEtaPhiBC0PileUpSPD(0),
 149 fhEtaPhiBCPlusPileUpSPD(0),           fhEtaPhiBCMinusPileUpSPD(0),
 150 fhPtNPileUpSPDVtx(0),                 fhPtNPileUpTrkVtx(0),
 151 fhPtNPileUpSPDVtxTimeCut(0),          fhPtNPileUpTrkVtxTimeCut(0),
 152 fhPtNPileUpSPDVtxTimeCut2(0),         fhPtNPileUpTrkVtxTimeCut2(0),
 153 fhPtPhotonNPileUpSPDVtx(0),           fhPtPhotonNPileUpTrkVtx(0),
 154 fhPtPhotonNPileUpSPDVtxTimeCut(0),    fhPtPhotonNPileUpTrkVtxTimeCut(0),
 155 fhPtPhotonNPileUpSPDVtxTimeCut2(0),   fhPtPhotonNPileUpTrkVtxTimeCut2(0),
 156 fhEClusterSM(0),                      fhEPhotonSM(0),
 157 fhPtClusterSM(0),                     fhPtPhotonSM(0)
 158 {
 159   //default ctor
 160
 161   for(Int_t i = 0; i < 14; i++)
 162   {
 163     fhMCPt     [i] = 0;
 164     fhMCE      [i] = 0;
 165     fhMCPhi    [i] = 0;
 166     fhMCEta    [i] = 0;
 167     fhMCDeltaE [i] = 0;
 168     fhMCDeltaPt[i] = 0;
 169     fhMC2E     [i] = 0;
 170     fhMC2Pt    [i] = 0;
 171   }
 172
 173   for(Int_t i = 0; i < 7; i++)
 174   {
 175     fhPtPrimMC [i] = 0;
 176     fhEPrimMC  [i] = 0;
 177     fhPhiPrimMC[i] = 0;
 178     fhEtaPrimMC[i] = 0;
 179     fhYPrimMC  [i] = 0;
 180
 181     fhPtPrimMCAcc [i] = 0;
 182     fhEPrimMCAcc  [i] = 0;
 183     fhPhiPrimMCAcc[i] = 0;
 184     fhEtaPrimMCAcc[i] = 0;
 185     fhYPrimMCAcc  [i] = 0;
 186
 187     fhDispEtaDispPhi[i] = 0;
 188     fhLambda0DispPhi[i] = 0;
 189     fhLambda0DispEta[i] = 0;
 190
 191     fhPtPileUp       [i] = 0;
 192     fhPtChargedPileUp[i] = 0;
 193     fhPtPhotonPileUp [i] = 0;
 194
 195     fhLambda0PileUp       [i] = 0;
 196     fhLambda0ChargedPileUp[i] = 0;
 197
 198     fhClusterEFracLongTimePileUp  [i] = 0;
 199
 200     fhClusterCellTimePileUp       [i] = 0;
 201     fhClusterTimeDiffPileUp       [i] = 0;
 202     fhClusterTimeDiffChargedPileUp[i] = 0;
 203     fhClusterTimeDiffPhotonPileUp [i] = 0;
 204
 205     for(Int_t j = 0; j < 6; j++)
 206     {
 207       fhMCDispEtaDispPhi[i][j] = 0;
 208       fhMCLambda0DispEta[i][j] = 0;
 209       fhMCLambda0DispPhi[i][j] = 0;
 210     }
 211   }
 212
 213   for(Int_t i = 0; i < 6; i++)
 214   {
 215     fhMCELambda0    [i]                  = 0;
 216     fhMCELambda1    [i]                  = 0;
 217     fhMCEDispersion [i]                  = 0;
 218     fhMCNCellsE     [i]                  = 0;
 219     fhMCMaxCellDiffClusterE[i]           = 0;
 220     fhLambda0DispEta[i]                  = 0;
 221     fhLambda0DispPhi[i]                  = 0;
 222
 223     fhMCLambda0vsClusterMaxCellDiffE0[i] = 0;
 224     fhMCLambda0vsClusterMaxCellDiffE2[i] = 0;
 225     fhMCLambda0vsClusterMaxCellDiffE6[i] = 0;
 226     fhMCNCellsvsClusterMaxCellDiffE0 [i] = 0;
 227     fhMCNCellsvsClusterMaxCellDiffE2 [i] = 0;
 228     fhMCNCellsvsClusterMaxCellDiffE6 [i] = 0;
 229
 230     fhMCEDispEta       [i]               = 0;
 231     fhMCEDispPhi       [i]               = 0;
 232     fhMCESumEtaPhi     [i]               = 0;
 233     fhMCEDispEtaPhiDiff[i]               = 0;
 234     fhMCESphericity    [i]               = 0;
 235   }
 236
 237   for(Int_t i = 0; i < 5; i++)
 238   {
 239     fhClusterCuts[i] = 0;
 240   }
 241
 242   // Track matching residuals
 243   for(Int_t i = 0; i < 2; i++)
 244   {
 245     fhTrackMatchedDEta   [i] = 0;             fhTrackMatchedDPhi   [i] = 0;         fhTrackMatchedDEtaDPhi   [i] = 0;
 246     fhTrackMatchedDEtaNeg[i] = 0;             fhTrackMatchedDPhiNeg[i] = 0;         fhTrackMatchedDEtaDPhiNeg[i] = 0;
 247     fhTrackMatchedDEtaPos[i] = 0;             fhTrackMatchedDPhiPos[i] = 0;         fhTrackMatchedDEtaDPhiPos[i] = 0;
 248     fhTrackMatchedDEtaTRD[i] = 0;             fhTrackMatchedDPhiTRD[i] = 0;
 249     fhTrackMatchedDEtaMCOverlap[i] = 0;       fhTrackMatchedDPhiMCOverlap[i] = 0;
 250     fhTrackMatchedDEtaMCNoOverlap[i] = 0;     fhTrackMatchedDPhiMCNoOverlap[i] = 0;
 251     fhTrackMatchedDEtaMCConversion[i] = 0;    fhTrackMatchedDPhiMCConversion[i] = 0;
 252     fhTrackMatchedMCParticle[i] = 0;          fhTrackMatchedMCParticle[i] = 0;
 253     fhdEdx[i] = 0;                            fhEOverP[i] = 0;
 254     fhEOverPTRD[i] = 0;
 255   }
 256
 257   for(Int_t i = 0; i < 4; i++)
 258   {
 259     fhClusterMultSPDPileUp[i] = 0;
 260     fhClusterMultNoPileUp [i] = 0;
 261   }
 262
 263   for(Int_t i = 0; i < 11; i++)
 264   {
 265     fhEtaPhiTriggerEMCALBC             [i] = 0 ;
 266     fhTimeTriggerEMCALBC               [i] = 0 ;
 267     fhEtaPhiTriggerEMCALBCUM           [i] = 0 ;
 268     fhTimeTriggerEMCALBCUM             [i] = 0 ;
 269
 270     fhEtaPhiPhotonTriggerEMCALBC       [i] = 0 ;
 271     fhTimePhotonTriggerEMCALBC         [i] = 0 ;
 272     fhEtaPhiPhotonTriggerEMCALBCUM     [i] = 0 ;
 273     fhTimePhotonTriggerEMCALBCUM       [i] = 0 ;
 274
 275     fhTimePhotonTriggerEMCALBCPileUpSPD[i] = 0 ;
 276     fhTimeTriggerEMCALBCPileUpSPD      [i] = 0 ;
 277
 278     fhEtaPhiTriggerEMCALBCCluster      [i] = 0 ;
 279     fhEtaPhiTriggerEMCALBCUMCluster    [i] = 0 ;
 280   }
 281
 282   //Initialize parameters
 283   InitParameters();
 284
 285 }
 286
 287 //_________________________________________________________________________________________
 288 Bool_t  AliAnaPhoton::ClusterSelected(AliVCluster* calo, TLorentzVector mom, Int_t nMaxima)
 289 {
 290   //Select clusters if they pass different cuts
 291
 292   Float_t ptcluster  = mom.Pt();
 293   Float_t ecluster   = mom.E();
 294   Float_t l0cluster  = calo->GetM02();
 295   Float_t etacluster = mom.Eta();
 296   Float_t phicluster = mom.Phi();
 297   if(phicluster < 0) phicluster+=TMath::TwoPi();
 298   Float_t tofcluster   = calo->GetTOF()*1.e9;
 299
 300   Bool_t matched = IsTrackMatched(calo,GetReader()->GetInputEvent());
 301
 302   if(GetDebug() > 2)
 303     printf("AliAnaPhoton::ClusterSelected() - Current Event %d; Before selection : E %2.2f, pT %2.2f, phi %2.2f, eta %2.2f\n",
 304            GetReader()->GetEventNumber(),
 305            ecluster,ptcluster, phicluster*TMath::RadToDeg(),etacluster);
 306
 307   fhClusterCuts[1]->Fill(ecluster);
 308
 309   if(ecluster > 0.5) fhEtaPhi->Fill(etacluster, phicluster);
 310
 311   Int_t   nSM  = GetModuleNumber(calo);
 312   if(nSM < GetCaloUtils()->GetNumberOfSuperModulesUsed() && nSM >=0)
 313   {
 314     fhEClusterSM ->Fill(ecluster ,nSM);
 315     fhPtClusterSM->Fill(ptcluster,nSM);
 316   }
 317
 318   FillEMCALTriggerClusterBCHistograms(calo->GetID(),ecluster,tofcluster,etacluster,phicluster);
 319
 320   //.......................................
 321   //If too small or big energy, skip it
 322   if(ecluster < GetMinEnergy() || ecluster > GetMaxEnergy() ) return kFALSE ;
 323
 324   if(GetDebug() > 2) printf("\t Cluster %d Pass E Cut \n",calo->GetID());
 325
 326   fhClusterCuts[2]->Fill(ecluster);
 327
 328   FillClusterPileUpHistograms(calo,matched,ptcluster,etacluster,phicluster,l0cluster);
 329
 330   //.......................................
 331   // TOF cut, BE CAREFUL WITH THIS CUT
 332   Double_t tof = calo->GetTOF()*1e9;
 333   if(tof < fTimeCutMin || tof > fTimeCutMax) return kFALSE;
 334
 335   if(GetDebug() > 2)  printf("\t Cluster %d Pass Time Cut \n",calo->GetID());
 336
 337   fhClusterCuts[3]->Fill(ecluster);
 338
 339   //.......................................
 340   if(calo->GetNCells() <= fNCellsCut && GetReader()->GetDataType() != AliCaloTrackReader::kMC) return kFALSE;
 341
 342   if(GetDebug() > 2) printf("\t Cluster %d Pass NCell Cut \n",calo->GetID());
 343
 344   fhClusterCuts[4]->Fill(ecluster);
 345
 346   if(nMaxima < fNLMCutMin || nMaxima > fNLMCutMax) return kFALSE ;
 347   if(GetDebug() > 2) printf(" \t Cluster %d pass NLM %d of out of range \n",calo->GetID(), nMaxima);
 348
 349   fhClusterCuts[5]->Fill(ecluster);
 350
 351   //.......................................
 352   //Check acceptance selection
 353   if(IsFiducialCutOn())
 354   {
 355     Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter) ;
 356     if(! in ) return kFALSE ;
 357   }
 358
 359   if(GetDebug() > 2) printf("\t Fiducial cut passed \n");
 360
 361   fhClusterCuts[6]->Fill(ecluster);
 362
 363   //.......................................
 364   //Skip matched clusters with tracks
 365
 366   // Fill matching residual histograms before PID cuts
 367   if(fFillTMHisto) FillTrackMatchingResidualHistograms(calo,0);
 368
 369   if(fRejectTrackMatch)
 370   {
 371     if(matched)
 372     {
 373       if(GetDebug() > 2) printf("\t Reject track-matched clusters\n");
 374       return kFALSE ;
 375     }
 376     else
 377       if(GetDebug() > 2)  printf(" Track-matching cut passed \n");
 378   }// reject matched clusters
 379
 380   fhClusterCuts[7]->Fill(ecluster);
 381
 382   if(fFillPileUpHistograms)
 383   {
 384     if(GetReader()->IsPileUpFromSPD())               {fhPtChargedPileUp[0]->Fill(ptcluster); fhLambda0ChargedPileUp[0]->Fill(ecluster,l0cluster); }
 385     if(GetReader()->IsPileUpFromEMCal())             {fhPtChargedPileUp[1]->Fill(ptcluster); fhLambda0ChargedPileUp[1]->Fill(ecluster,l0cluster); }
 386     if(GetReader()->IsPileUpFromSPDOrEMCal())        {fhPtChargedPileUp[2]->Fill(ptcluster); fhLambda0ChargedPileUp[2]->Fill(ecluster,l0cluster); }
 387     if(GetReader()->IsPileUpFromSPDAndEMCal())       {fhPtChargedPileUp[3]->Fill(ptcluster); fhLambda0ChargedPileUp[3]->Fill(ecluster,l0cluster); }
 388     if(GetReader()->IsPileUpFromSPDAndNotEMCal())    {fhPtChargedPileUp[4]->Fill(ptcluster); fhLambda0ChargedPileUp[4]->Fill(ecluster,l0cluster); }
 389     if(GetReader()->IsPileUpFromEMCalAndNotSPD())    {fhPtChargedPileUp[5]->Fill(ptcluster); fhLambda0ChargedPileUp[5]->Fill(ecluster,l0cluster); }
 390     if(GetReader()->IsPileUpFromNotSPDAndNotEMCal()) {fhPtChargedPileUp[6]->Fill(ptcluster); fhLambda0ChargedPileUp[6]->Fill(ecluster,l0cluster); }
 391   }
 392
 393   //.......................................
 394   //Check Distance to Bad channel, set bit.
 395   Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
 396   if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;
 397   if(distBad < fMinDist)
 398   {//In bad channel (PHOS cristal size 2.2x2.2 cm), EMCAL ( cell units )
 399     return kFALSE ;
 400   }
 401   else if(GetDebug() > 2) printf("\t Bad channel cut passed %4.2f > %2.2f \n",distBad, fMinDist);
 402
 403   fhClusterCuts[8]->Fill(ecluster);
 404
 405   if(GetDebug() > 0)
 406     printf("AliAnaPhoton::ClusterSelected() Current Event %d; After  selection : E %2.2f, pT %2.2f, phi %2.2f, eta %2.2f\n",
 407            GetReader()->GetEventNumber(),
 408            ecluster, ptcluster,mom.Phi()*TMath::RadToDeg(),mom.Eta());
 409
 410   //All checks passed, cluster selected
 411   return kTRUE;
 412
 413 }
 414
 415 //___________________________________________
 416 void AliAnaPhoton::FillAcceptanceHistograms()
 417 {
 418   //Fill acceptance histograms if MC data is available
 419
 420   Double_t photonY   = -100 ;
 421   Double_t photonE   = -1 ;
 422   Double_t photonPt  = -1 ;
 423   Double_t photonPhi =  100 ;
 424   Double_t photonEta = -1 ;
 425
 426   Int_t    pdg       =  0 ;
 427   Int_t    tag       =  0 ;
 428   Int_t    mcIndex   =  0 ;
 429   Bool_t   inacceptance = kFALSE;
 430
 431   if(GetReader()->ReadStack())
 432   {
 433     AliStack * stack = GetMCStack();
 434     if(stack)
 435     {
 436       for(Int_t i=0 ; i<stack->GetNtrack(); i++)
 437       {
 438         if(GetReader()->AcceptOnlyHIJINGLabels() && !GetReader()->IsHIJINGLabel(i)) continue ;
 439
 440         TParticle * prim = stack->Particle(i) ;
 441         pdg = prim->GetPdgCode();
 442         //printf("i %d, %s %d  %s %d \n",i, stack->Particle(i)->GetName(), stack->Particle(i)->GetPdgCode(),
 443         //                             prim->GetName(), prim->GetPdgCode());
 444
 445         if(pdg == 22)
 446         {
 447           // Get tag of this particle photon from fragmentation, decay, prompt ...
 448           tag = GetMCAnalysisUtils()->CheckOrigin(i,GetReader());
 449           if(!GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton))
 450           {
 451             //A conversion photon from a hadron, skip this kind of photon
 452             // printf("AliAnaPhoton::FillAcceptanceHistograms() - not a photon, weird!\n ");
 453             // GetMCAnalysisUtils()->PrintMCTag(tag);
 454
 455             return;
 456           }
 457
 458           //Get photon kinematics
 459           if(prim->Energy() == TMath::Abs(prim->Pz()))  continue ; //Protection against floating point exception
 460
 461           photonY   = 0.5*TMath::Log((prim->Energy()-prim->Pz())/(prim->Energy()+prim->Pz())) ;
 462           photonE   = prim->Energy() ;
 463           photonPt  = prim->Pt() ;
 464           photonPhi = TMath::RadToDeg()*prim->Phi() ;
 465           if(photonPhi < 0) photonPhi+=TMath::TwoPi();
 466           photonEta = prim->Eta() ;
 467
 468           //Check if photons hit the Calorimeter
 469           TLorentzVector lv;
 470           prim->Momentum(lv);
 471           inacceptance = kFALSE;
 472           if     (fCalorimeter == "PHOS")
 473           {
 474             if(GetPHOSGeometry() && GetCaloUtils()->IsPHOSGeoMatrixSet())
 475             {
 476               Int_t mod ;
 477               Double_t x,z ;
 478               if(GetPHOSGeometry()->ImpactOnEmc(prim,mod,z,x))
 479                 inacceptance = kTRUE;
 480               if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
 481             }
 482             else
 483             {
 484               if(GetFiducialCut()->IsInFiducialCut(lv,fCalorimeter))
 485                 inacceptance = kTRUE ;
 486               if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
 487             }
 488           }
 489           else if(fCalorimeter == "EMCAL" && GetCaloUtils()->IsEMCALGeoMatrixSet())
 490           {
 491             if(GetEMCALGeometry())
 492             {
 493               Int_t absID=0;
 494
 495               GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(prim->Eta(),prim->Phi(),absID);
 496
 497               if( absID >= 0)
 498                 inacceptance = kTRUE;
 499
 500               //                  if(GetEMCALGeometry()->Impact(phot1) && GetEMCALGeometry()->Impact(phot2))
 501               //                    inacceptance = kTRUE;
 502               if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
 503             }
 504             else
 505             {
 506               if(GetFiducialCut()->IsInFiducialCut(lv,fCalorimeter))
 507                 inacceptance = kTRUE ;
 508               if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
 509             }
 510           }       //In EMCAL
 511
 512           //Fill histograms
 513           fhYPrimMC[kmcPPhoton]->Fill(photonPt, photonY) ;
 514           if(TMath::Abs(photonY) < 1.0)
 515           {
 516             fhEPrimMC  [kmcPPhoton]->Fill(photonE ) ;
 517             fhPtPrimMC [kmcPPhoton]->Fill(photonPt) ;
 518             fhPhiPrimMC[kmcPPhoton]->Fill(photonE , photonPhi) ;
 519             fhEtaPrimMC[kmcPPhoton]->Fill(photonE , photonEta) ;
 520           }
 521           if(inacceptance)
 522           {
 523             fhEPrimMCAcc  [kmcPPhoton]->Fill(photonE ) ;
 524             fhPtPrimMCAcc [kmcPPhoton]->Fill(photonPt) ;
 525             fhPhiPrimMCAcc[kmcPPhoton]->Fill(photonE , photonPhi) ;
 526             fhEtaPrimMCAcc[kmcPPhoton]->Fill(photonE , photonEta) ;
 527             fhYPrimMCAcc  [kmcPPhoton]->Fill(photonE , photonY) ;
 528           }//Accepted
 529
 530           //Origin of photon
 531           if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt) && fhEPrimMC[kmcPPrompt])
 532           {
 533             mcIndex = kmcPPrompt;
 534           }
 535           else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation) && fhEPrimMC[kmcPFragmentation])
 536           {
 537             mcIndex = kmcPFragmentation ;
 538           }
 539           else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCISR) && fhEPrimMC[kmcPISR])
 540           {
 541             mcIndex = kmcPISR;
 542           }
 543           else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay)&& fhEPrimMC[kmcPPi0Decay])
 544           {
 545             mcIndex = kmcPPi0Decay;
 546           }
 547           else if( (GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay) ||
 548                     GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay)) && fhEPrimMC[kmcPOtherDecay])
 549           {
 550             mcIndex = kmcPOtherDecay;
 551           }
 552           else if(fhEPrimMC[kmcPOther])
 553           {
 554             mcIndex = kmcPOther;
 555           }//Other origin
 556
 557           fhYPrimMC[mcIndex]->Fill(photonPt, photonY) ;
 558           if(TMath::Abs(photonY) < 1.0)
 559           {
 560             fhEPrimMC  [mcIndex]->Fill(photonE ) ;
 561             fhPtPrimMC [mcIndex]->Fill(photonPt) ;
 562             fhPhiPrimMC[mcIndex]->Fill(photonE , photonPhi) ;
 563             fhEtaPrimMC[mcIndex]->Fill(photonE , photonEta) ;
 564           }
 565
 566           if(inacceptance)
 567           {
 568             fhEPrimMCAcc  [mcIndex]->Fill(photonE ) ;
 569             fhPtPrimMCAcc [mcIndex]->Fill(photonPt) ;
 570             fhPhiPrimMCAcc[mcIndex]->Fill(photonE , photonPhi) ;
 571             fhEtaPrimMCAcc[mcIndex]->Fill(photonE , photonEta) ;
 572             fhYPrimMCAcc  [mcIndex]->Fill(photonE , photonY) ;
 573           }//Accepted
 574
 575         }// Primary photon
 576       }//loop on primaries
 577     }//stack exists and data is MC
 578   }//read stack
 579   else if(GetReader()->ReadAODMCParticles())
 580   {
 581     TClonesArray * mcparticles = GetReader()->GetAODMCParticles();
 582     if(mcparticles)
 583     {
 584       Int_t nprim = mcparticles->GetEntriesFast();
 585
 586       for(Int_t i=0; i < nprim; i++)
 587       {
 588         if(GetReader()->AcceptOnlyHIJINGLabels() && !GetReader()->IsHIJINGLabel(i)) continue ;
 589
 590         AliAODMCParticle * prim = (AliAODMCParticle *) mcparticles->At(i);
 591
 592         pdg = prim->GetPdgCode();
 593
 594         if(pdg == 22)
 595         {
 596           // Get tag of this particle photon from fragmentation, decay, prompt ...
 597           tag = GetMCAnalysisUtils()->CheckOrigin(i,GetReader());
 598           if(!GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton))
 599           {
 600             //A conversion photon from a hadron, skip this kind of photon
 601             //            printf("AliAnaPhoton::FillAcceptanceHistograms() - not a photon, weird!\n ");
 602             //            GetMCAnalysisUtils()->PrintMCTag(tag);
 603
 604             return;
 605           }
 606
 607           //Get photon kinematics
 608           if(prim->E() == TMath::Abs(prim->Pz()))  continue ; //Protection against floating point exception
 609
 610           photonY   = 0.5*TMath::Log((prim->E()-prim->Pz())/(prim->E()+prim->Pz())) ;
 611           photonE   = prim->E() ;
 612           photonPt  = prim->Pt() ;
 613           photonPhi = prim->Phi() ;
 614           if(photonPhi < 0) photonPhi+=TMath::TwoPi();
 615           photonEta = prim->Eta() ;
 616
 617           //Check if photons hit the Calorimeter
 618           TLorentzVector lv;
 619           lv.SetPxPyPzE(prim->Px(),prim->Py(),prim->Pz(),prim->E());
 620           inacceptance = kFALSE;
 621           if     (fCalorimeter == "PHOS")
 622           {
 623             if(GetPHOSGeometry() && GetCaloUtils()->IsPHOSGeoMatrixSet())
 624             {
 625               Int_t mod ;
 626               Double_t x,z ;
 627               Double_t vtx[]={prim->Xv(),prim->Yv(),prim->Zv()};
 628               if(GetPHOSGeometry()->ImpactOnEmc(vtx, prim->Theta(),prim->Phi(),mod,z,x))
 629                 inacceptance = kTRUE;
 630               if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
 631             }
 632             else
 633             {
 634               if(GetFiducialCut()->IsInFiducialCut(lv,fCalorimeter))
 635                 inacceptance = kTRUE ;
 636               if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
 637             }
 638           }
 639           else if(fCalorimeter == "EMCAL" && GetCaloUtils()->IsEMCALGeoMatrixSet())
 640           {
 641             if(GetEMCALGeometry())
 642             {
 643               Int_t absID=0;
 644
 645               GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(prim->Eta(),prim->Phi(),absID);
 646
 647               if( absID >= 0)
 648                 inacceptance = kTRUE;
 649
 650               if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
 651             }
 652             else
 653             {
 654               if(GetFiducialCut()->IsInFiducialCut(lv,fCalorimeter))
 655                 inacceptance = kTRUE ;
 656               if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
 657             }
 658           }       //In EMCAL
 659
 660           //Fill histograms
 661
 662           fhYPrimMC[kmcPPhoton]->Fill(photonPt, photonY) ;
 663           if(TMath::Abs(photonY) < 1.0)
 664           {
 665             fhEPrimMC  [kmcPPhoton]->Fill(photonE ) ;
 666             fhPtPrimMC [kmcPPhoton]->Fill(photonPt) ;
 667             fhPhiPrimMC[kmcPPhoton]->Fill(photonE , photonPhi) ;
 668             fhEtaPrimMC[kmcPPhoton]->Fill(photonE , photonEta) ;
 669           }
 670
 671           if(inacceptance)
 672           {
 673             fhEPrimMCAcc[kmcPPhoton]  ->Fill(photonE ) ;
 674             fhPtPrimMCAcc[kmcPPhoton] ->Fill(photonPt) ;
 675             fhPhiPrimMCAcc[kmcPPhoton]->Fill(photonE , photonPhi) ;
 676             fhEtaPrimMCAcc[kmcPPhoton]->Fill(photonE , photonEta) ;
 677             fhYPrimMCAcc[kmcPPhoton]  ->Fill(photonE , photonY) ;
 678           }//Accepted
 679
 680           //Origin of photon
 681           if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt) && fhEPrimMC[kmcPPrompt])
 682           {
 683             mcIndex = kmcPPrompt;
 684           }
 685           else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation) && fhEPrimMC[kmcPFragmentation])
 686           {
 687             mcIndex = kmcPFragmentation ;
 688           }
 689           else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCISR) && fhEPrimMC[kmcPISR])
 690           {
 691             mcIndex = kmcPISR;
 692           }
 693           else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay)&& fhEPrimMC[kmcPPi0Decay])
 694           {
 695             mcIndex = kmcPPi0Decay;
 696           }
 697           else if( (GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay) ||
 698                     GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay)) && fhEPrimMC[kmcPOtherDecay])
 699           {
 700             mcIndex = kmcPOtherDecay;
 701           }
 702           else if(fhEPrimMC[kmcPOther])
 703           {
 704             mcIndex = kmcPOther;
 705           }//Other origin
 706
 707           fhYPrimMC[mcIndex]->Fill(photonPt, photonY) ;
 708           if(TMath::Abs(photonY) < 1.0)
 709           {
 710             fhEPrimMC  [mcIndex]->Fill(photonE ) ;
 711             fhPtPrimMC [mcIndex]->Fill(photonPt) ;
 712             fhPhiPrimMC[mcIndex]->Fill(photonE , photonPhi) ;
 713             fhEtaPrimMC[mcIndex]->Fill(photonE , photonEta) ;
 714           }
 715           if(inacceptance)
 716           {
 717             fhEPrimMCAcc  [mcIndex]->Fill(photonE ) ;
 718             fhPtPrimMCAcc [mcIndex]->Fill(photonPt) ;
 719             fhPhiPrimMCAcc[mcIndex]->Fill(photonE , photonPhi) ;
 720             fhEtaPrimMCAcc[mcIndex]->Fill(photonE , photonEta) ;
 721             fhYPrimMCAcc  [mcIndex]->Fill(photonE , photonY) ;
 722           }//Accepted
 723
 724         }// Primary photon
 725       }//loop on primaries
 726
 727     }//kmc array exists and data is MC
 728   }     // read AOD MC
 729 }
 730
 731 //________________________________________________________________________________________________________________
 732 void  AliAnaPhoton::FillEMCALTriggerClusterBCHistograms(Int_t idcalo,       Float_t ecluster,  Float_t tofcluster,
 733                                                         Float_t etacluster, Float_t phicluster)
 734
 735 {
 736   // Fill trigger related histograms
 737
 738   if(!fFillEMCALBCHistograms || fCalorimeter!="EMCAL") return ;
 739
 740   Float_t tofclusterUS = TMath::Abs(tofcluster);
 741
 742   if(ecluster > 2)
 743   {
 744     if      (tofclusterUS < 25) fhEtaPhiEMCALBC0->Fill(etacluster, phicluster);
 745     else if (tofclusterUS < 75) fhEtaPhiEMCALBC1->Fill(etacluster, phicluster);
 746     else                        fhEtaPhiEMCALBCN->Fill(etacluster, phicluster);
 747   }
 748
 749   Int_t  bc     = GetReader()->GetTriggerClusterBC();
 750   Int_t  id     = GetReader()->GetTriggerClusterId();
 751   Bool_t badMax = GetReader()->IsBadMaxCellTriggerEvent();
 752
 753   Int_t histoBC = bc+5;
 754   if(GetReader()->AreBadTriggerEventsRemoved()) histoBC=0; // histograms created only for one BC since the others where rejected
 755
 756   if(id==-2)
 757   {
 758     //printf("AliAnaPhoton::ClusterSelected() - No trigger found bc=%d\n",bc);
 759     fhEtaPhiNoTrigger->Fill(etacluster, phicluster);
 760     fhTimeNoTrigger  ->Fill(ecluster, tofcluster);
 761   }
 762   else if(TMath::Abs(bc) < 6)
 763   {
 764     if(!GetReader()->IsBadCellTriggerEvent() && !GetReader()->IsExoticEvent() )
 765     {
 766       if(GetReader()->IsTriggerMatched())
 767       {
 768         if(ecluster > 2) fhEtaPhiTriggerEMCALBC[histoBC]->Fill(etacluster, phicluster);
 769         fhTimeTriggerEMCALBC[histoBC]->Fill(ecluster, tofcluster);
 770         if(GetReader()->IsPileUpFromSPD()) fhTimeTriggerEMCALBCPileUpSPD[histoBC]->Fill(ecluster, tofcluster);
 771
 772         if(idcalo ==  GetReader()->GetTriggerClusterId())
 773         {
 774           fhEtaPhiTriggerEMCALBCCluster[histoBC]->Fill(etacluster, phicluster);
 775           fhTimeTriggerEMCALBCCluster        ->Fill(ecluster, tofcluster);
 776
 777           if(bc==0)
 778           {
 779             Float_t threshold = GetReader()->GetEventTriggerL1Threshold() ;
 780             if(GetReader()->IsEventEMCALL0()) threshold = GetReader()->GetEventTriggerL0Threshold() ;
 781
 782             if(ecluster > threshold)
 783               fhEtaPhiTriggerEMCALBCClusterOverTh->Fill(etacluster, phicluster);
 784             else if(ecluster > threshold-1)
 785               fhEtaPhiTriggerEMCALBCClusterBelowTh1->Fill(etacluster, phicluster);
 786             else
 787               fhEtaPhiTriggerEMCALBCClusterBelowTh2->Fill(etacluster, phicluster);
 788           }
 789         }
 790       }
 791       else
 792       {
 793         if(ecluster > 2) fhEtaPhiTriggerEMCALBCUM[histoBC]->Fill(etacluster, phicluster);
 794         fhTimeTriggerEMCALBCUM[histoBC]->Fill(ecluster, tofcluster);
 795
 796         if(bc==0)
 797         {
 798           if(GetReader()->IsTriggerMatchedOpenCuts(0)) fhTimeTriggerEMCALBC0UMReMatchOpenTime   ->Fill(ecluster, tofcluster);
 799           if(GetReader()->IsTriggerMatchedOpenCuts(1)) fhTimeTriggerEMCALBC0UMReMatchCheckNeigh ->Fill(ecluster, tofcluster);
 800           if(GetReader()->IsTriggerMatchedOpenCuts(2)) fhTimeTriggerEMCALBC0UMReMatchBoth       ->Fill(ecluster, tofcluster);
 801         }
 802
 803         if(idcalo ==  GetReader()->GetTriggerClusterId())
 804         {
 805           fhEtaPhiTriggerEMCALBCUMCluster[histoBC]->Fill(etacluster, phicluster);
 806           fhTimeTriggerEMCALBCUMCluster->Fill(ecluster, tofcluster);
 807           if(bc==0)
 808           {
 809             Float_t threshold = GetReader()->GetEventTriggerL1Threshold() ;
 810             if(GetReader()->IsEventEMCALL0()) threshold = GetReader()->GetEventTriggerL0Threshold() ;
 811
 812             if(ecluster > threshold)
 813               fhEtaPhiTriggerEMCALBCUMClusterOverTh->Fill(etacluster, phicluster);
 814             else if(ecluster > threshold-1)
 815               fhEtaPhiTriggerEMCALBCUMClusterBelowTh1->Fill(etacluster, phicluster);
 816             else
 817               fhEtaPhiTriggerEMCALBCUMClusterBelowTh2->Fill(etacluster, phicluster);
 818
 819             if(GetReader()->IsTriggerMatchedOpenCuts(0))
 820             {
 821               fhEtaPhiTriggerEMCALBCUMReMatchOpenTimeCluster->Fill(etacluster, phicluster);
 822               fhTimeTriggerEMCALBCUMReMatchOpenTimeCluster  ->Fill(ecluster, tofcluster);
 823             }
 824             if(GetReader()->IsTriggerMatchedOpenCuts(1))
 825             {
 826               fhEtaPhiTriggerEMCALBCUMReMatchCheckNeighCluster->Fill(etacluster, phicluster);
 827               fhTimeTriggerEMCALBCUMReMatchCheckNeighCluster  ->Fill(ecluster, tofcluster);
 828             }
 829             if(GetReader()->IsTriggerMatchedOpenCuts(2))
 830             {
 831               fhEtaPhiTriggerEMCALBCUMReMatchBothCluster->Fill(etacluster, phicluster);
 832               fhTimeTriggerEMCALBCUMReMatchBothCluster  ->Fill(ecluster, tofcluster);
 833             }
 834
 835           }
 836         }
 837       }
 838     }// neither bad nor exotic
 839     else if(GetReader()->IsBadCellTriggerEvent() && GetReader()->IsExoticEvent())
 840     {
 841       if(GetReader()->IsTriggerMatched())
 842       {
 843         if(ecluster > 2) fhEtaPhiTriggerEMCALBCBadExotic->Fill(etacluster, phicluster);
 844         fhTimeTriggerEMCALBCBadExotic->Fill(ecluster, tofcluster);
 845         if(badMax)  fhTimeTriggerEMCALBCBadMaxCellExotic->Fill(ecluster, tofcluster);
 846       }
 847       else
 848       {
 849         if(ecluster > 2) fhEtaPhiTriggerEMCALBCUMBadExotic->Fill(etacluster, phicluster);
 850         fhTimeTriggerEMCALBCUMBadExotic->Fill(ecluster, tofcluster);
 851         if(badMax)  fhTimeTriggerEMCALBCUMBadMaxCellExotic->Fill(ecluster, tofcluster);
 852
 853       }
 854     }// Bad and exotic cluster trigger
 855     else if(GetReader()->IsBadCellTriggerEvent() )
 856     {
 857       if(GetReader()->IsTriggerMatched())
 858       {
 859         if(ecluster > 2) fhEtaPhiTriggerEMCALBCBad->Fill(etacluster, phicluster);
 860         fhTimeTriggerEMCALBCBad->Fill(ecluster, tofcluster);
 861         if(badMax)  fhTimeTriggerEMCALBCBadMaxCell->Fill(ecluster, tofcluster);
 862       }
 863       else
 864       {
 865         if(ecluster > 2) fhEtaPhiTriggerEMCALBCUMBad->Fill(etacluster, phicluster);
 866         fhTimeTriggerEMCALBCUMBad->Fill(ecluster, tofcluster);
 867         if(badMax)  fhTimeTriggerEMCALBCUMBadMaxCell->Fill(ecluster, tofcluster);
 868       }
 869     }// Bad cluster trigger
 870     else if(GetReader()->IsExoticEvent() )
 871     {
 872       if(GetReader()->IsTriggerMatched())
 873       {
 874         if(ecluster > 2) fhEtaPhiTriggerEMCALBCExotic->Fill(etacluster, phicluster);
 875         fhTimeTriggerEMCALBCExotic->Fill(ecluster, tofcluster);
 876       }
 877       else
 878       {
 879         if(ecluster > 2) fhEtaPhiTriggerEMCALBCUMExotic->Fill(etacluster, phicluster);
 880         fhTimeTriggerEMCALBCUMExotic->Fill(ecluster, tofcluster);
 881       }
 882     }
 883   }
 884   else if(TMath::Abs(bc) >= 6)
 885     printf("AliAnaPhoton::ClusterSelected() - Trigger BC not expected = %d\n",bc);
 886
 887 }
 888
 889 //_________________________________________________________________________________________________________
 890 void  AliAnaPhoton::FillClusterPileUpHistograms(AliVCluster * calo, Bool_t matched,     Float_t ptcluster,
 891                                                 Float_t etacluster, Float_t phicluster, Float_t l0cluster)
 892 {
 893   // Fill some histograms related to pile up before any cluster cut is applied
 894
 895   if(!fFillPileUpHistograms) return ;
 896
 897   // Get the fraction of the cluster energy that carries the cell with highest energy and its absId
 898   AliVCaloCells* cells = 0;
 899   if(fCalorimeter == "EMCAL") cells = GetEMCALCells();
 900   else                        cells = GetPHOSCells();
 901
 902   Float_t maxCellFraction = 0.;
 903   Int_t absIdMax = GetCaloUtils()->GetMaxEnergyCell(cells, calo,maxCellFraction);
 904
 905   Double_t tmax  = cells->GetCellTime(absIdMax);
 906   GetCaloUtils()->RecalibrateCellTime(tmax, fCalorimeter, absIdMax,GetReader()->GetInputEvent()->GetBunchCrossNumber());
 907   tmax*=1.e9;
 908
 909   Bool_t okPhoton = kFALSE;
 910   if( GetCaloPID()->GetIdentifiedParticleType(calo)== AliCaloPID::kPhoton) okPhoton = kTRUE;
 911
 912   Float_t clusterLongTimePt = 0;
 913   Float_t clusterOKTimePt   = 0;
 914
 915   //Loop on cells inside cluster, max cell must be over 100 MeV and time in BC=0
 916   if(cells->GetCellAmplitude(absIdMax) > 0.1 && TMath::Abs(tmax) < 30)
 917   {
 918     for (Int_t ipos = 0; ipos < calo->GetNCells(); ipos++)
 919     {
 920       Int_t absId  = calo->GetCellsAbsId()[ipos];
 921
 922       if( absId == absIdMax ) continue ;
 923
 924       Double_t time  = cells->GetCellTime(absId);
 925       Float_t  amp   = cells->GetCellAmplitude(absId);
 926       Int_t    bc    = GetReader()->GetInputEvent()->GetBunchCrossNumber();
 927       GetCaloUtils()->GetEMCALRecoUtils()->AcceptCalibrateCell(absId,bc,amp,time,cells);
 928       time*=1e9;
 929
 930       Float_t diff = (tmax-time);
 931
 932       if(GetReader()->IsInTimeWindow(time,amp)) clusterOKTimePt   += amp;
 933       else                                      clusterLongTimePt += amp;
 934
 935       if( cells->GetCellAmplitude(absIdMax) < 0.1 ) continue ;
 936
 937       if(GetReader()->IsPileUpFromSPD())
 938       {
 939         fhClusterCellTimePileUp[0]->Fill(ptcluster, time);
 940         fhClusterTimeDiffPileUp[0]->Fill(ptcluster, diff);
 941         if(!matched)
 942         {
 943           fhClusterTimeDiffChargedPileUp[0]->Fill(ptcluster, diff);
 944           if(okPhoton)  fhClusterTimeDiffPhotonPileUp[0]->Fill(ptcluster, diff);
 945         }
 946       }
 947
 948       if(GetReader()->IsPileUpFromEMCal())
 949       {
 950         fhClusterCellTimePileUp[1]->Fill(ptcluster, time);
 951         fhClusterTimeDiffPileUp[1]->Fill(ptcluster, diff);
 952         if(!matched)
 953         {
 954           fhClusterTimeDiffChargedPileUp[1]->Fill(ptcluster, diff);
 955           if(okPhoton)  fhClusterTimeDiffPhotonPileUp[1]->Fill(ptcluster, diff);
 956         }
 957       }
 958
 959       if(GetReader()->IsPileUpFromSPDOrEMCal())
 960       {
 961         fhClusterCellTimePileUp[2]->Fill(ptcluster, time);
 962         fhClusterTimeDiffPileUp[2]->Fill(ptcluster, diff);
 963         if(!matched)
 964         {
 965           fhClusterTimeDiffChargedPileUp[2]->Fill(ptcluster, diff);
 966           if(okPhoton)  fhClusterTimeDiffPhotonPileUp[2]->Fill(ptcluster, diff);
 967         }
 968       }
 969
 970       if(GetReader()->IsPileUpFromSPDAndEMCal())
 971       {
 972         fhClusterCellTimePileUp[3]->Fill(ptcluster, time);
 973         fhClusterTimeDiffPileUp[3]->Fill(ptcluster, diff);
 974         if(!matched)
 975         {
 976           fhClusterTimeDiffChargedPileUp[3]->Fill(ptcluster, diff);
 977           if(okPhoton)  fhClusterTimeDiffPhotonPileUp[3]->Fill(ptcluster, diff);
 978         }
 979       }
 980
 981       if(GetReader()->IsPileUpFromSPDAndNotEMCal())
 982       {
 983         fhClusterCellTimePileUp[4]->Fill(ptcluster, time);
 984         fhClusterTimeDiffPileUp[4]->Fill(ptcluster, diff);
 985         if(!matched)
 986         {
 987           fhClusterTimeDiffChargedPileUp[4]->Fill(ptcluster, diff);
 988           if(okPhoton)  fhClusterTimeDiffPhotonPileUp[4]->Fill(ptcluster, diff);
 989         }
 990       }
 991
 992       if(GetReader()->IsPileUpFromEMCalAndNotSPD())
 993       {
 994         fhClusterCellTimePileUp[5]->Fill(ptcluster, time);
 995         fhClusterTimeDiffPileUp[5]->Fill(ptcluster, diff);
 996         if(!matched)
 997         {
 998           fhClusterTimeDiffChargedPileUp[5]->Fill(ptcluster, diff);
 999           if(okPhoton)  fhClusterTimeDiffPhotonPileUp[5]->Fill(ptcluster, diff);
1000         }
1001       }
1002
1003       if(GetReader()->IsPileUpFromNotSPDAndNotEMCal())
1004       {
1005         fhClusterCellTimePileUp[6]->Fill(ptcluster, time);
1006         fhClusterTimeDiffPileUp[6]->Fill(ptcluster, diff);
1007         if(!matched)
1008         {
1009           fhClusterTimeDiffChargedPileUp[6]->Fill(ptcluster, diff);
1010           if(okPhoton)  fhClusterTimeDiffPhotonPileUp[6]->Fill(ptcluster, diff);
1011         }
1012       }
1013     }//loop
1014
1015
1016     Float_t frac = 0;
1017     if(clusterLongTimePt+clusterOKTimePt > 0.001)
1018       frac = clusterLongTimePt/(clusterLongTimePt+clusterOKTimePt);
1019     //printf("E long %f, E OK %f, Fraction large time %f, E %f\n",clusterLongTimePt,clusterOKTimePt,frac,ptcluster);
1020
1021     if(GetReader()->IsPileUpFromSPD())               {fhPtPileUp[0]->Fill(ptcluster); fhLambda0PileUp[0]->Fill(ptcluster,l0cluster); fhClusterEFracLongTimePileUp[0]->Fill(ptcluster,frac);}
1022     if(GetReader()->IsPileUpFromEMCal())             {fhPtPileUp[1]->Fill(ptcluster); fhLambda0PileUp[1]->Fill(ptcluster,l0cluster); fhClusterEFracLongTimePileUp[1]->Fill(ptcluster,frac);}
1023     if(GetReader()->IsPileUpFromSPDOrEMCal())        {fhPtPileUp[2]->Fill(ptcluster); fhLambda0PileUp[2]->Fill(ptcluster,l0cluster); fhClusterEFracLongTimePileUp[2]->Fill(ptcluster,frac);}
1024     if(GetReader()->IsPileUpFromSPDAndEMCal())       {fhPtPileUp[3]->Fill(ptcluster); fhLambda0PileUp[3]->Fill(ptcluster,l0cluster); fhClusterEFracLongTimePileUp[3]->Fill(ptcluster,frac);}
1025     if(GetReader()->IsPileUpFromSPDAndNotEMCal())    {fhPtPileUp[4]->Fill(ptcluster); fhLambda0PileUp[4]->Fill(ptcluster,l0cluster); fhClusterEFracLongTimePileUp[4]->Fill(ptcluster,frac);}
1026     if(GetReader()->IsPileUpFromEMCalAndNotSPD())    {fhPtPileUp[5]->Fill(ptcluster); fhLambda0PileUp[5]->Fill(ptcluster,l0cluster); fhClusterEFracLongTimePileUp[5]->Fill(ptcluster,frac);}
1027     if(GetReader()->IsPileUpFromNotSPDAndNotEMCal()) {fhPtPileUp[6]->Fill(ptcluster); fhLambda0PileUp[6]->Fill(ptcluster,l0cluster); fhClusterEFracLongTimePileUp[6]->Fill(ptcluster,frac);}
1028
1029     fhEtaPhiBC0->Fill(etacluster,phicluster);
1030     if(GetReader()->IsPileUpFromSPD()) fhEtaPhiBC0PileUpSPD    ->Fill(etacluster,phicluster);
1031   }
1032
1033   else if (tmax > 25)         {fhEtaPhiBCPlus ->Fill(etacluster,phicluster); if(GetReader()->IsPileUpFromSPD()) fhEtaPhiBCPlusPileUpSPD ->Fill(etacluster,phicluster); }
1034   else if (tmax <-25)         {fhEtaPhiBCMinus->Fill(etacluster,phicluster); if(GetReader()->IsPileUpFromSPD()) fhEtaPhiBCMinusPileUpSPD->Fill(etacluster,phicluster); }
1035 }
1036
1037 //_______________________________________________
1038 void AliAnaPhoton::FillPileUpHistogramsPerEvent()
1039 {
1040   // Fill some histograms per event to understand pile-up
1041   // Open the time cut in the reader to be more meaningful
1042
1043   if(!fFillPileUpHistograms) return;
1044
1045   AliVEvent * event = GetReader()->GetInputEvent();
1046
1047   AliESDEvent* esdEv = dynamic_cast<AliESDEvent*> (event);
1048   AliAODEvent* aodEv = dynamic_cast<AliAODEvent*> (event);
1049
1050   // N pile up vertices
1051   Int_t nVtxSPD = -1;
1052   Int_t nVtxTrk = -1;
1053   TLorentzVector mom;
1054
1055   if      (esdEv)
1056   {
1057                 nVtxSPD = esdEv->GetNumberOfPileupVerticesSPD();
1058                 nVtxTrk = esdEv->GetNumberOfPileupVerticesTracks();
1059   }//ESD
1060   else if (aodEv)
1061   {
1062                 nVtxSPD = aodEv->GetNumberOfPileupVerticesSPD();
1063                 nVtxTrk = aodEv->GetNumberOfPileupVerticesTracks();
1064   }//AOD
1065
1066
1067         // Get the appropriate list of clusters
1068         TClonesArray * clusterList = 0;
1069         TString  clusterListName   = GetReader()->GetEMCALClusterListName();
1070         if     (event->FindListObject(clusterListName))
1071                 clusterList = dynamic_cast<TClonesArray*> (event->FindListObject(clusterListName));
1072         else if(GetReader()->GetOutputEvent())
1073                 clusterList = dynamic_cast<TClonesArray*> (GetReader()->GetOutputEvent()->FindListObject(clusterListName));
1074
1075         // Loop on clusters, get the maximum energy cluster as reference
1076   Int_t nclusters = 0;
1077         if(clusterList) nclusters = clusterList->GetEntriesFast();
1078         else            nclusters = event->GetNumberOfCaloClusters();
1079
1080   Int_t   idMax = 0;
1081   Float_t  eMax = 0;
1082   Float_t  tMax = 0;
1083   for(Int_t iclus = 0; iclus < nclusters ; iclus++)
1084   {
1085                 AliVCluster * clus = 0;
1086                 if(clusterList) clus = (AliVCluster*) (clusterList->At(iclus));
1087     else            clus = GetReader()->GetInputEvent()->GetCaloCluster(iclus);
1088
1089                 if(!clus)            continue;
1090
1091                 if(!clus->IsEMCAL()) continue;
1092
1093                 Float_t tof = clus->GetTOF()*1e9;
1094                 if(clus->E() > eMax && TMath::Abs(tof) < 30)
1095     {
1096       eMax  = clus->E();
1097                         tMax  = tof;
1098       idMax = iclus;
1099     }
1100
1101                 clus->GetMomentum(mom,GetVertex(0));
1102                 Float_t pt = mom.Pt();
1103
1104                 fhPtNPileUpSPDVtx->Fill(pt,nVtxSPD);
1105                 fhPtNPileUpTrkVtx->Fill(pt,nVtxTrk);
1106
1107                 if(TMath::Abs(tof) < 30)
1108                 {
1109                         fhPtNPileUpSPDVtxTimeCut->Fill(pt,nVtxSPD);
1110                         fhPtNPileUpTrkVtxTimeCut->Fill(pt,nVtxTrk);
1111                 }
1112
1113     if(tof < 75 && tof > -30)
1114     {
1115       fhPtNPileUpSPDVtxTimeCut2->Fill(pt,nVtxSPD);
1116       fhPtNPileUpTrkVtxTimeCut2->Fill(pt,nVtxTrk);
1117     }
1118
1119     fhTimePtNoCut->Fill(pt,tof);
1120     if(GetReader()->IsPileUpFromSPD()) fhTimePtSPD->Fill(pt,tof);
1121
1122   }
1123
1124   if(eMax < 5) return;
1125
1126   // Loop again on clusters to compare this max cluster t and the rest of the clusters, if E > 0.3
1127   Int_t n20  = 0;
1128   Int_t n40  = 0;
1129   Int_t n    = 0;
1130   Int_t nOK  = 0;
1131
1132   for(Int_t iclus = 0; iclus < nclusters ; iclus++)
1133   {
1134                 AliVCluster * clus = 0;
1135                 if(clusterList) clus = (AliVCluster*) (clusterList->At(iclus));
1136     else            clus = GetReader()->GetInputEvent()->GetCaloCluster(iclus);
1137
1138                 if(!clus)            continue;
1139
1140                 if(!clus->IsEMCAL()) continue;
1141
1142     if(clus->E() < 0.3 || iclus==idMax) continue;
1143
1144     Float_t tdiff = TMath::Abs(tMax-clus->GetTOF()*1e9);
1145     n++;
1146     if(tdiff < 25) nOK++;
1147     else
1148     {
1149       n20++;
1150       if(tdiff > 40 ) n40++;
1151     }
1152   }
1153
1154   // Check pile-up and fill histograms depending on the different cluster multiplicities
1155   if(GetReader()->IsPileUpFromSPD())
1156   {
1157     fhClusterMultSPDPileUp[0]->Fill(eMax,n  );
1158     fhClusterMultSPDPileUp[1]->Fill(eMax,nOK);
1159     fhClusterMultSPDPileUp[2]->Fill(eMax,n20);
1160     fhClusterMultSPDPileUp[3]->Fill(eMax,n40);
1161   }
1162   else
1163   {
1164     fhClusterMultNoPileUp[0]->Fill(eMax,n  );
1165     fhClusterMultNoPileUp[1]->Fill(eMax,nOK);
1166     fhClusterMultNoPileUp[2]->Fill(eMax,n20);
1167     fhClusterMultNoPileUp[3]->Fill(eMax,n40);
1168   }
1169
1170 }
1171
1172
1173 //_________________________________________________________________________________________________
1174 void AliAnaPhoton::FillPileUpHistograms(Float_t energy, Float_t pt, Float_t time)
1175 {
1176   // Fill some histograms to understand pile-up
1177   if(!fFillPileUpHistograms) return;
1178
1179   //printf("E %f, time %f\n",energy,time);
1180   AliVEvent * event = GetReader()->GetInputEvent();
1181
1182   if(GetReader()->IsPileUpFromSPD())               fhPtPhotonPileUp[0]->Fill(pt);
1183   if(GetReader()->IsPileUpFromEMCal())             fhPtPhotonPileUp[1]->Fill(pt);
1184   if(GetReader()->IsPileUpFromSPDOrEMCal())        fhPtPhotonPileUp[2]->Fill(pt);
1185   if(GetReader()->IsPileUpFromSPDAndEMCal())       fhPtPhotonPileUp[3]->Fill(pt);
1186   if(GetReader()->IsPileUpFromSPDAndNotEMCal())    fhPtPhotonPileUp[4]->Fill(pt);
1187   if(GetReader()->IsPileUpFromEMCalAndNotSPD())    fhPtPhotonPileUp[5]->Fill(pt);
1188   if(GetReader()->IsPileUpFromNotSPDAndNotEMCal()) fhPtPhotonPileUp[6]->Fill(pt);
1189
1190   fhTimePtPhotonNoCut->Fill(pt,time);
1191   if(GetReader()->IsPileUpFromSPD()) fhTimePtPhotonSPD->Fill(pt,time);
1192
1193   if(energy < 8) return; // Fill time figures for high energy clusters not too close to trigger threshold
1194
1195   AliESDEvent* esdEv = dynamic_cast<AliESDEvent*> (event);
1196   AliAODEvent* aodEv = dynamic_cast<AliAODEvent*> (event);
1197
1198   // N pile up vertices
1199   Int_t nVtxSPD = -1;
1200   Int_t nVtxTrk = -1;
1201
1202   if      (esdEv)
1203   {
1204     nVtxSPD = esdEv->GetNumberOfPileupVerticesSPD();
1205     nVtxTrk = esdEv->GetNumberOfPileupVerticesTracks();
1206
1207   }//ESD
1208   else if (aodEv)
1209   {
1210     nVtxSPD = aodEv->GetNumberOfPileupVerticesSPD();
1211     nVtxTrk = aodEv->GetNumberOfPileupVerticesTracks();
1212   }//AOD
1213
1214   fhTimeNPileUpVertSPD  ->Fill(time,nVtxSPD);
1215   fhTimeNPileUpVertTrack->Fill(time,nVtxTrk);
1216
1217   fhPtPhotonNPileUpSPDVtx->Fill(pt,nVtxSPD);
1218   fhPtPhotonNPileUpTrkVtx->Fill(pt,nVtxTrk);
1219
1220   if(TMath::Abs(time) < 25)
1221   {
1222     fhPtPhotonNPileUpSPDVtxTimeCut->Fill(pt,nVtxSPD);
1223     fhPtPhotonNPileUpTrkVtxTimeCut->Fill(pt,nVtxTrk);
1224   }
1225
1226   if(time < 75 && time > -25)
1227   {
1228     fhPtPhotonNPileUpSPDVtxTimeCut2->Fill(pt,nVtxSPD);
1229     fhPtPhotonNPileUpTrkVtxTimeCut2->Fill(pt,nVtxTrk);
1230   }
1231
1232   //printf("Is SPD %d, Is SPD Multi %d, n spd %d, n track %d\n",
1233   //       GetReader()->IsPileUpFromSPD(),event->IsPileupFromSPDInMultBins(),nVtxSPD,nVtxTrk);
1234
1235   Int_t ncont = -1;
1236   Float_t z1 = -1, z2 = -1;
1237   Float_t diamZ = -1;
1238   for(Int_t iVert=0; iVert<nVtxSPD;iVert++)
1239   {
1240     if      (esdEv)
1241     {
1242       const AliESDVertex* pv=esdEv->GetPileupVertexSPD(iVert);
1243       ncont=pv->GetNContributors();
1244       z1 = esdEv->GetPrimaryVertexSPD()->GetZ();
1245       z2 = pv->GetZ();
1246       diamZ = esdEv->GetDiamondZ();
1247     }//ESD
1248     else if (aodEv)
1249     {
1250       AliAODVertex *pv=aodEv->GetVertex(iVert);
1251       if(pv->GetType()!=AliAODVertex::kPileupSPD) continue;
1252       ncont=pv->GetNContributors();
1253       z1=aodEv->GetPrimaryVertexSPD()->GetZ();
1254       z2=pv->GetZ();
1255       diamZ = aodEv->GetDiamondZ();
1256     }// AOD
1257
1258     Double_t distZ  = TMath::Abs(z2-z1);
1259     diamZ  = TMath::Abs(z2-diamZ);
1260
1261     fhTimeNPileUpVertContributors  ->Fill(time,ncont);
1262     fhTimePileUpMainVertexZDistance->Fill(time,distZ);
1263     fhTimePileUpMainVertexZDiamond ->Fill(time,diamZ);
1264
1265   }// loop
1266 }
1267
1268 //____________________________________________________________________________________
1269 void  AliAnaPhoton::FillShowerShapeHistograms(AliVCluster* cluster, Int_t mcTag)
1270 {
1271   //Fill cluster Shower Shape histograms
1272
1273   if(!fFillSSHistograms || GetMixedEvent()) return;
1274
1275   Float_t energy  = cluster->E();
1276   Int_t   ncells  = cluster->GetNCells();
1277   Float_t lambda0 = cluster->GetM02();
1278   Float_t lambda1 = cluster->GetM20();
1279   Float_t disp    = cluster->GetDispersion()*cluster->GetDispersion();
1280
1281   TLorentzVector mom;
1282   if(GetReader()->GetDataType() != AliCaloTrackReader::kMC)
1283   {
1284     cluster->GetMomentum(mom,GetVertex(0)) ;
1285   }//Assume that come from vertex in straight line
1286   else
1287   {
1288     Double_t vertex[]={0,0,0};
1289     cluster->GetMomentum(mom,vertex) ;
1290   }
1291
1292   Float_t eta = mom.Eta();
1293   Float_t phi = mom.Phi();
1294   if(phi < 0) phi+=TMath::TwoPi();
1295
1296   fhLam0E ->Fill(energy,lambda0);
1297   fhLam1E ->Fill(energy,lambda1);
1298   fhDispE ->Fill(energy,disp);
1299
1300   if(fCalorimeter == "EMCAL" && GetModuleNumber(cluster) > 5)
1301   {
1302     fhLam0ETRD->Fill(energy,lambda0);
1303     fhLam1ETRD->Fill(energy,lambda1);
1304     fhDispETRD->Fill(energy,disp);
1305   }
1306
1307   Float_t l0   = 0., l1   = 0.;
1308   Float_t dispp= 0., dEta = 0., dPhi    = 0.;
1309   Float_t sEta = 0., sPhi = 0., sEtaPhi = 0.;
1310   if(fCalorimeter == "EMCAL" && !fFillOnlySimpleSSHisto)
1311   {
1312     GetCaloUtils()->GetEMCALRecoUtils()->RecalculateClusterShowerShapeParameters(GetEMCALGeometry(), GetReader()->GetInputEvent()->GetEMCALCells(), cluster,
1313                                                                                  l0, l1, dispp, dEta, dPhi, sEta, sPhi, sEtaPhi);
1314     //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",
1315     //       l0, l1, dispp, dEta, dPhi, sEta, sPhi, sEtaPhi );
1316     //printf("AliAnaPhoton::FillShowerShapeHistogram - dispersion %f, dispersion eta+phi %f \n",
1317     //       disp, dPhi+dEta );
1318     fhDispEtaE        -> Fill(energy,dEta);
1319     fhDispPhiE        -> Fill(energy,dPhi);
1320     fhSumEtaE         -> Fill(energy,sEta);
1321     fhSumPhiE         -> Fill(energy,sPhi);
1322     fhSumEtaPhiE      -> Fill(energy,sEtaPhi);
1323     fhDispEtaPhiDiffE -> Fill(energy,dPhi-dEta);
1324     if(dEta+dPhi>0)fhSphericityE     -> Fill(energy,(dPhi-dEta)/(dEta+dPhi));
1325     if(dEta+sEta>0)fhDispSumEtaDiffE -> Fill(energy,(dEta-sEta)/((dEta+sEta)/2.));
1326     if(dPhi+sPhi>0)fhDispSumPhiDiffE -> Fill(energy,(dPhi-sPhi)/((dPhi+sPhi)/2.));
1327
1328     Int_t ebin = -1;
1329     if      (energy < 2 ) ebin = 0;
1330     else if (energy < 4 ) ebin = 1;
1331     else if (energy < 6 ) ebin = 2;
1332     else if (energy < 10) ebin = 3;
1333     else if (energy < 15) ebin = 4;
1334     else if (energy < 20) ebin = 5;
1335     else                  ebin = 6;
1336
1337     fhDispEtaDispPhi[ebin]->Fill(dEta   ,dPhi);
1338     fhLambda0DispEta[ebin]->Fill(lambda0,dEta);
1339     fhLambda0DispPhi[ebin]->Fill(lambda0,dPhi);
1340
1341   }
1342
1343   // if track-matching was of, check effect of track-matching residual cut
1344
1345   if(!fRejectTrackMatch)
1346   {
1347     Float_t dZ  = cluster->GetTrackDz();
1348     Float_t dR  = cluster->GetTrackDx();
1349     if(cluster->IsEMCAL() && GetCaloUtils()->IsRecalculationOfClusterTrackMatchingOn())
1350     {
1351       dR = 2000., dZ = 2000.;
1352       GetCaloUtils()->GetEMCALRecoUtils()->GetMatchedResiduals(cluster->GetID(),dZ,dR);
1353     }
1354
1355     if(TMath::Abs(dZ) < 0.05 && TMath::Abs(dR) < 0.05)
1356     {
1357       fhLam0ETM ->Fill(energy,lambda0);
1358       fhLam1ETM ->Fill(energy,lambda1);
1359       fhDispETM ->Fill(energy,disp);
1360
1361       if(fCalorimeter == "EMCAL" && GetModuleNumber(cluster) > 5)
1362       {
1363         fhLam0ETMTRD->Fill(energy,lambda0);
1364         fhLam1ETMTRD->Fill(energy,lambda1);
1365         fhDispETMTRD->Fill(energy,disp);
1366       }
1367     }
1368   }// if track-matching was of, check effect of matching residual cut
1369
1370
1371   if(!fFillOnlySimpleSSHisto)
1372   {
1373     if(energy < 2)
1374     {
1375       fhNCellsLam0LowE ->Fill(ncells,lambda0);
1376       fhNCellsLam1LowE ->Fill(ncells,lambda1);
1377       fhNCellsDispLowE ->Fill(ncells,disp);
1378
1379       fhLam1Lam0LowE  ->Fill(lambda1,lambda0);
1380       fhLam0DispLowE  ->Fill(lambda0,disp);
1381       fhDispLam1LowE  ->Fill(disp,lambda1);
1382       fhEtaLam0LowE   ->Fill(eta,lambda0);
1383       fhPhiLam0LowE   ->Fill(phi,lambda0);
1384     }
1385     else
1386     {
1387       fhNCellsLam0HighE ->Fill(ncells,lambda0);
1388       fhNCellsLam1HighE ->Fill(ncells,lambda1);
1389       fhNCellsDispHighE ->Fill(ncells,disp);
1390
1391       fhLam1Lam0HighE  ->Fill(lambda1,lambda0);
1392       fhLam0DispHighE  ->Fill(lambda0,disp);
1393       fhDispLam1HighE  ->Fill(disp,lambda1);
1394       fhEtaLam0HighE   ->Fill(eta, lambda0);
1395       fhPhiLam0HighE   ->Fill(phi, lambda0);
1396     }
1397   }
1398
1399   if(IsDataMC())
1400   {
1401     AliVCaloCells* cells = 0;
1402     if(fCalorimeter == "EMCAL") cells = GetEMCALCells();
1403     else                        cells = GetPHOSCells();
1404
1405     //Fill histograms to check shape of embedded clusters
1406     Float_t fraction = 0;
1407     // printf("check embedding %i\n",GetReader()->IsEmbeddedClusterSelectionOn());
1408
1409     if(GetReader()->IsEmbeddedClusterSelectionOn())
1410     {//Only working for EMCAL
1411       //        printf("embedded\n");
1412       Float_t clusterE = 0; // recalculate in case corrections applied.
1413       Float_t cellE    = 0;
1414       for(Int_t icell  = 0; icell < cluster->GetNCells(); icell++)
1415       {
1416         cellE    = cells->GetCellAmplitude(cluster->GetCellAbsId(icell));
1417         clusterE+=cellE;
1418         fraction+=cellE*cluster->GetCellAmplitudeFraction(icell);
1419       }
1420
1421       //Fraction of total energy due to the embedded signal
1422       fraction/=clusterE;
1423
1424       if(GetDebug() > 1 )
1425         printf("AliAnaPhoton::FillShowerShapeHistogram() - Energy fraction of embedded signal %2.3f, Energy %2.3f\n",fraction, clusterE);
1426
1427       fhEmbeddedSignalFractionEnergy->Fill(clusterE,fraction);
1428
1429     }  // embedded fraction
1430
1431     // Get the fraction of the cluster energy that carries the cell with highest energy
1432     Float_t maxCellFraction = 0.;
1433     Int_t absID = GetCaloUtils()->GetMaxEnergyCell(cells, cluster,maxCellFraction);
1434
1435     if( absID < 0 ) AliFatal("Wrong absID");
1436
1437     // Check the origin and fill histograms
1438
1439     Int_t mcIndex = -1;
1440
1441     if( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton)     &&
1442        !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCConversion) &&
1443        !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0)        &&
1444        !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEta))
1445     {
1446       mcIndex = kmcssPhoton ;
1447
1448       if(!GetReader()->IsEmbeddedClusterSelectionOn())
1449       {
1450         //Check particle overlaps in cluster
1451
1452         // Compare the primary depositing more energy with the rest,
1453         // if no photon/electron as comon ancestor (conversions), count as other particle
1454         const UInt_t nlabels = cluster->GetNLabels();
1455         Int_t overpdg[nlabels];
1456         Int_t noverlaps = GetMCAnalysisUtils()->GetNOverlaps(cluster->GetLabels(), nlabels,mcTag,-1,GetReader(),overpdg);
1457
1458         //printf("N overlaps %d \n",noverlaps);
1459
1460         if(noverlaps == 0)
1461         {
1462           fhMCPhotonELambda0NoOverlap  ->Fill(energy, lambda0);
1463         }
1464         else if(noverlaps == 1)
1465         {
1466           fhMCPhotonELambda0TwoOverlap ->Fill(energy, lambda0);
1467         }
1468         else if(noverlaps > 1)
1469         {
1470           fhMCPhotonELambda0NOverlap   ->Fill(energy, lambda0);
1471         }
1472         else
1473         {
1474           printf("AliAnaPhoton::FillShowerShapeHistogram() - n overlaps = %d!!\n", noverlaps);
1475         }
1476       }//No embedding
1477
1478       //Fill histograms to check shape of embedded clusters
1479       if(GetReader()->IsEmbeddedClusterSelectionOn())
1480       {
1481         if     (fraction > 0.9)
1482         {
1483           fhEmbedPhotonELambda0FullSignal   ->Fill(energy, lambda0);
1484         }
1485         else if(fraction > 0.5)
1486         {
1487           fhEmbedPhotonELambda0MostlySignal ->Fill(energy, lambda0);
1488         }
1489         else if(fraction > 0.1)
1490         {
1491           fhEmbedPhotonELambda0MostlyBkg    ->Fill(energy, lambda0);
1492         }
1493         else
1494         {
1495           fhEmbedPhotonELambda0FullBkg      ->Fill(energy, lambda0);
1496         }
1497       } // embedded
1498
1499     }//photon   no conversion
1500     else if  ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton)     &&
1501                GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCConversion) &&
1502               !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0)        &&
1503               !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEta))
1504     {
1505       mcIndex = kmcssConversion ;
1506     }//conversion photon
1507
1508     else if  ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCElectron))
1509     {
1510       mcIndex = kmcssElectron ;
1511     }//electron
1512     else if  ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0)  )
1513     {
1514       mcIndex = kmcssPi0 ;
1515
1516       //Fill histograms to check shape of embedded clusters
1517       if(GetReader()->IsEmbeddedClusterSelectionOn())
1518       {
1519         if     (fraction > 0.9)
1520         {
1521           fhEmbedPi0ELambda0FullSignal   ->Fill(energy, lambda0);
1522         }
1523         else if(fraction > 0.5)
1524         {
1525           fhEmbedPi0ELambda0MostlySignal ->Fill(energy, lambda0);
1526         }
1527         else if(fraction > 0.1)
1528         {
1529           fhEmbedPi0ELambda0MostlyBkg    ->Fill(energy, lambda0);
1530         }
1531         else
1532         {
1533           fhEmbedPi0ELambda0FullBkg      ->Fill(energy, lambda0);
1534         }
1535       } // embedded
1536
1537     }//pi0
1538     else if  ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEta)  )
1539     {
1540       mcIndex = kmcssEta ;
1541     }//eta
1542     else
1543     {
1544       mcIndex = kmcssOther ;
1545     }//other particles
1546
1547     fhMCELambda0           [mcIndex]->Fill(energy, lambda0);
1548     fhMCELambda1           [mcIndex]->Fill(energy, lambda1);
1549     fhMCEDispersion        [mcIndex]->Fill(energy, disp);
1550     fhMCNCellsE            [mcIndex]->Fill(energy, ncells);
1551     fhMCMaxCellDiffClusterE[mcIndex]->Fill(energy, maxCellFraction);
1552
1553     if(!fFillOnlySimpleSSHisto)
1554     {
1555       if     (energy < 2.)
1556       {
1557         fhMCLambda0vsClusterMaxCellDiffE0[mcIndex]->Fill(lambda0, maxCellFraction);
1558         fhMCNCellsvsClusterMaxCellDiffE0 [mcIndex]->Fill(ncells,  maxCellFraction);
1559       }
1560       else if(energy < 6.)
1561       {
1562         fhMCLambda0vsClusterMaxCellDiffE2[mcIndex]->Fill(lambda0, maxCellFraction);
1563         fhMCNCellsvsClusterMaxCellDiffE2 [mcIndex]->Fill(ncells,  maxCellFraction);
1564       }
1565       else
1566       {
1567         fhMCLambda0vsClusterMaxCellDiffE6[mcIndex]->Fill(lambda0, maxCellFraction);
1568         fhMCNCellsvsClusterMaxCellDiffE6 [mcIndex]->Fill(ncells,  maxCellFraction);
1569       }
1570
1571       if(fCalorimeter == "EMCAL")
1572       {
1573         fhMCEDispEta        [mcIndex]-> Fill(energy,dEta);
1574         fhMCEDispPhi        [mcIndex]-> Fill(energy,dPhi);
1575         fhMCESumEtaPhi      [mcIndex]-> Fill(energy,sEtaPhi);
1576         fhMCEDispEtaPhiDiff [mcIndex]-> Fill(energy,dPhi-dEta);
1577         if(dEta+dPhi>0)fhMCESphericity[mcIndex]-> Fill(energy,(dPhi-dEta)/(dEta+dPhi));
1578
1579         Int_t ebin = -1;
1580         if      (energy < 2 ) ebin = 0;
1581         else if (energy < 4 ) ebin = 1;
1582         else if (energy < 6 ) ebin = 2;
1583         else if (energy < 10) ebin = 3;
1584         else if (energy < 15) ebin = 4;
1585         else if (energy < 20) ebin = 5;
1586         else                  ebin = 6;
1587
1588         fhMCDispEtaDispPhi[ebin][mcIndex]->Fill(dEta   ,dPhi);
1589         fhMCLambda0DispEta[ebin][mcIndex]->Fill(lambda0,dEta);
1590         fhMCLambda0DispPhi[ebin][mcIndex]->Fill(lambda0,dPhi);
1591       }
1592     }
1593   }//MC data
1594
1595 }
1596
1597 //__________________________________________________________________________
1598 void AliAnaPhoton::FillTrackMatchingResidualHistograms(AliVCluster* cluster,
1599                                                        Int_t cut)
1600 {
1601   // If selected, fill histograms with residuals of matched clusters, help to define track matching cut
1602   // Residual filled for different cuts 0 (No cut), after 1 PID cut
1603
1604   Float_t dZ = cluster->GetTrackDz();
1605   Float_t dR = cluster->GetTrackDx();
1606
1607   if(cluster->IsEMCAL() && GetCaloUtils()->IsRecalculationOfClusterTrackMatchingOn())
1608   {
1609     dR = 2000., dZ = 2000.;
1610     GetCaloUtils()->GetEMCALRecoUtils()->GetMatchedResiduals(cluster->GetID(),dZ,dR);
1611   }
1612
1613   AliVTrack *track = GetCaloUtils()->GetMatchedTrack(cluster, GetReader()->GetInputEvent());
1614
1615   Bool_t positive = kFALSE;
1616   if(track) positive = (track->Charge()>0);
1617
1618   if(fhTrackMatchedDEta[cut] && TMath::Abs(dR) < 999)
1619   {
1620     fhTrackMatchedDEta[cut]->Fill(cluster->E(),dZ);
1621     fhTrackMatchedDPhi[cut]->Fill(cluster->E(),dR);
1622     if(cluster->E() > 0.5) fhTrackMatchedDEtaDPhi[cut]->Fill(dZ,dR);
1623
1624     if(track)
1625     {
1626       if(positive)
1627       {
1628         fhTrackMatchedDEtaPos[cut]->Fill(cluster->E(),dZ);
1629         fhTrackMatchedDPhiPos[cut]->Fill(cluster->E(),dR);
1630         if(cluster->E() > 0.5) fhTrackMatchedDEtaDPhiPos[cut]->Fill(dZ,dR);
1631       }
1632       else
1633       {
1634         fhTrackMatchedDEtaNeg[cut]->Fill(cluster->E(),dZ);
1635         fhTrackMatchedDPhiNeg[cut]->Fill(cluster->E(),dR);
1636         if(cluster->E() > 0.5) fhTrackMatchedDEtaDPhiNeg[cut]->Fill(dZ,dR);
1637       }
1638     }
1639
1640     Int_t nSMod = GetModuleNumber(cluster);
1641
1642     if(fCalorimeter=="EMCAL" &&  nSMod > 5)
1643     {
1644       fhTrackMatchedDEtaTRD[cut]->Fill(cluster->E(),dZ);
1645       fhTrackMatchedDPhiTRD[cut]->Fill(cluster->E(),dR);
1646     }
1647
1648     // Check dEdx and E/p of matched clusters
1649
1650     if(TMath::Abs(dZ) < 0.05 && TMath::Abs(dR) < 0.05)
1651     {
1652       if(track)
1653       {
1654         Float_t dEdx   = track->GetTPCsignal();
1655         Float_t eOverp = cluster->E()/track->P();
1656
1657         fhdEdx[cut]  ->Fill(cluster->E(), dEdx);
1658         fhEOverP[cut]->Fill(cluster->E(), eOverp);
1659
1660         if(fCalorimeter=="EMCAL" && nSMod > 5)
1661           fhEOverPTRD[cut]->Fill(cluster->E(), eOverp);
1662
1663
1664       }
1665       else
1666         printf("AliAnaPhoton::FillTrackMatchingResidualHistograms() - Residual OK but (dR, dZ)= (%2.4f,%2.4f) no track associated WHAT? \n", dR,dZ);
1667
1668
1669
1670       if(IsDataMC())
1671       {
1672
1673         Int_t tag = GetMCAnalysisUtils()->CheckOrigin(cluster->GetLabels(),cluster->GetNLabels(),GetReader());
1674
1675         if  ( !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion)  )
1676         {
1677           if       ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0)      ||
1678                     GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta)       ) fhTrackMatchedMCParticle[cut]->Fill(cluster->E(), 2.5 );
1679           else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton)    ) fhTrackMatchedMCParticle[cut]->Fill(cluster->E(), 0.5 );
1680           else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron)  ) fhTrackMatchedMCParticle[cut]->Fill(cluster->E(), 1.5 );
1681           else                                                                                 fhTrackMatchedMCParticle[cut]->Fill(cluster->E(), 3.5 );
1682
1683           // Check if several particles contributed to cluster and discard overlapped mesons
1684           if(!GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) ||
1685              !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta))
1686           {
1687             if(cluster->GetNLabels()==1)
1688             {
1689               fhTrackMatchedDEtaMCNoOverlap[cut]->Fill(cluster->E(),dZ);
1690               fhTrackMatchedDPhiMCNoOverlap[cut]->Fill(cluster->E(),dR);
1691             }
1692             else
1693             {
1694               fhTrackMatchedDEtaMCOverlap[cut]->Fill(cluster->E(),dZ);
1695               fhTrackMatchedDPhiMCOverlap[cut]->Fill(cluster->E(),dR);
1696             }
1697
1698           }// Check overlaps
1699
1700         }
1701         else
1702         {
1703           if       ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0)      ||
1704                     GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta)       ) fhTrackMatchedMCParticle[cut]->Fill(cluster->E(), 6.5 );
1705           else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton)    ) fhTrackMatchedMCParticle[cut]->Fill(cluster->E(), 4.5 );
1706           else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron)  ) fhTrackMatchedMCParticle[cut]->Fill(cluster->E(), 5.5 );
1707           else                                                                                 fhTrackMatchedMCParticle[cut]->Fill(cluster->E(), 7.5 );
1708
1709           // Check if several particles contributed to cluster
1710           if(!GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) ||
1711              !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta))
1712           {
1713             fhTrackMatchedDEtaMCConversion[cut]->Fill(cluster->E(),dZ);
1714             fhTrackMatchedDPhiMCConversion[cut]->Fill(cluster->E(),dR);
1715
1716           }// Check overlaps
1717
1718         }
1719
1720       } // MC
1721
1722     } // residuals window
1723
1724   } // Small residual
1725
1726 }
1727
1728 //___________________________________________
1729 TObjString *  AliAnaPhoton::GetAnalysisCuts()
1730 {
1731   //Save parameters used for analysis
1732   TString parList ; //this will be list of parameters used for this analysis.
1733   const Int_t buffersize = 255;
1734   char onePar[buffersize] ;
1735
1736   snprintf(onePar,buffersize,"--- AliAnaPhoton ---\n") ;
1737   parList+=onePar ;
1738   snprintf(onePar,buffersize,"Calorimeter: %s\n",fCalorimeter.Data()) ;
1739   parList+=onePar ;
1740   snprintf(onePar,buffersize,"fMinDist =%2.2f (Minimal distance to bad channel to accept cluster) \n",fMinDist) ;
1741   parList+=onePar ;
1742   snprintf(onePar,buffersize,"fMinDist2=%2.2f (Cuts on Minimal distance to study acceptance evaluation) \n",fMinDist2) ;
1743   parList+=onePar ;
1744   snprintf(onePar,buffersize,"fMinDist3=%2.2f (One more cut on distance used for acceptance-efficiency study) \n",fMinDist3) ;
1745   parList+=onePar ;
1746   snprintf(onePar,buffersize,"fRejectTrackMatch: %d\n",fRejectTrackMatch) ;
1747   parList+=onePar ;
1748
1749   //Get parameters set in base class.
1750   parList += GetBaseParametersList() ;
1751
1752   //Get parameters set in PID class.
1753   parList += GetCaloPID()->GetPIDParametersList() ;
1754
1755   //Get parameters set in FiducialCut class (not available yet)
1756   //parlist += GetFidCut()->GetFidCutParametersList()
1757
1758   return new TObjString(parList) ;
1759 }
1760
1761 //________________________________________________________________________
1762 TList *  AliAnaPhoton::GetCreateOutputObjects()
1763 {
1764   // Create histograms to be saved in output file and
1765   // store them in outputContainer
1766   TList * outputContainer = new TList() ;
1767   outputContainer->SetName("PhotonHistos") ;
1768
1769   Int_t nptbins  = GetHistogramRanges()->GetHistoPtBins();  Float_t ptmax  = GetHistogramRanges()->GetHistoPtMax();  Float_t ptmin  = GetHistogramRanges()->GetHistoPtMin();
1770   Int_t nphibins = GetHistogramRanges()->GetHistoPhiBins(); Float_t phimax = GetHistogramRanges()->GetHistoPhiMax(); Float_t phimin = GetHistogramRanges()->GetHistoPhiMin();
1771   Int_t netabins = GetHistogramRanges()->GetHistoEtaBins(); Float_t etamax = GetHistogramRanges()->GetHistoEtaMax(); Float_t etamin = GetHistogramRanges()->GetHistoEtaMin();
1772   Int_t ssbins   = GetHistogramRanges()->GetHistoShowerShapeBins();  Float_t ssmax   = GetHistogramRanges()->GetHistoShowerShapeMax();  Float_t ssmin   = GetHistogramRanges()->GetHistoShowerShapeMin();
1773   Int_t nbins    = GetHistogramRanges()->GetHistoNClusterCellBins(); Int_t   nmax    = GetHistogramRanges()->GetHistoNClusterCellMax(); Int_t   nmin    = GetHistogramRanges()->GetHistoNClusterCellMin();
1774   Int_t ntimebins= GetHistogramRanges()->GetHistoTimeBins();         Float_t timemax = GetHistogramRanges()->GetHistoTimeMax();         Float_t timemin = GetHistogramRanges()->GetHistoTimeMin();
1775
1776   Int_t   nresetabins = GetHistogramRanges()->GetHistoTrackResidualEtaBins();
1777   Float_t resetamax   = GetHistogramRanges()->GetHistoTrackResidualEtaMax();
1778   Float_t resetamin   = GetHistogramRanges()->GetHistoTrackResidualEtaMin();
1779   Int_t   nresphibins = GetHistogramRanges()->GetHistoTrackResidualPhiBins();
1780   Float_t resphimax   = GetHistogramRanges()->GetHistoTrackResidualPhiMax();
1781   Float_t resphimin   = GetHistogramRanges()->GetHistoTrackResidualPhiMin();
1782
1783   Int_t   ndedxbins   = GetHistogramRanges()->GetHistodEdxBins();
1784   Float_t dedxmax     = GetHistogramRanges()->GetHistodEdxMax();
1785   Float_t dedxmin     = GetHistogramRanges()->GetHistodEdxMin();
1786   Int_t   nPoverEbins = GetHistogramRanges()->GetHistoPOverEBins();
1787   Float_t pOverEmax   = GetHistogramRanges()->GetHistoPOverEMax();
1788   Float_t pOverEmin   = GetHistogramRanges()->GetHistoPOverEMin();
1789
1790   Int_t bin[] = {0,2,4,6,10,15,20,100}; // energy bins for SS studies
1791
1792   Int_t nTrigBC  = 1;
1793   Int_t iBCShift = 0;
1794   if(!GetReader()->AreBadTriggerEventsRemoved())
1795   {
1796     nTrigBC = 11;
1797     iBCShift = 5;
1798   }
1799
1800   TString cut[] = {"Open","Reader","E","Time","NCells","NLM","Fidutial","Matching","Bad","PID"};
1801   for (Int_t i = 0; i < 10 ;  i++)
1802   {
1803     fhClusterCuts[i] = new TH1F(Form("hCut_%d_%s", i, cut[i].Data()),
1804                                 Form("Number of clusters that pass cuts <= %d, %s", i, cut[i].Data()),
1805                                 nptbins,ptmin,ptmax);
1806     fhClusterCuts[i]->SetYTitle("d#it{N}/d#it{E} ");
1807     fhClusterCuts[i]->SetXTitle("#it{E} (GeV)");
1808     outputContainer->Add(fhClusterCuts[i]) ;
1809   }
1810
1811   fhEClusterSM = new TH2F("hEClusterSM","Raw clusters E and super-module number",
1812                           nptbins,ptmin,ptmax,
1813                           GetCaloUtils()->GetNumberOfSuperModulesUsed(),0,GetCaloUtils()->GetNumberOfSuperModulesUsed());
1814   fhEClusterSM->SetYTitle("SuperModule ");
1815   fhEClusterSM->SetXTitle("#it{E} (GeV)");
1816   outputContainer->Add(fhEClusterSM) ;
1817
1818   fhPtClusterSM = new TH2F("hPtClusterSM","Raw clusters #it{p}_[T] and super-module number",
1819                           nptbins,ptmin,ptmax,
1820                           GetCaloUtils()->GetNumberOfSuperModulesUsed(),0,GetCaloUtils()->GetNumberOfSuperModulesUsed());
1821   fhPtClusterSM->SetYTitle("SuperModule ");
1822   fhPtClusterSM->SetXTitle("#it{E} (GeV)");
1823   outputContainer->Add(fhPtClusterSM) ;
1824
1825   fhEPhotonSM = new TH2F("hEPhotonSM","Selected clusters E and super-module number",
1826                           nptbins,ptmin,ptmax,
1827                           GetCaloUtils()->GetNumberOfSuperModulesUsed(),0,GetCaloUtils()->GetNumberOfSuperModulesUsed());
1828   fhEPhotonSM->SetYTitle("SuperModule ");
1829   fhEPhotonSM->SetXTitle("#it{E} (GeV)");
1830   outputContainer->Add(fhEPhotonSM) ;
1831
1832   fhPtPhotonSM = new TH2F("hPtPhotonSM","Selected clusters #it{p}_[T] and super-module number",
1833                            nptbins,ptmin,ptmax,
1834                            GetCaloUtils()->GetNumberOfSuperModulesUsed(),0,GetCaloUtils()->GetNumberOfSuperModulesUsed());
1835   fhPtPhotonSM->SetYTitle("SuperModule ");
1836   fhPtPhotonSM->SetXTitle("#it{E} (GeV)");
1837   outputContainer->Add(fhPtPhotonSM) ;
1838
1839   fhNCellsE  = new TH2F ("hNCellsE","# of cells in cluster vs E of clusters", nptbins,ptmin,ptmax, nbins,nmin,nmax);
1840   fhNCellsE->SetXTitle("#it{E} (GeV)");
1841   fhNCellsE->SetYTitle("# of cells in cluster");
1842   outputContainer->Add(fhNCellsE);
1843
1844   fhCellsE  = new TH2F ("hCellsE","energy of cells in cluster vs E of clusters", nptbins,ptmin,ptmax, nptbins*2,ptmin,ptmax);
1845   fhCellsE->SetXTitle("#it{E}_{cluster} (GeV)");
1846   fhCellsE->SetYTitle("#it{E}_{cell} (GeV)");
1847   outputContainer->Add(fhCellsE);
1848
1849   fhTimePt  = new TH2F ("hTimePt","time of cluster vs pT of clusters", nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
1850   fhTimePt->SetXTitle("#it{p}_[T] (GeV/#it{c})");
1851   fhTimePt->SetYTitle("#it{time} (ns)");
1852   outputContainer->Add(fhTimePt);
1853
1854   fhMaxCellDiffClusterE  = new TH2F ("hMaxCellDiffClusterE","energy vs difference of cluster energy - max cell energy / cluster energy, good clusters",
1855                                      nptbins,ptmin,ptmax, 500,0,1.);
1856   fhMaxCellDiffClusterE->SetXTitle("#it{E}_{cluster} (GeV) ");
1857   fhMaxCellDiffClusterE->SetYTitle("(#it{E}_{cluster} - #it{E}_{cell max})/ #it{E}_{cluster}");
1858   outputContainer->Add(fhMaxCellDiffClusterE);
1859
1860   fhEPhoton  = new TH1F("hEPhoton","Number of #gamma over calorimeter vs energy",nptbins,ptmin,ptmax);
1861   fhEPhoton->SetYTitle("#it{counts}");
1862   fhEPhoton->SetXTitle("#it{E}_{#gamma}(GeV)");
1863   outputContainer->Add(fhEPhoton) ;
1864
1865   fhPtPhoton  = new TH1F("hPtPhoton","Number of #gamma over calorimeter vs #it{p}_[T]",nptbins,ptmin,ptmax);
1866   fhPtPhoton->SetYTitle("#it{counts}");
1867   fhPtPhoton->SetXTitle("p_{T #gamma}(GeV/#it{c})");
1868   outputContainer->Add(fhPtPhoton) ;
1869
1870   fhPtCentralityPhoton  = new TH2F("hPtCentralityPhoton","centrality vs #it{p}_[T]",nptbins,ptmin,ptmax, 100,0,100);
1871   fhPtCentralityPhoton->SetYTitle("Centrality");
1872   fhPtCentralityPhoton->SetXTitle("#it{p}_[T](GeV/#it{c})");
1873   outputContainer->Add(fhPtCentralityPhoton) ;
1874
1875   fhPtEventPlanePhoton  = new TH2F("hPtEventPlanePhoton","centrality vs #it{p}_[T]",nptbins,ptmin,ptmax, 100,0,TMath::Pi());
1876   fhPtEventPlanePhoton->SetYTitle("Event plane angle (rad)");
1877   fhPtEventPlanePhoton->SetXTitle("#it{p}_[T] (GeV/#it{c})");
1878   outputContainer->Add(fhPtEventPlanePhoton) ;
1879
1880   fhEtaPhi  = new TH2F
1881   ("hEtaPhi","cluster,#it{E} > 0.5 GeV, #eta vs #phi",netabins,etamin,etamax,nphibins,phimin,phimax);
1882   fhEtaPhi->SetYTitle("#phi (rad)");
1883   fhEtaPhi->SetXTitle("#eta");
1884   outputContainer->Add(fhEtaPhi) ;
1885
1886   if(fCalorimeter=="EMCAL" && fFillEMCALBCHistograms)
1887   {
1888     fhEtaPhiEMCALBC0  = new TH2F
1889     ("hEtaPhiEMCALBC0","cluster,#it{E} > 2 GeV, #eta vs #phi, for clusters with |time| < 25 ns, EMCAL-BC=0",netabins,etamin,etamax,nphibins,phimin,phimax);
1890     fhEtaPhiEMCALBC0->SetYTitle("#phi (rad)");
1891     fhEtaPhiEMCALBC0->SetXTitle("#eta");
1892     outputContainer->Add(fhEtaPhiEMCALBC0) ;
1893
1894     fhEtaPhiEMCALBC1  = new TH2F
1895     ("hEtaPhiEMCALBC1","cluster,#it{E} > 2 GeV, #eta vs #phi, for clusters with 25 < |time| < 75 ns, EMCAL-BC=1",netabins,etamin,etamax,nphibins,phimin,phimax);
1896     fhEtaPhiEMCALBC1->SetYTitle("#phi (rad)");
1897     fhEtaPhiEMCALBC1->SetXTitle("#eta");
1898     outputContainer->Add(fhEtaPhiEMCALBC1) ;
1899
1900     fhEtaPhiEMCALBCN  = new TH2F
1901     ("hEtaPhiEMCALBCN","cluster,#it{E} > 2 GeV, #eta vs #phi, for clusters with |time| > 75 ns, EMCAL-BC>1",netabins,etamin,etamax,nphibins,phimin,phimax);
1902     fhEtaPhiEMCALBCN->SetYTitle("#phi (rad)");
1903     fhEtaPhiEMCALBCN->SetXTitle("#eta");
1904     outputContainer->Add(fhEtaPhiEMCALBCN) ;
1905
1906     for(Int_t i = 0; i < nTrigBC; i++)
1907     {
1908       fhEtaPhiTriggerEMCALBC[i] = new TH2F
1909       (Form("hEtaPhiTriggerEMCALBC%d",i-iBCShift),
1910        Form("cluster #it{E} > 2 GeV, #eta vs #phi, Trigger EMCAL-BC=%d",i-iBCShift),
1911        netabins,etamin,etamax,nphibins,phimin,phimax);
1912       fhEtaPhiTriggerEMCALBC[i]->SetYTitle("#phi (rad)");
1913       fhEtaPhiTriggerEMCALBC[i]->SetXTitle("#eta");
1914       outputContainer->Add(fhEtaPhiTriggerEMCALBC[i]) ;
1915
1916       fhTimeTriggerEMCALBC[i] = new TH2F
1917       (Form("hTimeTriggerEMCALBC%d",i-iBCShift),
1918        Form("cluster #it{time} vs #it{E} of clusters, Trigger EMCAL-BC=%d",i-iBCShift),
1919        nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
1920       fhTimeTriggerEMCALBC[i]->SetXTitle("#it{E} (GeV)");
1921       fhTimeTriggerEMCALBC[i]->SetYTitle("#it{time} (ns)");
1922       outputContainer->Add(fhTimeTriggerEMCALBC[i]);
1923
1924       fhTimeTriggerEMCALBCPileUpSPD[i] = new TH2F
1925       (Form("hTimeTriggerEMCALBC%dPileUpSPD",i-iBCShift),
1926        Form("cluster #it{time} vs #it{E} of clusters, Trigger EMCAL-BC=%d",i-iBCShift),
1927        nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
1928       fhTimeTriggerEMCALBCPileUpSPD[i]->SetXTitle("#it{E} (GeV)");
1929       fhTimeTriggerEMCALBCPileUpSPD[i]->SetYTitle("#it{time} (ns)");
1930       outputContainer->Add(fhTimeTriggerEMCALBCPileUpSPD[i]);
1931
1932       fhEtaPhiTriggerEMCALBCUM[i] = new TH2F
1933       (Form("hEtaPhiTriggerEMCALBC%d_UnMatch",i-iBCShift),
1934        Form("cluster #it{E} > 2 GeV, #eta vs #phi, unmatched trigger EMCAL-BC=%d",i-iBCShift),
1935        netabins,etamin,etamax,nphibins,phimin,phimax);
1936       fhEtaPhiTriggerEMCALBCUM[i]->SetYTitle("#phi (rad)");
1937       fhEtaPhiTriggerEMCALBCUM[i]->SetXTitle("#eta");
1938       outputContainer->Add(fhEtaPhiTriggerEMCALBCUM[i]) ;
1939
1940       fhTimeTriggerEMCALBCUM[i] = new TH2F
1941       (Form("hTimeTriggerEMCALBC%d_UnMatch",i-iBCShift),
1942        Form("cluster #it{time} vs #it{E} of clusters, unmatched trigger EMCAL-BC=%d",i-iBCShift),
1943        nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
1944       fhTimeTriggerEMCALBCUM[i]->SetXTitle("#it{E} (GeV)");
1945       fhTimeTriggerEMCALBCUM[i]->SetYTitle("#it{time} (ns)");
1946       outputContainer->Add(fhTimeTriggerEMCALBCUM[i]);
1947
1948       fhEtaPhiTriggerEMCALBCCluster[i] = new TH2F
1949       (Form("hEtaPhiTriggerEMCALBC%d_OnlyTrigger",i-iBCShift),
1950        Form("trigger cluster, #eta vs #phi, Trigger EMCAL-BC=%d",i-iBCShift),
1951        netabins,etamin,etamax,nphibins,phimin,phimax);
1952       fhEtaPhiTriggerEMCALBCCluster[i]->SetYTitle("#phi (rad)");
1953       fhEtaPhiTriggerEMCALBCCluster[i]->SetXTitle("#eta");
1954       outputContainer->Add(fhEtaPhiTriggerEMCALBCCluster[i]) ;
1955
1956       fhEtaPhiTriggerEMCALBCUMCluster[i] = new TH2F
1957       (Form("hEtaPhiTriggerEMCALBC%d_OnlyTrigger_UnMatch",i-iBCShift),
1958        Form("trigger cluster, #eta vs #phi, unmatched trigger EMCAL-BC=%d",i-iBCShift),
1959        netabins,etamin,etamax,nphibins,phimin,phimax);
1960       fhEtaPhiTriggerEMCALBCUMCluster[i]->SetYTitle("#phi (rad)");
1961       fhEtaPhiTriggerEMCALBCUMCluster[i]->SetXTitle("#eta");
1962       outputContainer->Add(fhEtaPhiTriggerEMCALBCUMCluster[i]) ;
1963     }
1964
1965     fhTimeTriggerEMCALBCCluster = new TH2F("hTimeTriggerEMCALBC_OnlyTrigger",
1966                                            "trigger cluster #it{time} vs #it{E} of clusters",
1967                                            nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
1968     fhTimeTriggerEMCALBCCluster->SetXTitle("#it{E} (GeV)");
1969     fhTimeTriggerEMCALBCCluster->SetYTitle("#it{time} (ns)");
1970     outputContainer->Add(fhTimeTriggerEMCALBCCluster);
1971
1972     fhTimeTriggerEMCALBCUMCluster = new TH2F("hTimeTriggerEMCALBC_OnlyTrigger_UnMatch",
1973                                              "trigger cluster #it{time} vs #it{E} of clusters, unmatched trigger",
1974                                              nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
1975     fhTimeTriggerEMCALBCUMCluster->SetXTitle("#it{E} (GeV)");
1976     fhTimeTriggerEMCALBCUMCluster->SetYTitle("#it{time} (ns)");
1977     outputContainer->Add(fhTimeTriggerEMCALBCUMCluster);
1978
1979     fhEtaPhiTriggerEMCALBCClusterOverTh = new TH2F
1980     ("hEtaPhiTriggerEMCALBC0_OnlyTrigger_OverThreshold",
1981      "trigger cluster #it{E} > trigger threshold, #eta vs #phi, Trigger EMCAL-BC=0",
1982      netabins,etamin,etamax,nphibins,phimin,phimax);
1983     fhEtaPhiTriggerEMCALBCClusterOverTh->SetYTitle("#phi (rad)");
1984     fhEtaPhiTriggerEMCALBCClusterOverTh->SetXTitle("#eta");
1985     outputContainer->Add(fhEtaPhiTriggerEMCALBCClusterOverTh) ;
1986
1987     fhEtaPhiTriggerEMCALBCUMClusterOverTh = new TH2F
1988     ("hEtaPhiTriggerEMCALBC0_OnlyTrigger_OverThreshold_UnMatch",
1989      "trigger cluster #it{E} > trigger threshold, #eta vs #phi, unmatched trigger EMCAL-BC=0",
1990      netabins,etamin,etamax,nphibins,phimin,phimax);
1991     fhEtaPhiTriggerEMCALBCUMClusterOverTh->SetYTitle("#phi (rad)");
1992     fhEtaPhiTriggerEMCALBCUMClusterOverTh->SetXTitle("#eta");
1993     outputContainer->Add(fhEtaPhiTriggerEMCALBCUMClusterOverTh) ;
1994
1995     fhEtaPhiTriggerEMCALBCClusterBelowTh1 = new TH2F
1996     ("hEtaPhiTriggerEMCALBC0_OnlyTrigger_BelowThreshold1",
1997      "trigger cluster thresh-1 < #it{E} < thres, #eta vs #phi, Trigger EMCAL-BC=0",
1998      netabins,etamin,etamax,nphibins,phimin,phimax);
1999     fhEtaPhiTriggerEMCALBCClusterBelowTh1->SetYTitle("#phi (rad)");
2000     fhEtaPhiTriggerEMCALBCClusterBelowTh1->SetXTitle("#eta");
2001     outputContainer->Add(fhEtaPhiTriggerEMCALBCClusterBelowTh1) ;
2002
2003     fhEtaPhiTriggerEMCALBCUMClusterBelowTh1 = new TH2F
2004     ("hEtaPhiTriggerEMCALBC0_OnlyTrigger_BelowThreshold1_UnMatch",
2005      "trigger cluster thresh-1 < #it{E} < thres, #eta vs #phi, unmatched trigger EMCAL-BC=0",
2006      netabins,etamin,etamax,nphibins,phimin,phimax);
2007     fhEtaPhiTriggerEMCALBCUMClusterBelowTh1->SetYTitle("#phi (rad)");
2008     fhEtaPhiTriggerEMCALBCUMClusterBelowTh1->SetXTitle("#eta");
2009     outputContainer->Add(fhEtaPhiTriggerEMCALBCUMClusterBelowTh1) ;
2010
2011     fhEtaPhiTriggerEMCALBCClusterBelowTh2 = new TH2F
2012     ("hEtaPhiTriggerEMCALBC0_OnlyTrigger_BelowThreshold2",
2013      "trigger cluster thresh-2 < #it{E} < thres, #eta vs #phi, Trigger EMCAL-BC=0",
2014      netabins,etamin,etamax,nphibins,phimin,phimax);
2015     fhEtaPhiTriggerEMCALBCClusterBelowTh2->SetYTitle("#phi (rad)");
2016     fhEtaPhiTriggerEMCALBCClusterBelowTh2->SetXTitle("#eta");
2017     outputContainer->Add(fhEtaPhiTriggerEMCALBCClusterBelowTh2) ;
2018
2019     fhEtaPhiTriggerEMCALBCUMClusterBelowTh2 = new TH2F
2020     ("hEtaPhiTriggerEMCALBC0_OnlyTrigger_BelowThreshold2_UnMatch",
2021      "trigger cluster thresh-2 < #it{E} < thres, #eta vs #phi, unmatched trigger EMCAL-BC=0",
2022      netabins,etamin,etamax,nphibins,phimin,phimax);
2023     fhEtaPhiTriggerEMCALBCUMClusterBelowTh2->SetYTitle("#phi (rad)");
2024     fhEtaPhiTriggerEMCALBCUMClusterBelowTh2->SetXTitle("#eta");
2025     outputContainer->Add(fhEtaPhiTriggerEMCALBCUMClusterBelowTh2) ;
2026
2027     if(!GetReader()->AreBadTriggerEventsRemoved())
2028     {
2029       fhEtaPhiTriggerEMCALBCExotic = new TH2F
2030       ("hEtaPhiTriggerExotic",
2031        "cluster #it{E} > 2 GeV, #eta vs #phi, Trigger Exotic",
2032        netabins,etamin,etamax,nphibins,phimin,phimax);
2033       fhEtaPhiTriggerEMCALBCExotic->SetYTitle("#phi (rad)");
2034       fhEtaPhiTriggerEMCALBCExotic->SetXTitle("#eta");
2035       outputContainer->Add(fhEtaPhiTriggerEMCALBCExotic) ;
2036
2037       fhTimeTriggerEMCALBCExotic = new TH2F
2038       ("hTimeTriggerExotic",
2039        "cluster #it{time} vs #it{E} of clusters, Trigger Exotic ",
2040        nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
2041       fhTimeTriggerEMCALBCExotic->SetXTitle("#it{E} (GeV)");
2042       fhTimeTriggerEMCALBCExotic->SetYTitle("#it{time} (ns)");
2043       outputContainer->Add(fhTimeTriggerEMCALBCExotic);
2044
2045       fhEtaPhiTriggerEMCALBCUMExotic = new TH2F
2046       ("hEtaPhiTriggerExotic_UnMatch",
2047        "cluster #it{E} > 2 GeV, #eta vs #phi, unmatched trigger Exotic",
2048        netabins,etamin,etamax,nphibins,phimin,phimax);
2049       fhEtaPhiTriggerEMCALBCUMExotic->SetYTitle("#phi (rad)");
2050       fhEtaPhiTriggerEMCALBCUMExotic->SetXTitle("#eta");
2051       outputContainer->Add(fhEtaPhiTriggerEMCALBCUMExotic) ;
2052
2053       fhTimeTriggerEMCALBCUMExotic = new TH2F
2054       ("hTimeTriggerExotic_UnMatch",
2055        "cluster #it{time} vs #it{E} of clusters, unmatched trigger Exotic",
2056        nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
2057       fhTimeTriggerEMCALBCUMExotic->SetXTitle("#it{E} (GeV)");
2058       fhTimeTriggerEMCALBCUMExotic->SetYTitle("#it{time} (ns)");
2059       outputContainer->Add(fhTimeTriggerEMCALBCUMExotic);
2060
2061       fhEtaPhiTriggerEMCALBCExoticCluster = new TH2F
2062       ("hEtaPhiTriggerExotic_OnlyTrigger",
2063        "trigger cluster #it{E} > 2 GeV, #eta vs #phi, Trigger Exotic",
2064        netabins,etamin,etamax,nphibins,phimin,phimax);
2065       fhEtaPhiTriggerEMCALBCExoticCluster->SetYTitle("#phi (rad)");
2066       fhEtaPhiTriggerEMCALBCExoticCluster->SetXTitle("#eta");
2067       outputContainer->Add(fhEtaPhiTriggerEMCALBCExoticCluster) ;
2068
2069       fhTimeTriggerEMCALBCExoticCluster = new TH2F
2070       ("hTimeTriggerExotic_OnlyTrigger",
2071        "trigger cluster #it{time} vs #it{E} of clusters, Trigger Exotic",
2072        nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
2073       fhTimeTriggerEMCALBCExoticCluster->SetXTitle("#it{E} (GeV)");
2074       fhTimeTriggerEMCALBCExoticCluster->SetYTitle("#it{time} (ns)");
2075       outputContainer->Add(fhTimeTriggerEMCALBCExoticCluster);
2076
2077       fhEtaPhiTriggerEMCALBCUMExoticCluster = new TH2F
2078       ("hEtaPhiTriggerExotic_OnlyTrigger_UnMatch",
2079        "trigger cluster #it{E} > 2 GeV, #eta vs #phi, unmatched trigger Exotic",
2080        netabins,etamin,etamax,nphibins,phimin,phimax);
2081       fhEtaPhiTriggerEMCALBCUMExoticCluster->SetYTitle("#phi (rad)");
2082       fhEtaPhiTriggerEMCALBCUMExoticCluster->SetXTitle("#eta");
2083       outputContainer->Add(fhEtaPhiTriggerEMCALBCUMExoticCluster) ;
2084
2085       fhTimeTriggerEMCALBCUMExoticCluster = new TH2F
2086       ("hTimeTriggerExotic_OnlyTrigger_UnMatch",
2087        "trigger cluster #it{time} vs #it{E} of clusters, unmatched trigger Exotic",
2088        nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
2089       fhTimeTriggerEMCALBCUMExoticCluster->SetXTitle("#it{E} (GeV)");
2090       fhTimeTriggerEMCALBCUMExoticCluster->SetYTitle("#it{time} (ns)");
2091       outputContainer->Add(fhTimeTriggerEMCALBCUMExoticCluster);
2092
2093       fhEtaPhiTriggerEMCALBCBad = new TH2F
2094       ("hEtaPhiTriggerBad",
2095        "cluster #it{E} > 2 GeV, #eta vs #phi, Trigger Bad",
2096        netabins,etamin,etamax,nphibins,phimin,phimax);
2097       fhEtaPhiTriggerEMCALBCBad->SetYTitle("#phi (rad)");
2098       fhEtaPhiTriggerEMCALBCBad->SetXTitle("#eta");
2099       outputContainer->Add(fhEtaPhiTriggerEMCALBCBad) ;
2100
2101       fhTimeTriggerEMCALBCBad = new TH2F
2102       ("hTimeTriggerBad",
2103        "cluster #it{time} vs #it{E} of clusters, Trigger Bad ",
2104        nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
2105       fhTimeTriggerEMCALBCBad->SetXTitle("#it{E} (GeV)");
2106       fhTimeTriggerEMCALBCBad->SetYTitle("#it{time} (ns)");
2107       outputContainer->Add(fhTimeTriggerEMCALBCBad);
2108
2109       fhEtaPhiTriggerEMCALBCUMBad = new TH2F
2110       ("hEtaPhiTriggerBad_UnMatch",
2111        "cluster #it{E} > 2 GeV, #eta vs #phi, unmatched trigger Bad",
2112        netabins,etamin,etamax,nphibins,phimin,phimax);
2113       fhEtaPhiTriggerEMCALBCUMBad->SetYTitle("#phi (rad)");
2114       fhEtaPhiTriggerEMCALBCUMBad->SetXTitle("#eta");
2115       outputContainer->Add(fhEtaPhiTriggerEMCALBCUMBad) ;
2116
2117       fhTimeTriggerEMCALBCUMBad = new TH2F
2118       ("hTimeTriggerBad_UnMatch",
2119        "cluster #it{time} vs #it{E} of clusters, unmatched trigger Bad",
2120        nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
2121       fhTimeTriggerEMCALBCUMBad->SetXTitle("#it{E} (GeV)");
2122       fhTimeTriggerEMCALBCUMBad->SetYTitle("#it{time} (ns)");
2123       outputContainer->Add(fhTimeTriggerEMCALBCUMBad);
2124
2125       fhEtaPhiTriggerEMCALBCBadCluster = new TH2F
2126       ("hEtaPhiTriggerBad_OnlyTrigger",
2127        "trigger cluster #it{E} > 2 GeV, #eta vs #phi, Trigger Bad",
2128        netabins,etamin,etamax,nphibins,phimin,phimax);
2129       fhEtaPhiTriggerEMCALBCBadCluster->SetYTitle("#phi (rad)");
2130       fhEtaPhiTriggerEMCALBCBadCluster->SetXTitle("#eta");
2131       outputContainer->Add(fhEtaPhiTriggerEMCALBCBadCluster) ;
2132
2133       fhTimeTriggerEMCALBCBadCluster = new TH2F
2134       ("hTimeTriggerBad_OnlyTrigger",
2135        "trigger cluster #it{time} vs #it{E} of clusters, Trigger Bad",
2136        nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
2137       fhTimeTriggerEMCALBCBadCluster->SetXTitle("#it{E} (GeV)");
2138       fhTimeTriggerEMCALBCBadCluster->SetYTitle("#it{time} (ns)");
2139       outputContainer->Add(fhTimeTriggerEMCALBCBadCluster);
2140
2141       fhEtaPhiTriggerEMCALBCUMBadCluster = new TH2F
2142       ("hEtaPhiTriggerBad_OnlyTrigger_UnMatch",
2143        "trigger cluster #it{E} > 2 GeV, #eta vs #phi, unmatched trigger Bad",
2144        netabins,etamin,etamax,nphibins,phimin,phimax);
2145       fhEtaPhiTriggerEMCALBCUMBadCluster->SetYTitle("#phi (rad)");
2146       fhEtaPhiTriggerEMCALBCUMBadCluster->SetXTitle("#eta");
2147       outputContainer->Add(fhEtaPhiTriggerEMCALBCUMBadCluster) ;
2148
2149       fhTimeTriggerEMCALBCUMBadCluster = new TH2F
2150       ("hTimeTriggerBad_OnlyTrigger_UnMatch",
2151        "trigger cluster time vs #it{E} of clusters, unmatched trigger Bad",
2152        nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
2153       fhTimeTriggerEMCALBCUMBadCluster->SetXTitle("#it{E} (GeV)");
2154       fhTimeTriggerEMCALBCUMBadCluster->SetYTitle("#it{time} (ns)");
2155       outputContainer->Add(fhTimeTriggerEMCALBCUMBadCluster);
2156
2157       fhEtaPhiTriggerEMCALBCBadExotic = new TH2F
2158       ("hEtaPhiTriggerBadExotic",
2159        "cluster #it{E} > 2 GeV, #eta vs #phi, Trigger Bad&Exotic",
2160        netabins,etamin,etamax,nphibins,phimin,phimax);
2161       fhEtaPhiTriggerEMCALBCBadExotic->SetYTitle("#phi (rad)");
2162       fhEtaPhiTriggerEMCALBCBadExotic->SetXTitle("#eta");
2163       outputContainer->Add(fhEtaPhiTriggerEMCALBCBadExotic) ;
2164
2165       fhTimeTriggerEMCALBCBadExotic = new TH2F
2166       ("hTimeTriggerBadExotic",
2167        "cluster #it{time} vs #it{E} of clusters, Trigger Bad&Exotic ",
2168        nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
2169       fhTimeTriggerEMCALBCBadExotic->SetXTitle("#it{E} (GeV)");
2170       fhTimeTriggerEMCALBCBadExotic->SetYTitle("#it{time} (ns)");
2171       outputContainer->Add(fhTimeTriggerEMCALBCBadExotic);
2172
2173       fhEtaPhiTriggerEMCALBCUMBadExotic = new TH2F
2174       ("hEtaPhiTriggerBadExotic_UnMatch",
2175        "cluster #it{E} > 2 GeV, #eta vs #phi, unmatched trigger Bad&Exotic",
2176        netabins,etamin,etamax,nphibins,phimin,phimax);
2177       fhEtaPhiTriggerEMCALBCUMBadExotic->SetYTitle("#phi (rad)");
2178       fhEtaPhiTriggerEMCALBCUMBadExotic->SetXTitle("#eta");
2179       outputContainer->Add(fhEtaPhiTriggerEMCALBCUMBadExotic) ;
2180
2181       fhTimeTriggerEMCALBCUMBadExotic = new TH2F
2182       ("hTimeTriggerBadExotic_UnMatch",
2183        "cluster #it{time} vs #it{E} of clusters, unmatched trigger Bad&Exotic",
2184        nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
2185       fhTimeTriggerEMCALBCUMBadExotic->SetXTitle("#it{E} (GeV)");
2186       fhTimeTriggerEMCALBCUMBadExotic->SetYTitle("#it{time} (ns)");
2187       outputContainer->Add(fhTimeTriggerEMCALBCUMBadExotic);
2188
2189       fhEtaPhiTriggerEMCALBCBadExoticCluster = new TH2F
2190       ("hEtaPhiTriggerBadExotic_OnlyTrigger",
2191        "trigger cluster #it{E} > 2 GeV, #eta vs #phi, Trigger Bad&Exotic",
2192        netabins,etamin,etamax,nphibins,phimin,phimax);
2193       fhEtaPhiTriggerEMCALBCBadExoticCluster->SetYTitle("#phi (rad)");
2194       fhEtaPhiTriggerEMCALBCBadExoticCluster->SetXTitle("#eta");
2195       outputContainer->Add(fhEtaPhiTriggerEMCALBCBadExoticCluster) ;
2196
2197       fhTimeTriggerEMCALBCBadExoticCluster = new TH2F
2198       ("hTimeTriggerBadExotic_OnlyTrigger",
2199        "trigger cluster #it{time} vs #it{E} of clusters, Trigger Bad&Exotic",
2200        nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
2201       fhTimeTriggerEMCALBCBadExoticCluster->SetXTitle("#it{E} (GeV)");
2202       fhTimeTriggerEMCALBCBadExoticCluster->SetYTitle("#it{time} (ns)");
2203       outputContainer->Add(fhTimeTriggerEMCALBCBadExoticCluster);
2204
2205       fhEtaPhiTriggerEMCALBCUMBadExoticCluster = new TH2F
2206       ("hEtaPhiTriggerBadExotic_OnlyTrigger_UnMatch",
2207        "trigger cluster #it{E} > 2 GeV, #eta vs #phi, unmatched trigger Bad&Exotic",
2208        netabins,etamin,etamax,nphibins,phimin,phimax);
2209       fhEtaPhiTriggerEMCALBCUMBadExoticCluster->SetYTitle("#phi (rad)");
2210       fhEtaPhiTriggerEMCALBCUMBadExoticCluster->SetXTitle("#eta");
2211       outputContainer->Add(fhEtaPhiTriggerEMCALBCUMBadExoticCluster) ;
2212
2213       fhTimeTriggerEMCALBCUMBadExoticCluster = new TH2F
2214       ("hTimeTriggerBadExotic_OnlyTrigger_UnMatch",
2215        "trigger cluster #it{time} vs #it{E} of clusters, unmatched trigger Bad&Exotic",
2216        nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
2217       fhTimeTriggerEMCALBCUMBadExoticCluster->SetXTitle("#it{E} (GeV)");
2218       fhTimeTriggerEMCALBCUMBadExoticCluster->SetYTitle("#it{time} (ns)");
2219       outputContainer->Add(fhTimeTriggerEMCALBCUMBadExoticCluster);
2220
2221       fhTimeTriggerEMCALBCBadMaxCell = new TH2F
2222       ("hTimeTriggerBadMaxCell",
2223        "cluster #it{time} vs #it{E} of clusters, Trigger BadMaxCell",
2224        nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
2225       fhTimeTriggerEMCALBCBadMaxCell->SetXTitle("#it{E} (GeV)");
2226       fhTimeTriggerEMCALBCBadMaxCell->SetYTitle("#it{time} (ns)");
2227       outputContainer->Add(fhTimeTriggerEMCALBCBadMaxCell);
2228
2229       fhTimeTriggerEMCALBCUMBadMaxCell = new TH2F
2230       ("hTimeTriggerBadMaxCell_UnMatch",
2231        "cluster #it{time} vs #it{E} of clusters, unmatched trigger BadMaxCell",
2232        nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
2233       fhTimeTriggerEMCALBCUMBadMaxCell->SetXTitle("#it{E} (GeV)");
2234       fhTimeTriggerEMCALBCUMBadMaxCell->SetYTitle("#it{time} (ns)");
2235       outputContainer->Add(fhTimeTriggerEMCALBCUMBadMaxCell);
2236
2237
2238       fhTimeTriggerEMCALBCBadMaxCellExotic = new TH2F
2239       ("hTimeTriggerBadMaxCellExotic",
2240        "cluster #it{time} vs #it{E} of clusters, Trigger BadMaxCell&Exotic",
2241        nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
2242       fhTimeTriggerEMCALBCBadMaxCellExotic->SetXTitle("#it{E} (GeV)");
2243       fhTimeTriggerEMCALBCBadMaxCellExotic->SetYTitle("#it{time} (ns)");
2244       outputContainer->Add(fhTimeTriggerEMCALBCBadMaxCellExotic);
2245
2246       fhTimeTriggerEMCALBCUMBadMaxCellExotic = new TH2F
2247       ("hTimeTriggerBadMaxCellExotic_UnMatch",
2248        "cluster #it{time} vs #it{E} of clusters, unmatched trigger BadMaxCell&Exotic",
2249        nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
2250       fhTimeTriggerEMCALBCUMBadMaxCellExotic->SetXTitle("#it{E} (GeV)");
2251       fhTimeTriggerEMCALBCUMBadMaxCellExotic->SetYTitle("#it{time} (ns)");
2252       outputContainer->Add(fhTimeTriggerEMCALBCUMBadMaxCellExotic);
2253
2254       fhTimeNoTrigger = new TH2F
2255       ("hTimeNoTrigger",
2256        "events with no foundable trigger, time vs e of clusters",
2257        nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
2258       fhTimeNoTrigger->SetXTitle("#it{E} (GeV)");
2259       fhTimeNoTrigger->SetYTitle("#it{time} (ns)");
2260       outputContainer->Add(fhTimeNoTrigger);
2261
2262       fhEtaPhiNoTrigger = new TH2F
2263       ("hEtaPhiNoTrigger",
2264        "events with no foundable trigger, eta vs phi of clusters",
2265        netabins,etamin,etamax,nphibins,phimin,phimax);
2266       fhEtaPhiNoTrigger->SetYTitle("#phi (rad)");
2267       fhEtaPhiNoTrigger->SetXTitle("#eta");
2268       outputContainer->Add(fhEtaPhiNoTrigger) ;
2269     }
2270
2271     fhEtaPhiTriggerEMCALBCUMReMatchOpenTimeCluster = new TH2F("hEtaPhiTriggerEMCALBC0_OnlyTrigger_UnMatch_ReMatch_OpenTime",
2272                                                        "cluster #it{E} > 2 GeV, #eta vs #phi, Trigger EMCAL-BC=0, un match, rematch open time",
2273                                                        netabins,etamin,etamax,nphibins,phimin,phimax);
2274     fhEtaPhiTriggerEMCALBCUMReMatchOpenTimeCluster->SetYTitle("#phi (rad)");
2275     fhEtaPhiTriggerEMCALBCUMReMatchOpenTimeCluster->SetXTitle("#eta");
2276     outputContainer->Add(fhEtaPhiTriggerEMCALBCUMReMatchOpenTimeCluster) ;
2277
2278     fhTimeTriggerEMCALBCUMReMatchOpenTimeCluster = new TH2F("hTimeTrigger_OnlyTrigger_UnMatch_ReMatch_OpenTime",
2279                                                      "cluster #it{time} vs #it{E} of clusters, no match, rematch open time",
2280                                                      nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
2281     fhTimeTriggerEMCALBCUMReMatchOpenTimeCluster->SetXTitle("#it{E} (GeV)");
2282     fhTimeTriggerEMCALBCUMReMatchOpenTimeCluster->SetYTitle("#it{time} (ns)");
2283     outputContainer->Add(fhTimeTriggerEMCALBCUMReMatchOpenTimeCluster);
2284
2285
2286     fhEtaPhiTriggerEMCALBCUMReMatchCheckNeighCluster = new TH2F("hEtaPhiTriggerEMCALBC0_OnlyTrigger_UnMatch_ReMatch_CheckNeighbours",
2287                                                          "cluster #it{E} > 2 GeV, #eta vs #phi, Trigger EMCAL-BC=0, un match, rematch with neighbour patches",
2288                                                          netabins,etamin,etamax,nphibins,phimin,phimax);
2289     fhEtaPhiTriggerEMCALBCUMReMatchCheckNeighCluster->SetYTitle("#phi (rad)");
2290     fhEtaPhiTriggerEMCALBCUMReMatchCheckNeighCluster->SetXTitle("#eta");
2291     outputContainer->Add(fhEtaPhiTriggerEMCALBCUMReMatchCheckNeighCluster) ;
2292
2293     fhTimeTriggerEMCALBCUMReMatchCheckNeighCluster = new TH2F("hTimeTrigger_OnlyTrigger_UnMatch_ReMatch_CheckNeighbours",
2294                                                        "cluster #it{time} vs #it{E} of clusters, no match, rematch with neigbour parches",
2295                                                        nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
2296     fhTimeTriggerEMCALBCUMReMatchCheckNeighCluster->SetXTitle("#it{E} (GeV)");
2297     fhTimeTriggerEMCALBCUMReMatchCheckNeighCluster->SetYTitle("#it{time} (ns)");
2298     outputContainer->Add(fhTimeTriggerEMCALBCUMReMatchCheckNeighCluster);
2299
2300     fhEtaPhiTriggerEMCALBCUMReMatchBothCluster = new TH2F("hEtaPhiTriggerEMCALBC0_OnlyTrigger_UnMatch_ReMatch_Both",
2301                                                    "cluster #it{E} > 2 GeV, #eta vs #phi, Trigger EMCAL-BC=0, un match, rematch open time and neighbour",
2302                                                    netabins,etamin,etamax,nphibins,phimin,phimax);
2303     fhEtaPhiTriggerEMCALBCUMReMatchBothCluster->SetYTitle("#phi (rad)");
2304     fhEtaPhiTriggerEMCALBCUMReMatchBothCluster->SetXTitle("#eta");
2305     outputContainer->Add(fhEtaPhiTriggerEMCALBCUMReMatchBothCluster) ;
2306
2307     fhTimeTriggerEMCALBCUMReMatchBothCluster = new TH2F("hTimeTrigger_OnlyTrigger_UnMatch_ReMatch_Both",
2308                                                  "cluster #it{time} vs #it{E} of clusters, no match, rematch open time and neigbour",
2309                                                  nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
2310     fhTimeTriggerEMCALBCUMReMatchBothCluster->SetXTitle("#it{E} (GeV)");
2311     fhTimeTriggerEMCALBCUMReMatchBothCluster->SetYTitle("#it{time} (ns)");
2312     outputContainer->Add(fhTimeTriggerEMCALBCUMReMatchBothCluster);
2313
2314     fhTimeTriggerEMCALBC0UMReMatchOpenTime = new TH2F("hTimeTriggerBC0_UnMatch_ReMatch_OpenTime",
2315                                                             "cluster #it{time} vs #it{E} of clusters, no match, rematch open time",
2316                                                             nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
2317     fhTimeTriggerEMCALBC0UMReMatchOpenTime->SetXTitle("#it{E} (GeV)");
2318     fhTimeTriggerEMCALBC0UMReMatchOpenTime->SetYTitle("#it{time} (ns)");
2319     outputContainer->Add(fhTimeTriggerEMCALBC0UMReMatchOpenTime);
2320
2321
2322     fhTimeTriggerEMCALBC0UMReMatchCheckNeigh = new TH2F("hTimeTriggerBC0_UnMatch_ReMatch_CheckNeighbours",
2323                                                               "cluster #it{time} vs #it{E} of clusters, no match, rematch with neigbour parches",
2324                                                               nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
2325     fhTimeTriggerEMCALBC0UMReMatchCheckNeigh->SetXTitle("#it{E} (GeV)");
2326     fhTimeTriggerEMCALBC0UMReMatchCheckNeigh->SetYTitle("#it{time} (ns)");
2327     outputContainer->Add(fhTimeTriggerEMCALBC0UMReMatchCheckNeigh);
2328
2329     fhTimeTriggerEMCALBC0UMReMatchBoth = new TH2F("hTimeTriggerBC0_UnMatch_ReMatch_Both",
2330                                                         "cluster #it{time} vs #it{E} of clusters, no match, rematch open time and neigbour",
2331                                                         nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
2332     fhTimeTriggerEMCALBC0UMReMatchBoth->SetXTitle("#it{E} (GeV)");
2333     fhTimeTriggerEMCALBC0UMReMatchBoth->SetYTitle("#it{time} (ns)");
2334     outputContainer->Add(fhTimeTriggerEMCALBC0UMReMatchBoth);
2335
2336   }
2337
2338   fhPhiPhoton  = new TH2F
2339   ("hPhiPhoton","#phi_{#gamma} vs #it{p}_[T]",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
2340   fhPhiPhoton->SetYTitle("#phi (rad)");
2341   fhPhiPhoton->SetXTitle("p_{T #gamma} (GeV/#it{c})");
2342   outputContainer->Add(fhPhiPhoton) ;
2343
2344   fhEtaPhoton  = new TH2F
2345   ("hEtaPhoton","#eta_{#gamma} vs #it{p}_[T]",nptbins,ptmin,ptmax,netabins,etamin,etamax);
2346   fhEtaPhoton->SetYTitle("#eta");
2347   fhEtaPhoton->SetXTitle("p_{T #gamma} (GeV/#it{c})");
2348   outputContainer->Add(fhEtaPhoton) ;
2349
2350   fhEtaPhiPhoton  = new TH2F
2351   ("hEtaPhiPhoton","#eta vs #phi",netabins,etamin,etamax,nphibins,phimin,phimax);
2352   fhEtaPhiPhoton->SetYTitle("#phi (rad)");
2353   fhEtaPhiPhoton->SetXTitle("#eta");
2354   outputContainer->Add(fhEtaPhiPhoton) ;
2355   if(GetMinPt() < 0.5)
2356   {
2357     fhEtaPhi05Photon  = new TH2F
2358     ("hEtaPhi05Photon","#eta vs #phi, E < 0.5",netabins,etamin,etamax,nphibins,phimin,phimax);
2359     fhEtaPhi05Photon->SetYTitle("#phi (rad)");
2360     fhEtaPhi05Photon->SetXTitle("#eta");
2361     outputContainer->Add(fhEtaPhi05Photon) ;
2362   }
2363
2364   if(fCalorimeter=="EMCAL" && fFillEMCALBCHistograms)
2365   {
2366     fhEtaPhiPhotonEMCALBC0  = new TH2F
2367     ("hEtaPhiPhotonEMCALBC0","identified photon, #it{E} > 2 GeV, #eta vs #phi, for clusters with |time| < 25 ns, EMCAL-BC=0",netabins,etamin,etamax,nphibins,phimin,phimax);
2368     fhEtaPhiPhotonEMCALBC0->SetYTitle("#phi (rad)");
2369     fhEtaPhiPhotonEMCALBC0->SetXTitle("#eta");
2370     outputContainer->Add(fhEtaPhiPhotonEMCALBC0) ;
2371
2372     fhEtaPhiPhotonEMCALBC1  = new TH2F
2373     ("hEtaPhiPhotonEMCALBC1","identified photon, #it{E} > 2 GeV, #eta vs #phi, for clusters with 25 < |time| < 75 ns, EMCAL-BC=1",netabins,etamin,etamax,nphibins,phimin,phimax);
2374     fhEtaPhiPhotonEMCALBC1->SetYTitle("#phi (rad)");
2375     fhEtaPhiPhotonEMCALBC1->SetXTitle("#eta");
2376     outputContainer->Add(fhEtaPhiPhotonEMCALBC1) ;
2377
2378     fhEtaPhiPhotonEMCALBCN  = new TH2F
2379     ("hEtaPhiPhotonEMCALBCN","identified photon, #it{E} > 2 GeV, #eta vs #phi, for clusters with |time| > 75 ns, EMCAL-BC>1",netabins,etamin,etamax,nphibins,phimin,phimax);
2380     fhEtaPhiPhotonEMCALBCN->SetYTitle("#phi (rad)");
2381     fhEtaPhiPhotonEMCALBCN->SetXTitle("#eta");
2382     outputContainer->Add(fhEtaPhiPhotonEMCALBCN) ;
2383
2384     for(Int_t i = 0; i < nTrigBC; i++)
2385     {
2386       fhEtaPhiPhotonTriggerEMCALBC[i] = new TH2F
2387       (Form("hEtaPhiPhotonTriggerEMCALBC%d",i-iBCShift),
2388        Form("photon #it{E} > 2 GeV, #eta vs #phi, PhotonTrigger EMCAL-BC=%d",i-iBCShift),
2389        netabins,etamin,etamax,nphibins,phimin,phimax);
2390       fhEtaPhiPhotonTriggerEMCALBC[i]->SetYTitle("#phi (rad)");
2391       fhEtaPhiPhotonTriggerEMCALBC[i]->SetXTitle("#eta");
2392       outputContainer->Add(fhEtaPhiPhotonTriggerEMCALBC[i]) ;
2393
2394       fhTimePhotonTriggerEMCALBC[i] = new TH2F
2395       (Form("hTimePhotonTriggerEMCALBC%d",i-iBCShift),
2396        Form("photon #it{time} vs #it{E} of clusters, PhotonTrigger EMCAL-BC=%d",i-iBCShift),
2397        nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
2398       fhTimePhotonTriggerEMCALBC[i]->SetXTitle("#it{E} (GeV)");
2399       fhTimePhotonTriggerEMCALBC[i]->SetYTitle("#it{time} (ns)");
2400       outputContainer->Add(fhTimePhotonTriggerEMCALBC[i]);
2401
2402       fhTimePhotonTriggerEMCALBCPileUpSPD[i] = new TH2F
2403       (Form("hTimePhotonTriggerEMCALBC%dPileUpSPD",i-iBCShift),
2404        Form("photon #it{time} vs #it{E}, PhotonTrigger EMCAL-BC=%d",i-iBCShift),
2405        nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
2406       fhTimePhotonTriggerEMCALBCPileUpSPD[i]->SetXTitle("#it{E} (GeV)");
2407       fhTimePhotonTriggerEMCALBCPileUpSPD[i]->SetYTitle("#it{time} (ns)");
2408       outputContainer->Add(fhTimePhotonTriggerEMCALBCPileUpSPD[i]);
2409
2410       fhEtaPhiPhotonTriggerEMCALBCUM[i] = new TH2F
2411       (Form("hEtaPhiPhotonTriggerEMCALBC%d_UnMatch",i-iBCShift),
2412        Form("photon #it{E} > 2 GeV, #eta vs #phi, unmatched trigger EMCAL-BC=%d",i-iBCShift),
2413        netabins,etamin,etamax,nphibins,phimin,phimax);
2414       fhEtaPhiPhotonTriggerEMCALBCUM[i]->SetYTitle("#phi (rad)");
2415       fhEtaPhiPhotonTriggerEMCALBCUM[i]->SetXTitle("#eta");
2416       outputContainer->Add(fhEtaPhiPhotonTriggerEMCALBCUM[i]) ;
2417
2418       fhTimePhotonTriggerEMCALBCUM[i] = new TH2F
2419       (Form("hTimePhotonTriggerEMCALBC%d_UnMatch",i-iBCShift),
2420        Form("photon #it{time} vs #it{E}, unmatched trigger EMCAL-BC=%d",i-iBCShift),
2421        nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
2422       fhTimePhotonTriggerEMCALBCUM[i]->SetXTitle("#it{E} (GeV)");
2423       fhTimePhotonTriggerEMCALBCUM[i]->SetYTitle("#it{time} (ns)");
2424       outputContainer->Add(fhTimePhotonTriggerEMCALBCUM[i]);
2425
2426     }
2427
2428     fhTimePhotonTriggerEMCALBC0UMReMatchOpenTime = new TH2F("hTimePhotonTriggerBC0_UnMatch_ReMatch_OpenTime",
2429                                                       "cluster #it{time} vs #it{E} of photons, no match, rematch open time",
2430                                                       nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
2431     fhTimePhotonTriggerEMCALBC0UMReMatchOpenTime->SetXTitle("#it{E} (GeV)");
2432     fhTimePhotonTriggerEMCALBC0UMReMatchOpenTime->SetYTitle("#it{time} (ns)");
2433     outputContainer->Add(fhTimePhotonTriggerEMCALBC0UMReMatchOpenTime);
2434
2435
2436     fhTimePhotonTriggerEMCALBC0UMReMatchCheckNeigh = new TH2F("hTimePhotonTriggerBC0_UnMatch_ReMatch_CheckNeighbours",
2437                                                         "cluster #it{time} vs #it{E} of photons, no match, rematch with neigbour parches",
2438                                                         nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
2439     fhTimePhotonTriggerEMCALBC0UMReMatchCheckNeigh->SetXTitle("#it{E} (GeV)");
2440     fhTimePhotonTriggerEMCALBC0UMReMatchCheckNeigh->SetYTitle("#it{time} (ns)");
2441     outputContainer->Add(fhTimePhotonTriggerEMCALBC0UMReMatchCheckNeigh);
2442
2443     fhTimePhotonTriggerEMCALBC0UMReMatchBoth = new TH2F("hTimePhotonTriggerBC0_UnMatch_ReMatch_Both",
2444                                                   "cluster #it{time} vs #it{E} of photons, no match, rematch open time and neigbour",
2445                                                   nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
2446     fhTimePhotonTriggerEMCALBC0UMReMatchBoth->SetXTitle("#it{E} (GeV)");
2447     fhTimePhotonTriggerEMCALBC0UMReMatchBoth->SetYTitle("#it{time} (ns)");
2448     outputContainer->Add(fhTimePhotonTriggerEMCALBC0UMReMatchBoth);
2449
2450   }
2451
2452   fhNLocMax = new TH2F("hNLocMax","Number of local maxima in cluster",
2453                        nptbins,ptmin,ptmax,10,0,10);
2454   fhNLocMax ->SetYTitle("N maxima");
2455   fhNLocMax ->SetXTitle("#it{E} (GeV)");
2456   outputContainer->Add(fhNLocMax) ;
2457
2458   //Shower shape
2459   if(fFillSSHistograms)
2460   {
2461     fhLam0E  = new TH2F ("hLam0E","#lambda_{0}^{2} vs E", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2462     fhLam0E->SetYTitle("#lambda_{0}^{2}");
2463     fhLam0E->SetXTitle("#it{E} (GeV)");
2464     outputContainer->Add(fhLam0E);
2465
2466     fhLam1E  = new TH2F ("hLam1E","#lambda_{1}^{2} vs E", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2467     fhLam1E->SetYTitle("#lambda_{1}^{2}");
2468     fhLam1E->SetXTitle("#it{E} (GeV)");
2469     outputContainer->Add(fhLam1E);
2470
2471     fhDispE  = new TH2F ("hDispE"," dispersion^{2} vs E", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2472     fhDispE->SetYTitle("D^{2}");
2473     fhDispE->SetXTitle("#it{E} (GeV) ");
2474     outputContainer->Add(fhDispE);
2475
2476     if(!fRejectTrackMatch)
2477     {
2478       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);
2479       fhLam0ETM->SetYTitle("#lambda_{0}^{2}");
2480       fhLam0ETM->SetXTitle("#it{E} (GeV)");
2481       outputContainer->Add(fhLam0ETM);
2482
2483       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);
2484       fhLam1ETM->SetYTitle("#lambda_{1}^{2}");
2485       fhLam1ETM->SetXTitle("#it{E} (GeV)");
2486       outputContainer->Add(fhLam1ETM);
2487
2488       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);
2489       fhDispETM->SetYTitle("D^{2}");
2490       fhDispETM->SetXTitle("#it{E} (GeV) ");
2491       outputContainer->Add(fhDispETM);
2492     }
2493
2494     if(fCalorimeter == "EMCAL")
2495     {
2496       fhLam0ETRD  = new TH2F ("hLam0ETRD","#lambda_{0}^{2} vs E, EMCAL SM covered by TRD", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2497       fhLam0ETRD->SetYTitle("#lambda_{0}^{2}");
2498       fhLam0ETRD->SetXTitle("#it{E} (GeV)");
2499       outputContainer->Add(fhLam0ETRD);
2500
2501       fhLam1ETRD  = new TH2F ("hLam1ETRD","#lambda_{1}^{2} vs E, EMCAL SM covered by TRD", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2502       fhLam1ETRD->SetYTitle("#lambda_{1}^{2}");
2503       fhLam1ETRD->SetXTitle("#it{E} (GeV)");
2504       outputContainer->Add(fhLam1ETRD);
2505
2506       fhDispETRD  = new TH2F ("hDispETRD"," dispersion^{2} vs E, EMCAL SM covered by TRD", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2507       fhDispETRD->SetYTitle("Dispersion^{2}");
2508       fhDispETRD->SetXTitle("#it{E} (GeV) ");
2509       outputContainer->Add(fhDispETRD);
2510
2511       if(!fRejectTrackMatch)
2512       {
2513         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);
2514         fhLam0ETMTRD->SetYTitle("#lambda_{0}^{2}");
2515         fhLam0ETMTRD->SetXTitle("#it{E} (GeV)");
2516         outputContainer->Add(fhLam0ETMTRD);
2517
2518         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);
2519         fhLam1ETMTRD->SetYTitle("#lambda_{1}^{2}");
2520         fhLam1ETMTRD->SetXTitle("#it{E} (GeV)");
2521         outputContainer->Add(fhLam1ETMTRD);
2522
2523         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);
2524         fhDispETMTRD->SetYTitle("Dispersion^{2}");
2525         fhDispETMTRD->SetXTitle("#it{E} (GeV) ");
2526         outputContainer->Add(fhDispETMTRD);
2527       }
2528     }
2529
2530     if(!fFillOnlySimpleSSHisto)
2531     {
2532       fhNCellsLam0LowE  = new TH2F ("hNCellsLam0LowE","N_{cells} in cluster vs #lambda_{0}^{2}, E < 2 GeV", nbins,nmin, nmax, ssbins,ssmin,ssmax);
2533       fhNCellsLam0LowE->SetXTitle("N_{cells}");
2534       fhNCellsLam0LowE->SetYTitle("#lambda_{0}^{2}");
2535       outputContainer->Add(fhNCellsLam0LowE);
2536
2537       fhNCellsLam0HighE  = new TH2F ("hNCellsLam0HighE","N_{cells} in cluster vs #lambda_{0}^{2}, #it{E} > 2 GeV", nbins,nmin, nmax, ssbins,ssmin,ssmax);
2538       fhNCellsLam0HighE->SetXTitle("N_{cells}");
2539       fhNCellsLam0HighE->SetYTitle("#lambda_{0}^{2}");
2540       outputContainer->Add(fhNCellsLam0HighE);
2541
2542       fhNCellsLam1LowE  = new TH2F ("hNCellsLam1LowE","N_{cells} in cluster vs #lambda_{1}^{2}, E < 2 GeV", nbins,nmin, nmax, ssbins,ssmin,ssmax);
2543       fhNCellsLam1LowE->SetXTitle("N_{cells}");
2544       fhNCellsLam1LowE->SetYTitle("#lambda_{0}^{2}");
2545       outputContainer->Add(fhNCellsLam1LowE);
2546
2547       fhNCellsLam1HighE  = new TH2F ("hNCellsLam1HighE","N_{cells} in cluster vs #lambda_{1}^{2}, #it{E} > 2 GeV", nbins,nmin, nmax, ssbins,ssmin,ssmax);
2548       fhNCellsLam1HighE->SetXTitle("N_{cells}");
2549       fhNCellsLam1HighE->SetYTitle("#lambda_{0}^{2}");
2550       outputContainer->Add(fhNCellsLam1HighE);
2551
2552       fhNCellsDispLowE  = new TH2F ("hNCellsDispLowE","N_{cells} in cluster vs dispersion^{2}, E < 2 GeV", nbins,nmin, nmax, ssbins,ssmin,ssmax);
2553       fhNCellsDispLowE->SetXTitle("N_{cells}");
2554       fhNCellsDispLowE->SetYTitle("D^{2}");
2555       outputContainer->Add(fhNCellsDispLowE);
2556
2557       fhNCellsDispHighE  = new TH2F ("hNCellsDispHighE","N_{cells} in cluster vs dispersion^{2}, E < 2 GeV", nbins,nmin, nmax, ssbins,ssmin,ssmax);
2558       fhNCellsDispHighE->SetXTitle("N_{cells}");
2559       fhNCellsDispHighE->SetYTitle("D^{2}");
2560       outputContainer->Add(fhNCellsDispHighE);
2561
2562       fhEtaLam0LowE  = new TH2F ("hEtaLam0LowE","#eta vs #lambda_{0}^{2}, E < 2 GeV", netabins,etamin,etamax, ssbins,ssmin,ssmax);
2563       fhEtaLam0LowE->SetYTitle("#lambda_{0}^{2}");
2564       fhEtaLam0LowE->SetXTitle("#eta");
2565       outputContainer->Add(fhEtaLam0LowE);
2566
2567       fhPhiLam0LowE  = new TH2F ("hPhiLam0LowE","#phi vs #lambda_{0}^{2}, E < 2 GeV", nphibins,phimin,phimax, ssbins,ssmin,ssmax);
2568       fhPhiLam0LowE->SetYTitle("#lambda_{0}^{2}");
2569       fhPhiLam0LowE->SetXTitle("#phi");
2570       outputContainer->Add(fhPhiLam0LowE);
2571
2572       fhEtaLam0HighE  = new TH2F ("hEtaLam0HighE","#eta vs #lambda_{0}^{2}, #it{E} > 2 GeV", netabins,etamin,etamax, ssbins,ssmin,ssmax);
2573       fhEtaLam0HighE->SetYTitle("#lambda_{0}^{2}");
2574       fhEtaLam0HighE->SetXTitle("#eta");
2575       outputContainer->Add(fhEtaLam0HighE);
2576
2577       fhPhiLam0HighE  = new TH2F ("hPhiLam0HighE","#phi vs #lambda_{0}^{2}, #it{E} > 2 GeV", nphibins,phimin,phimax, ssbins,ssmin,ssmax);
2578       fhPhiLam0HighE->SetYTitle("#lambda_{0}^{2}");
2579       fhPhiLam0HighE->SetXTitle("#phi");
2580       outputContainer->Add(fhPhiLam0HighE);
2581
2582       fhLam1Lam0LowE  = new TH2F ("hLam1Lam0LowE","#lambda_{0}^{2} vs #lambda_{1}^{2} in cluster of E < 2 GeV",  ssbins,ssmin,ssmax, ssbins,ssmin,ssmax);
2583       fhLam1Lam0LowE->SetYTitle("#lambda_{0}^{2}");
2584       fhLam1Lam0LowE->SetXTitle("#lambda_{1}^{2}");
2585       outputContainer->Add(fhLam1Lam0LowE);
2586
2587       fhLam1Lam0HighE  = new TH2F ("hLam1Lam0HighE","#lambda_{0}^{2} vs #lambda_{1}^{2} in cluster of #it{E} > 2 GeV",  ssbins,ssmin,ssmax, ssbins,ssmin,ssmax);
2588       fhLam1Lam0HighE->SetYTitle("#lambda_{0}^{2}");
2589       fhLam1Lam0HighE->SetXTitle("#lambda_{1}^{2}");
2590       outputContainer->Add(fhLam1Lam0HighE);
2591
2592       fhLam0DispLowE  = new TH2F ("hLam0DispLowE","#lambda_{0}^{2} vs dispersion^{2} in cluster of E < 2 GeV",  ssbins,ssmin,ssmax, ssbins,ssmin,ssmax);
2593       fhLam0DispLowE->SetXTitle("#lambda_{0}^{2}");
2594       fhLam0DispLowE->SetYTitle("D^{2}");
2595       outputContainer->Add(fhLam0DispLowE);
2596
2597       fhLam0DispHighE  = new TH2F ("hLam0DispHighE","#lambda_{0}^{2} vs dispersion^{2} in cluster of #it{E} > 2 GeV",  ssbins,ssmin,ssmax, ssbins,ssmin,ssmax);
2598       fhLam0DispHighE->SetXTitle("#lambda_{0}^{2}");
2599       fhLam0DispHighE->SetYTitle("D^{2}");
2600       outputContainer->Add(fhLam0DispHighE);
2601
2602       fhDispLam1LowE  = new TH2F ("hDispLam1LowE","Dispersion^{2} vs #lambda_{1}^{2} in cluster of E < 2 GeV",  ssbins,ssmin,ssmax, ssbins,ssmin,ssmax);
2603       fhDispLam1LowE->SetXTitle("D^{2}");
2604       fhDispLam1LowE->SetYTitle("#lambda_{1}^{2}");
2605       outputContainer->Add(fhDispLam1LowE);
2606
2607       fhDispLam1HighE  = new TH2F ("hDispLam1HighE","Dispersion^{2} vs #lambda_{1^{2}} in cluster of #it{E} > 2 GeV",  ssbins,ssmin,ssmax, ssbins,ssmin,ssmax);
2608       fhDispLam1HighE->SetXTitle("D^{2}");
2609       fhDispLam1HighE->SetYTitle("#lambda_{1}^{2}");
2610       outputContainer->Add(fhDispLam1HighE);
2611
2612       if(fCalorimeter == "EMCAL")
2613       {
2614         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);
2615         fhDispEtaE->SetXTitle("#it{E} (GeV)");
2616         fhDispEtaE->SetYTitle("#sigma^{2}_{#eta #eta}");
2617         outputContainer->Add(fhDispEtaE);
2618
2619         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);
2620         fhDispPhiE->SetXTitle("#it{E} (GeV)");
2621         fhDispPhiE->SetYTitle("#sigma^{2}_{#phi #phi}");
2622         outputContainer->Add(fhDispPhiE);
2623
2624         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);
2625         fhSumEtaE->SetXTitle("#it{E} (GeV)");
2626         fhSumEtaE->SetYTitle("#delta^{2}_{#eta #eta}");
2627         outputContainer->Add(fhSumEtaE);
2628
2629         fhSumPhiE  = new TH2F ("hSumPhiE","#delta^{2}_{#phi #phi} = #Sigma w_{i}(#phi_{i})^{2}/ #Sigma w_{i} - <#phi>^{2} vs E",
2630                                nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
2631         fhSumPhiE->SetXTitle("#it{E} (GeV)");
2632         fhSumPhiE->SetYTitle("#delta^{2}_{#phi #phi}");
2633         outputContainer->Add(fhSumPhiE);
2634
2635         fhSumEtaPhiE  = new TH2F ("hSumEtaPhiE","#delta^{2}_{#eta #phi} = #Sigma w_{i}(#phi_{i} #eta_{i} ) / #Sigma w_{i} - <#phi><#eta> vs E",
2636                                   nptbins,ptmin,ptmax, 2*ssbins,-ssmax,ssmax);
2637         fhSumEtaPhiE->SetXTitle("#it{E} (GeV)");
2638         fhSumEtaPhiE->SetYTitle("#delta^{2}_{#eta #phi}");
2639         outputContainer->Add(fhSumEtaPhiE);
2640
2641         fhDispEtaPhiDiffE  = new TH2F ("hDispEtaPhiDiffE","#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta} vs E",
2642                                        nptbins,ptmin,ptmax,200, -10,10);
2643         fhDispEtaPhiDiffE->SetXTitle("#it{E} (GeV)");
2644         fhDispEtaPhiDiffE->SetYTitle("#sigma^{2}_{#phi #phi}-#sigma^{2}_{#eta #eta}");
2645         outputContainer->Add(fhDispEtaPhiDiffE);
2646
2647         fhSphericityE  = new TH2F ("hSphericityE","(#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi}) vs E",
2648                                    nptbins,ptmin,ptmax, 200, -1,1);
2649         fhSphericityE->SetXTitle("#it{E} (GeV)");
2650         fhSphericityE->SetYTitle("s = (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi})");
2651         outputContainer->Add(fhSphericityE);
2652
2653         fhDispSumEtaDiffE  = new TH2F ("hDispSumEtaDiffE","#sigma^{2}_{#eta #eta} - #delta^{2}_{#eta #eta} / average vs E",  nptbins,ptmin,ptmax, 200,-0.01,0.01);
2654         fhDispSumEtaDiffE->SetXTitle("#it{E} (GeV)");
2655         fhDispSumEtaDiffE->SetYTitle("#sigma^{2}_{#eta #eta} - #delta^{2}_{#eta #eta} / average");
2656         outputContainer->Add(fhDispSumEtaDiffE);
2657
2658         fhDispSumPhiDiffE  = new TH2F ("hDispSumPhiDiffE","#sigma^{2}_{#phi #phi} - #delta^{2}_{#phi #phi} / average vs E",  nptbins,ptmin,ptmax, 200,-0.01,0.01);
2659         fhDispSumPhiDiffE->SetXTitle("#it{E} (GeV)");
2660         fhDispSumPhiDiffE->SetYTitle("#sigma^{2}_{#phi #phi} - #delta^{2}_{#phi #phi} / average");
2661         outputContainer->Add(fhDispSumPhiDiffE);
2662
2663         for(Int_t i = 0; i < 7; i++)
2664         {
2665           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]),
2666                                           ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
2667           fhDispEtaDispPhi[i]->SetXTitle("#sigma^{2}_{#eta #eta}");
2668           fhDispEtaDispPhi[i]->SetYTitle("#sigma^{2}_{#phi #phi}");
2669           outputContainer->Add(fhDispEtaDispPhi[i]);
2670
2671           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]),
2672                                           ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
2673           fhLambda0DispEta[i]->SetXTitle("#lambda^{2}_{0}");
2674           fhLambda0DispEta[i]->SetYTitle("#sigma^{2}_{#eta #eta}");
2675           outputContainer->Add(fhLambda0DispEta[i]);
2676
2677           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]),
2678                                           ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
2679           fhLambda0DispPhi[i]->SetXTitle("#lambda^{2}_{0}");
2680           fhLambda0DispPhi[i]->SetYTitle("#sigma^{2}_{#phi #phi}");
2681           outputContainer->Add(fhLambda0DispPhi[i]);
2682         }
2683       }
2684     }
2685   } // Shower shape
2686
2687   // Track Matching
2688
2689   if(fFillTMHisto)
2690   {
2691     TString cutTM [] = {"NoCut",""};
2692
2693     for(Int_t i = 0; i < 2; i++)
2694     {
2695       fhTrackMatchedDEta[i]  = new TH2F
2696       (Form("hTrackMatchedDEta%s",cutTM[i].Data()),
2697        Form("d#eta of cluster-track vs cluster energy, %s",cutTM[i].Data()),
2698        nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
2699       fhTrackMatchedDEta[i]->SetYTitle("d#eta");
2700       fhTrackMatchedDEta[i]->SetXTitle("#it{E}_{cluster} (GeV)");
2701
2702       fhTrackMatchedDPhi[i]  = new TH2F
2703       (Form("hTrackMatchedDPhi%s",cutTM[i].Data()),
2704        Form("d#phi of cluster-track vs cluster energy, %s",cutTM[i].Data()),
2705        nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
2706       fhTrackMatchedDPhi[i]->SetYTitle("d#phi (rad)");
2707       fhTrackMatchedDPhi[i]->SetXTitle("#it{E}_{cluster} (GeV)");
2708
2709       fhTrackMatchedDEtaDPhi[i]  = new TH2F
2710       (Form("hTrackMatchedDEtaDPhi%s",cutTM[i].Data()),
2711        Form("d#eta vs d#phi of cluster-track vs cluster energy, %s",cutTM[i].Data()),
2712        nresetabins,resetamin,resetamax,nresphibins,resphimin,resphimax);
2713       fhTrackMatchedDEtaDPhi[i]->SetYTitle("d#phi (rad)");
2714       fhTrackMatchedDEtaDPhi[i]->SetXTitle("d#eta");
2715
2716       fhTrackMatchedDEtaPos[i]  = new TH2F
2717       (Form("hTrackMatchedDEtaPos%s",cutTM[i].Data()),
2718        Form("d#eta of cluster-track vs cluster energy, %s",cutTM[i].Data()),
2719        nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
2720       fhTrackMatchedDEtaPos[i]->SetYTitle("d#eta");
2721       fhTrackMatchedDEtaPos[i]->SetXTitle("#it{E}_{cluster} (GeV)");
2722
2723       fhTrackMatchedDPhiPos[i]  = new TH2F
2724       (Form("hTrackMatchedDPhiPos%s",cutTM[i].Data()),
2725        Form("d#phi of cluster-track vs cluster energy, %s",cutTM[i].Data()),
2726        nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
2727       fhTrackMatchedDPhiPos[i]->SetYTitle("d#phi (rad)");
2728       fhTrackMatchedDPhiPos[i]->SetXTitle("#it{E}_{cluster} (GeV)");
2729
2730       fhTrackMatchedDEtaDPhiPos[i]  = new TH2F
2731       (Form("hTrackMatchedDEtaDPhiPos%s",cutTM[i].Data()),
2732        Form("d#eta vs d#phi of cluster-track vs cluster energy, %s",cutTM[i].Data()),
2733        nresetabins,resetamin,resetamax,nresphibins,resphimin,resphimax);
2734       fhTrackMatchedDEtaDPhiPos[i]->SetYTitle("d#phi (rad)");
2735       fhTrackMatchedDEtaDPhiPos[i]->SetXTitle("d#eta");
2736
2737       fhTrackMatchedDEtaNeg[i]  = new TH2F
2738       (Form("hTrackMatchedDEtaNeg%s",cutTM[i].Data()),
2739        Form("d#eta of cluster-track vs cluster energy, %s",cutTM[i].Data()),
2740        nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
2741       fhTrackMatchedDEtaNeg[i]->SetYTitle("d#eta");
2742       fhTrackMatchedDEtaNeg[i]->SetXTitle("#it{E}_{cluster} (GeV)");
2743
2744       fhTrackMatchedDPhiNeg[i]  = new TH2F
2745       (Form("hTrackMatchedDPhiNeg%s",cutTM[i].Data()),
2746        Form("d#phi of cluster-track vs cluster energy, %s",cutTM[i].Data()),
2747        nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
2748       fhTrackMatchedDPhiNeg[i]->SetYTitle("d#phi (rad)");
2749       fhTrackMatchedDPhiNeg[i]->SetXTitle("#it{E}_{cluster} (GeV)");
2750
2751       fhTrackMatchedDEtaDPhiNeg[i]  = new TH2F
2752       (Form("hTrackMatchedDEtaDPhiNeg%s",cutTM[i].Data()),
2753        Form("d#eta vs d#phi of cluster-track vs cluster energy, %s",cutTM[i].Data()),
2754        nresetabins,resetamin,resetamax,nresphibins,resphimin,resphimax);
2755       fhTrackMatchedDEtaDPhiNeg[i]->SetYTitle("d#phi (rad)");
2756       fhTrackMatchedDEtaDPhiNeg[i]->SetXTitle("d#eta");
2757
2758       fhdEdx[i]  = new TH2F (Form("hdEdx%s",cutTM[i].Data()),Form("matched track <dE/dx> vs cluster E, %s",cutTM[i].Data()),
2759                              nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
2760       fhdEdx[i]->SetXTitle("#it{E} (GeV)");
2761       fhdEdx[i]->SetYTitle("<dE/dx>");
2762
2763       fhEOverP[i]  = new TH2F (Form("hEOverP%s",cutTM[i].Data()),Form("matched track E/p vs cluster E, %s",cutTM[i].Data()),
2764                                nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
2765       fhEOverP[i]->SetXTitle("#it{E} (GeV)");
2766       fhEOverP[i]->SetYTitle("E/p");
2767
2768       outputContainer->Add(fhTrackMatchedDEta[i]) ;
2769       outputContainer->Add(fhTrackMatchedDPhi[i]) ;
2770       outputContainer->Add(fhTrackMatchedDEtaDPhi[i]) ;
2771       outputContainer->Add(fhTrackMatchedDEtaPos[i]) ;
2772       outputContainer->Add(fhTrackMatchedDPhiPos[i]) ;
2773       outputContainer->Add(fhTrackMatchedDEtaDPhiPos[i]) ;
2774       outputContainer->Add(fhTrackMatchedDEtaNeg[i]) ;
2775       outputContainer->Add(fhTrackMatchedDPhiNeg[i]) ;
2776       outputContainer->Add(fhTrackMatchedDEtaDPhiNeg[i]) ;
2777       outputContainer->Add(fhdEdx[i]);
2778       outputContainer->Add(fhEOverP[i]);
2779
2780       if(fCalorimeter=="EMCAL")
2781       {
2782         fhTrackMatchedDEtaTRD[i]  = new TH2F
2783         (Form("hTrackMatchedDEtaTRD%s",cutTM[i].Data()),
2784          Form("d#eta of cluster-track vs cluster energy, SM behind TRD, %s",cutTM[i].Data()),
2785          nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
2786         fhTrackMatchedDEtaTRD[i]->SetYTitle("d#eta");
2787         fhTrackMatchedDEtaTRD[i]->SetXTitle("#it{E}_{cluster} (GeV)");
2788
2789         fhTrackMatchedDPhiTRD[i]  = new TH2F
2790         (Form("hTrackMatchedDPhiTRD%s",cutTM[i].Data()),
2791          Form("d#phi of cluster-track vs cluster energy, SM behing TRD, %s",cutTM[i].Data()),
2792          nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
2793         fhTrackMatchedDPhiTRD[i]->SetYTitle("d#phi (rad)");
2794         fhTrackMatchedDPhiTRD[i]->SetXTitle("#it{E}_{cluster} (GeV)");
2795
2796         fhEOverPTRD[i]  = new TH2F
2797         (Form("hEOverPTRD%s",cutTM[i].Data()),
2798          Form("matched track E/p vs cluster E, behind TRD, %s",cutTM[i].Data()),
2799          nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
2800         fhEOverPTRD[i]->SetXTitle("#it{E} (GeV)");
2801         fhEOverPTRD[i]->SetYTitle("E/p");
2802
2803         outputContainer->Add(fhTrackMatchedDEtaTRD[i]) ;
2804         outputContainer->Add(fhTrackMatchedDPhiTRD[i]) ;
2805         outputContainer->Add(fhEOverPTRD[i]);
2806       }
2807
2808       if(IsDataMC())
2809       {
2810         fhTrackMatchedDEtaMCNoOverlap[i]  = new TH2F
2811         (Form("hTrackMatchedDEtaMCNoOverlap%s",cutTM[i].Data()),
2812          Form("d#eta of cluster-track vs cluster energy, no other MC particles overlap %s",cutTM[i].Data()),
2813          nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
2814         fhTrackMatchedDEtaMCNoOverlap[i]->SetYTitle("d#eta");
2815         fhTrackMatchedDEtaMCNoOverlap[i]->SetXTitle("#it{E}_{cluster} (GeV)");
2816
2817         fhTrackMatchedDPhiMCNoOverlap[i]  = new TH2F
2818         (Form("hTrackMatchedDPhiMCNoOverlap%s",cutTM[i].Data()),
2819          Form("d#phi of cluster-track vs cluster energy, no other MC particles overlap %s",cutTM[i].Data()),
2820          nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
2821         fhTrackMatchedDPhiMCNoOverlap[i]->SetYTitle("d#phi (rad)");
2822         fhTrackMatchedDPhiMCNoOverlap[i]->SetXTitle("#it{E}_{cluster} (GeV)");
2823
2824         outputContainer->Add(fhTrackMatchedDEtaMCNoOverlap[i]) ;
2825         outputContainer->Add(fhTrackMatchedDPhiMCNoOverlap[i]) ;
2826         fhTrackMatchedDEtaMCOverlap[i]  = new TH2F
2827         (Form("hTrackMatchedDEtaMCOverlap%s",cutTM[i].Data()),
2828          Form("d#eta of cluster-track vs cluster energy, several MC particles overlap %s",cutTM[i].Data()),
2829          nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
2830         fhTrackMatchedDEtaMCOverlap[i]->SetYTitle("d#eta");
2831         fhTrackMatchedDEtaMCOverlap[i]->SetXTitle("#it{E}_{cluster} (GeV)");
2832
2833         fhTrackMatchedDPhiMCOverlap[i]  = new TH2F
2834         (Form("hTrackMatchedDPhiMCOverlap%s",cutTM[i].Data()),
2835          Form("d#phi of cluster-track vs cluster energy, several MC particles overlap %s",cutTM[i].Data()),
2836          nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
2837         fhTrackMatchedDPhiMCOverlap[i]->SetYTitle("d#phi (rad)");
2838         fhTrackMatchedDPhiMCOverlap[i]->SetXTitle("#it{E}_{cluster} (GeV)");
2839
2840         outputContainer->Add(fhTrackMatchedDEtaMCOverlap[i]) ;
2841         outputContainer->Add(fhTrackMatchedDPhiMCOverlap[i]) ;
2842
2843         fhTrackMatchedDEtaMCConversion[i]  = new TH2F
2844         (Form("hTrackMatchedDEtaMCConversion%s",cutTM[i].Data()),
2845          Form("d#eta of cluster-track vs cluster energy, no other MC particles overlap appart from conversions %s",cutTM[i].Data()),
2846          nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
2847         fhTrackMatchedDEtaMCConversion[i]->SetYTitle("d#eta");
2848         fhTrackMatchedDEtaMCConversion[i]->SetXTitle("#it{E}_{cluster} (GeV)");
2849
2850         fhTrackMatchedDPhiMCConversion[i]  = new TH2F
2851         (Form("hTrackMatchedDPhiMCConversion%s",cutTM[i].Data()),
2852          Form("d#phi of cluster-track vs cluster energy, no other MC particles overlap appart from conversions %s",cutTM[i].Data()),
2853          nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
2854         fhTrackMatchedDPhiMCConversion[i]->SetYTitle("d#phi (rad)");
2855         fhTrackMatchedDPhiMCConversion[i]->SetXTitle("#it{E}_{cluster} (GeV)");
2856
2857         outputContainer->Add(fhTrackMatchedDEtaMCConversion[i]) ;
2858         outputContainer->Add(fhTrackMatchedDPhiMCConversion[i]) ;
2859
2860         fhTrackMatchedMCParticle[i]  = new TH2F
2861         (Form("hTrackMatchedMCParticle%s",cutTM[i].Data()),
2862          Form("Origin of particle vs energy %s",cutTM[i].Data()),
2863          nptbins,ptmin,ptmax,8,0,8);
2864         fhTrackMatchedMCParticle[i]->SetXTitle("#it{E} (GeV)");
2865         //fhTrackMatchedMCParticle[i]->SetYTitle("Particle type");
2866
2867         fhTrackMatchedMCParticle[i]->GetYaxis()->SetBinLabel(1 ,"Photon");
2868         fhTrackMatchedMCParticle[i]->GetYaxis()->SetBinLabel(2 ,"Electron");
2869         fhTrackMatchedMCParticle[i]->GetYaxis()->SetBinLabel(3 ,"Meson Merged");
2870         fhTrackMatchedMCParticle[i]->GetYaxis()->SetBinLabel(4 ,"Rest");
2871         fhTrackMatchedMCParticle[i]->GetYaxis()->SetBinLabel(5 ,"Conv. Photon");
2872         fhTrackMatchedMCParticle[i]->GetYaxis()->SetBinLabel(6 ,"Conv. Electron");
2873         fhTrackMatchedMCParticle[i]->GetYaxis()->SetBinLabel(7 ,"Conv. Merged");
2874         fhTrackMatchedMCParticle[i]->GetYaxis()->SetBinLabel(8 ,"Conv. Rest");
2875
2876         outputContainer->Add(fhTrackMatchedMCParticle[i]);
2877       }
2878     }
2879   }
2880
2881   if(fFillPileUpHistograms)
2882   {
2883
2884     TString pileUpName[] = {"SPD","EMCAL","SPDOrEMCAL","SPDAndEMCAL","SPDAndNotEMCAL","EMCALAndNotSPD","NotSPDAndNotEMCAL"} ;
2885
2886     for(Int_t i = 0 ; i < 7 ; i++)
2887     {
2888       fhPtPileUp[i]  = new TH1F(Form("hPtPileUp%s",pileUpName[i].Data()),
2889                                 Form("Cluster  #it{p}_[T] distribution, %s Pile-Up event",pileUpName[i].Data()), nptbins,ptmin,ptmax);
2890       fhPtPileUp[i]->SetXTitle("#it{p}_[T] (GeV/#it{c})");
2891       outputContainer->Add(fhPtPileUp[i]);
2892
2893       fhPtChargedPileUp[i]  = new TH1F(Form("hPtChargedPileUp%s",pileUpName[i].Data()),
2894                                        Form("Charged clusters #it{p}_[T] distribution, %s Pile-Up event",pileUpName[i].Data()), nptbins,ptmin,ptmax);
2895       fhPtChargedPileUp[i]->SetXTitle("#it{p}_[T] (GeV/#it{c})");
2896       outputContainer->Add(fhPtChargedPileUp[i]);
2897
2898       fhPtPhotonPileUp[i]  = new TH1F(Form("hPtPhotonPileUp%s",pileUpName[i].Data()),
2899                                       Form("Selected photon #it{p}_[T] distribution, %s Pile-Up event",pileUpName[i].Data()), nptbins,ptmin,ptmax);
2900       fhPtPhotonPileUp[i]->SetXTitle("#it{p}_[T] (GeV/#it{c})");
2901       outputContainer->Add(fhPtPhotonPileUp[i]);
2902
2903
2904       fhClusterEFracLongTimePileUp[i]  = new TH2F(Form("hClusterEFracLongTimePileUp%s",pileUpName[i].Data()),
2905                                                   Form("Cluster E vs fraction of cluster energy from large T cells, %s Pile-Up event",pileUpName[i].Data()),
2906                                                   nptbins,ptmin,ptmax,200,0,1);
2907       fhClusterEFracLongTimePileUp[i]->SetXTitle("#it{E} (GeV)");
2908       fhClusterEFracLongTimePileUp[i]->SetYTitle("E(large time) / E");
2909       outputContainer->Add(fhClusterEFracLongTimePileUp[i]);
2910
2911       fhClusterCellTimePileUp[i]  = new TH2F(Form("hClusterCellTimePileUp%s",pileUpName[i].Data()),
2912                                              Form("Cluster E vs cell time in cluster, %s Pile-Up event",pileUpName[i].Data()),
2913                                              nptbins,ptmin,ptmax,ntimebins,timemin,timemax);
2914       fhClusterCellTimePileUp[i]->SetXTitle("#it{E} (GeV)");
2915       fhClusterCellTimePileUp[i]->SetYTitle("t_{cell} (ns)");
2916       outputContainer->Add(fhClusterCellTimePileUp[i]);
2917
2918       fhClusterTimeDiffPileUp[i]  = new TH2F(Form("hClusterTimeDiffPileUp%s",pileUpName[i].Data()),
2919                                              Form("Cluster E vs t_{max}-t_{cell} in cluster, %s Pile-Up event",pileUpName[i].Data()),
2920                                              nptbins,ptmin,ptmax,400,-200,200);
2921       fhClusterTimeDiffPileUp[i]->SetXTitle("#it{E} (GeV)");
2922       fhClusterTimeDiffPileUp[i]->SetYTitle("t_{max}-t_{cell} (ns)");
2923       outputContainer->Add(fhClusterTimeDiffPileUp[i]);
2924
2925       fhClusterTimeDiffChargedPileUp[i]  = new TH2F(Form("hClusterTimeDiffChargedPileUp%s",pileUpName[i].Data()),
2926                                                     Form("Charged clusters E vs t_{max}-t_{cell} in cluster, %s Pile-Up event",pileUpName[i].Data()),
2927                                                     nptbins,ptmin,ptmax,400,-200,200);
2928       fhClusterTimeDiffChargedPileUp[i]->SetXTitle("#it{E} (GeV)");
2929       fhClusterTimeDiffChargedPileUp[i]->SetYTitle("t_{max}-t_{cell} (ns)");
2930       outputContainer->Add(fhClusterTimeDiffChargedPileUp[i]);
2931
2932       fhClusterTimeDiffPhotonPileUp[i]  = new TH2F(Form("hClusterTimeDiffPhotonPileUp%s",pileUpName[i].Data()),
2933                                                    Form("Selected photon E vs t_{max}-t_{cell} in cluster, %s Pile-Up event",pileUpName[i].Data()),
2934                                                    nptbins,ptmin,ptmax,400,-200,200);
2935       fhClusterTimeDiffPhotonPileUp[i]->SetXTitle("#it{E} (GeV)");
2936       fhClusterTimeDiffPhotonPileUp[i]->SetYTitle("t_{max}-t_{cell} (ns)");
2937       outputContainer->Add(fhClusterTimeDiffPhotonPileUp[i]);
2938
2939       fhLambda0PileUp[i]  = new TH2F(Form("hLambda0PileUp%s",pileUpName[i].Data()),
2940                                      Form("Cluster E vs #lambda^{2}_{0} in cluster, %s Pile-Up event",pileUpName[i].Data()),
2941                                      nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2942       fhLambda0PileUp[i]->SetXTitle("#it{E} (GeV)");
2943       fhLambda0PileUp[i]->SetYTitle("#lambda^{2}_{0}");
2944       outputContainer->Add(fhLambda0PileUp[i]);
2945
2946       fhLambda0ChargedPileUp[i]  = new TH2F(Form("hLambda0ChargedPileUp%s",pileUpName[i].Data()),
2947                                             Form("Charged clusters E vs #lambda^{2}_{0}in cluster, %s Pile-Up event",pileUpName[i].Data()), nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2948       fhLambda0ChargedPileUp[i]->SetXTitle("#it{E} (GeV)");
2949       fhLambda0ChargedPileUp[i]->SetYTitle("#lambda^{2}_{0}");
2950       outputContainer->Add(fhLambda0ChargedPileUp[i]);
2951
2952     }
2953
2954     fhEtaPhiBC0  = new TH2F ("hEtaPhiBC0","eta-phi for clusters tof corresponding to BC=0",netabins,etamin,etamax, nphibins,phimin,phimax);
2955     fhEtaPhiBC0->SetXTitle("#eta ");
2956     fhEtaPhiBC0->SetYTitle("#phi (rad)");
2957     outputContainer->Add(fhEtaPhiBC0);
2958
2959     fhEtaPhiBCPlus  = new TH2F ("hEtaPhiBCPlus","eta-phi for clusters tof corresponding to BC>0",netabins,etamin,etamax, nphibins,phimin,phimax);
2960     fhEtaPhiBCPlus->SetXTitle("#eta ");
2961     fhEtaPhiBCPlus->SetYTitle("#phi (rad)");
2962     outputContainer->Add(fhEtaPhiBCPlus);
2963
2964     fhEtaPhiBCMinus  = new TH2F ("hEtaPhiBCMinus","eta-phi for clusters tof corresponding to BC<0",netabins,etamin,etamax, nphibins,phimin,phimax);
2965     fhEtaPhiBCMinus->SetXTitle("#eta ");
2966     fhEtaPhiBCMinus->SetYTitle("#phi (rad)");
2967     outputContainer->Add(fhEtaPhiBCMinus);
2968
2969     fhEtaPhiBC0PileUpSPD  = new TH2F ("hEtaPhiBC0PileUpSPD","eta-phi for clusters tof corresponding to BC=0, SPD pile-up",netabins,etamin,etamax, nphibins,phimin,phimax);
2970     fhEtaPhiBC0PileUpSPD->SetXTitle("#eta ");
2971     fhEtaPhiBC0PileUpSPD->SetYTitle("#phi (rad)");
2972     outputContainer->Add(fhEtaPhiBC0PileUpSPD);
2973
2974     fhEtaPhiBCPlusPileUpSPD  = new TH2F ("hEtaPhiBCPlusPileUpSPD","eta-phi for clusters tof corresponding to BC>0, SPD pile-up",netabins,etamin,etamax, nphibins,phimin,phimax);
2975     fhEtaPhiBCPlusPileUpSPD->SetXTitle("#eta ");
2976     fhEtaPhiBCPlusPileUpSPD->SetYTitle("#phi (rad)");
2977     outputContainer->Add(fhEtaPhiBCPlusPileUpSPD);
2978
2979     fhEtaPhiBCMinusPileUpSPD  = new TH2F ("hEtaPhiBCMinusPileUpSPD","eta-phi for clusters tof corresponding to BC<0, SPD pile-up",netabins,etamin,etamax, nphibins,phimin,phimax);
2980     fhEtaPhiBCMinusPileUpSPD->SetXTitle("#eta ");
2981     fhEtaPhiBCMinusPileUpSPD->SetYTitle("#phi (rad)");
2982     outputContainer->Add(fhEtaPhiBCMinusPileUpSPD);
2983
2984     fhTimePtNoCut  = new TH2F ("hTimePt_NoCut","time of cluster vs pT of clusters, no cut", nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
2985     fhTimePtNoCut->SetXTitle("#it{p}_[T] (GeV/#it{c})");
2986     fhTimePtNoCut->SetYTitle("#it{time} (ns)");
2987     outputContainer->Add(fhTimePtNoCut);
2988
2989     fhTimePtSPD  = new TH2F ("hTimePt_SPD","time of cluster vs pT of clusters, SPD cut", nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
2990     fhTimePtSPD->SetXTitle("#it{p}_[T] (GeV/#it{c})");
2991     fhTimePtSPD->SetYTitle("#it{time} (ns)");
2992     outputContainer->Add(fhTimePtSPD);
2993
2994     fhTimePtPhotonNoCut  = new TH2F ("hTimePtPhoton_NoCut","time of photon cluster vs pT of clusters, no cut", nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
2995     fhTimePtPhotonNoCut->SetXTitle("#it{p}_[T] (GeV/#it{c})");
2996     fhTimePtPhotonNoCut->SetYTitle("#it{time} (ns)");
2997     outputContainer->Add(fhTimePtPhotonNoCut);
2998
2999     fhTimePtPhotonSPD  = new TH2F ("hTimePtPhoton_SPD","time of  photon cluster vs pT of clusters, SPD cut", nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
3000     fhTimePtPhotonSPD->SetXTitle("#it{p}_[T] (GeV/#it{c})");
3001     fhTimePtPhotonSPD->SetYTitle("#it{time} (ns)");
3002     outputContainer->Add(fhTimePtPhotonSPD);
3003
3004     fhTimeNPileUpVertSPD  = new TH2F ("hTime_NPileUpVertSPD","time of cluster vs N pile-up SPD vertex", ntimebins,timemin,timemax,20,0,20);
3005     fhTimeNPileUpVertSPD->SetYTitle("# vertex ");
3006     fhTimeNPileUpVertSPD->SetXTitle("#it{time} (ns)");
3007     outputContainer->Add(fhTimeNPileUpVertSPD);
3008
3009     fhTimeNPileUpVertTrack  = new TH2F ("hTime_NPileUpVertTracks","time of cluster vs N pile-up Tracks vertex", ntimebins,timemin,timemax, 20,0,20 );
3010     fhTimeNPileUpVertTrack->SetYTitle("# vertex ");
3011     fhTimeNPileUpVertTrack->SetXTitle("#it{time} (ns)");
3012     outputContainer->Add(fhTimeNPileUpVertTrack);
3013
3014     fhTimeNPileUpVertContributors  = new TH2F ("hTime_NPileUpVertContributors","time of cluster vs N constributors to pile-up SPD vertex", ntimebins,timemin,timemax,50,0,50);
3015     fhTimeNPileUpVertContributors->SetYTitle("# vertex ");
3016     fhTimeNPileUpVertContributors->SetXTitle("#it{time} (ns)");
3017     outputContainer->Add(fhTimeNPileUpVertContributors);
3018
3019     fhTimePileUpMainVertexZDistance  = new TH2F ("hTime_PileUpMainVertexZDistance","time of cluster vs distance in Z pile-up SPD vertex - main SPD vertex",ntimebins,timemin,timemax,100,0,50);
3020     fhTimePileUpMainVertexZDistance->SetYTitle("distance Z (cm) ");
3021     fhTimePileUpMainVertexZDistance->SetXTitle("#it{time} (ns)");
3022     outputContainer->Add(fhTimePileUpMainVertexZDistance);
3023
3024     fhTimePileUpMainVertexZDiamond  = new TH2F ("hTime_PileUpMainVertexZDiamond","time of cluster vs distance in Z pile-up SPD vertex - z diamond",ntimebins,timemin,timemax,100,0,50);
3025     fhTimePileUpMainVertexZDiamond->SetYTitle("diamond distance Z (cm) ");
3026     fhTimePileUpMainVertexZDiamond->SetXTitle("#it{time} (ns)");
3027     outputContainer->Add(fhTimePileUpMainVertexZDiamond);
3028
3029     TString title[] = {"no |t diff| cut","|t diff|<20 ns","|t diff|>20 ns","|t diff|>40 ns"};
3030     TString name [] = {"TDiffNoCut","TDiffSmaller20ns","TDiffLarger20ns","TDiffLarger40ns"};
3031     for(Int_t i = 0; i < 4; i++)
3032     {
3033       fhClusterMultSPDPileUp[i] = new TH2F(Form("fhClusterMultSPDPileUp_%s", name[i].Data()),
3034                                            Form("Number of clusters per pile up event with #it{E} > 0.5 and %s respect cluster max vs cluster max E ",title[i].Data()),
3035                                            nptbins,ptmin,ptmax,100,0,100);
3036       fhClusterMultSPDPileUp[i]->SetYTitle("n clusters ");
3037       fhClusterMultSPDPileUp[i]->SetXTitle("#it{E}_{cluster max} (GeV)");
3038       outputContainer->Add(fhClusterMultSPDPileUp[i]) ;
3039
3040       fhClusterMultNoPileUp[i] = new TH2F(Form("fhClusterMultNoPileUp_%s", name[i].Data()),
3041                                           Form("Number of clusters per non pile up event with #it{E} > 0.5 and %s respect cluster max vs cluster max E ",title[i].Data()),
3042                                           nptbins,ptmin,ptmax,100,0,100);
3043       fhClusterMultNoPileUp[i]->SetYTitle("n clusters ");
3044       fhClusterMultNoPileUp[i]->SetXTitle("#it{E}_{cluster max} (GeV)");
3045       outputContainer->Add(fhClusterMultNoPileUp[i]) ;
3046     }
3047
3048     fhPtNPileUpSPDVtx  = new TH2F ("hPt_NPileUpVertSPD","pT of cluster vs N pile-up SPD vertex",
3049                                    nptbins,ptmin,ptmax,20,0,20);
3050     fhPtNPileUpSPDVtx->SetYTitle("# vertex ");
3051     fhPtNPileUpSPDVtx->SetXTitle("#it{p}_[T] (GeV/#it{c})");
3052     outputContainer->Add(fhPtNPileUpSPDVtx);
3053
3054     fhPtNPileUpTrkVtx  = new TH2F ("hPt_NPileUpVertTracks","pT of cluster vs N pile-up Tracks vertex",
3055                                    nptbins,ptmin,ptmax, 20,0,20 );
3056     fhPtNPileUpTrkVtx->SetYTitle("# vertex ");
3057     fhPtNPileUpTrkVtx->SetXTitle("#it{p}_[T] (GeV/#it{c})");
3058     outputContainer->Add(fhPtNPileUpTrkVtx);
3059
3060     fhPtNPileUpSPDVtxTimeCut  = new TH2F ("hPt_NPileUpVertSPD_TimeCut","pT of cluster vs N pile-up SPD vertex, |tof| < 25 ns",
3061                                           nptbins,ptmin,ptmax,20,0,20);
3062     fhPtNPileUpSPDVtxTimeCut->SetYTitle("# vertex ");
3063     fhPtNPileUpSPDVtxTimeCut->SetXTitle("#it{p}_[T] (GeV/#it{c})");
3064     outputContainer->Add(fhPtNPileUpSPDVtxTimeCut);
3065
3066     fhPtNPileUpTrkVtxTimeCut  = new TH2F ("hPt_NPileUpVertTracks_TimeCut","pT of cluster vs N pile-up Tracks vertex, |tof| < 25 ns",
3067                                           nptbins,ptmin,ptmax, 20,0,20 );
3068     fhPtNPileUpTrkVtxTimeCut->SetYTitle("# vertex ");
3069     fhPtNPileUpTrkVtxTimeCut->SetXTitle("#it{p}_[T] (GeV/#it{c})");
3070     outputContainer->Add(fhPtNPileUpTrkVtxTimeCut);
3071
3072                 fhPtNPileUpSPDVtxTimeCut2  = new TH2F ("hPt_NPileUpVertSPD_TimeCut2","pT of cluster vs N pile-up SPD vertex, -25 < tof < 75 ns",
3073                                            nptbins,ptmin,ptmax,20,0,20);
3074     fhPtNPileUpSPDVtxTimeCut2->SetYTitle("# vertex ");
3075     fhPtNPileUpSPDVtxTimeCut2->SetXTitle("#it{p}_[T] (GeV/#it{c})");
3076     outputContainer->Add(fhPtNPileUpSPDVtxTimeCut2);
3077
3078     fhPtNPileUpTrkVtxTimeCut2  = new TH2F ("hPt_NPileUpVertTracks_TimeCut2","pT of cluster vs N pile-up Tracks vertex, -25 < tof < 75 ns",
3079                                            nptbins,ptmin,ptmax, 20,0,20 );
3080     fhPtNPileUpTrkVtxTimeCut2->SetYTitle("# vertex ");
3081     fhPtNPileUpTrkVtxTimeCut2->SetXTitle("#it{p}_[T] (GeV/#it{c})");
3082     outputContainer->Add(fhPtNPileUpTrkVtxTimeCut2);
3083
3084     fhPtPhotonNPileUpSPDVtx  = new TH2F ("hPtPhoton_NPileUpVertSPD","pT of cluster vs N pile-up SPD vertex",
3085                                          nptbins,ptmin,ptmax,20,0,20);
3086     fhPtPhotonNPileUpSPDVtx->SetYTitle("# vertex ");
3087     fhPtPhotonNPileUpSPDVtx->SetXTitle("#it{p}_[T] (GeV/#it{c})");
3088     outputContainer->Add(fhPtPhotonNPileUpSPDVtx);
3089
3090     fhPtPhotonNPileUpTrkVtx  = new TH2F ("hPtPhoton_NPileUpVertTracks","pT of cluster vs N pile-up Tracks vertex",
3091                                          nptbins,ptmin,ptmax, 20,0,20 );
3092     fhPtPhotonNPileUpTrkVtx->SetYTitle("# vertex ");
3093     fhPtPhotonNPileUpTrkVtx->SetXTitle("#it{p}_[T] (GeV/#it{c})");
3094     outputContainer->Add(fhPtPhotonNPileUpTrkVtx);
3095
3096     fhPtPhotonNPileUpSPDVtxTimeCut  = new TH2F ("hPtPhoton_NPileUpVertSPD_TimeCut","pT of cluster vs N pile-up SPD vertex, |tof| < 25 ns",
3097                                                 nptbins,ptmin,ptmax,20,0,20);
3098     fhPtPhotonNPileUpSPDVtxTimeCut->SetYTitle("# vertex ");
3099     fhPtPhotonNPileUpSPDVtxTimeCut->SetXTitle("#it{p}_[T] (GeV/#it{c})");
3100     outputContainer->Add(fhPtPhotonNPileUpSPDVtxTimeCut);
3101
3102     fhPtPhotonNPileUpTrkVtxTimeCut  = new TH2F ("hPtPhoton_NPileUpVertTracks_TimeCut","pT of cluster vs N pile-up Tracks vertex, |tof| < 25 ns",
3103                                                 nptbins,ptmin,ptmax, 20,0,20 );
3104     fhPtPhotonNPileUpTrkVtxTimeCut->SetYTitle("# vertex ");
3105     fhPtPhotonNPileUpTrkVtxTimeCut->SetXTitle("#it{p}_[T] (GeV/#it{c})");
3106     outputContainer->Add(fhPtPhotonNPileUpTrkVtxTimeCut);
3107
3108     fhPtPhotonNPileUpSPDVtxTimeCut2  = new TH2F ("hPtPhoton_NPileUpVertSPD_TimeCut2","pT of cluster vs N pile-up SPD vertex, -25 < tof < 75 ns",
3109                                                  nptbins,ptmin,ptmax,20,0,20);
3110     fhPtPhotonNPileUpSPDVtxTimeCut2->SetYTitle("# vertex ");
3111     fhPtPhotonNPileUpSPDVtxTimeCut2->SetXTitle("#it{p}_[T] (GeV/#it{c})");
3112     outputContainer->Add(fhPtPhotonNPileUpSPDVtxTimeCut2);
3113
3114     fhPtPhotonNPileUpTrkVtxTimeCut2  = new TH2F ("hPtPhoton_NPileUpVertTracks_TimeCut2","pT of cluster vs N pile-up Tracks vertex,  -25 < tof < 75 ns",
3115                                                  nptbins,ptmin,ptmax, 20,0,20 );
3116     fhPtPhotonNPileUpTrkVtxTimeCut2->SetYTitle("# vertex ");
3117     fhPtPhotonNPileUpTrkVtxTimeCut2->SetXTitle("#it{p}_[T] (GeV/#it{c})");
3118     outputContainer->Add(fhPtPhotonNPileUpTrkVtxTimeCut2);
3119
3120   }
3121
3122   if(IsDataMC())
3123   {
3124     TString ptype[] = { "#gamma", "#gamma_{#pi decay}","#gamma_{other decay}", "#pi^{0}","#eta",
3125       "e^{#pm}","#gamma->e^{#pm}","hadron?","Anti-N","Anti-P",
3126       "#gamma_{prompt}","#gamma_{fragmentation}","#gamma_{ISR}","String"                                } ;
3127
3128     TString pname[] = { "Photon","PhotonPi0Decay","PhotonOtherDecay","Pi0","Eta","Electron",
3129       "Conversion", "Hadron", "AntiNeutron","AntiProton",
3130       "PhotonPrompt","PhotonFragmentation","PhotonISR","String" } ;
3131
3132     for(Int_t i = 0; i < fNOriginHistograms; i++)
3133     {
3134       fhMCE[i]  = new TH1F(Form("hE_MC%s",pname[i].Data()),
3135                            Form("cluster from %s : E ",ptype[i].Data()),
3136                            nptbins,ptmin,ptmax);
3137       fhMCE[i]->SetXTitle("#it{E} (GeV)");
3138       outputContainer->Add(fhMCE[i]) ;
3139
3140       fhMCPt[i]  = new TH1F(Form("hPt_MC%s",pname[i].Data()),
3141                             Form("cluster from %s : #it{p}_[T] ",ptype[i].Data()),
3142                             nptbins,ptmin,ptmax);
3143       fhMCPt[i]->SetXTitle("#it{p}_[T] (GeV/#it{c})");
3144       outputContainer->Add(fhMCPt[i]) ;
3145
3146       fhMCEta[i]  = new TH2F(Form("hEta_MC%s",pname[i].Data()),
3147                              Form("cluster from %s : #eta ",ptype[i].Data()),
3148                              nptbins,ptmin,ptmax,netabins,etamin,etamax);
3149       fhMCEta[i]->SetYTitle("#eta");
3150       fhMCEta[i]->SetXTitle("#it{E} (GeV)");
3151       outputContainer->Add(fhMCEta[i]) ;
3152
3153       fhMCPhi[i]  = new TH2F(Form("hPhi_MC%s",pname[i].Data()),
3154                              Form("cluster from %s : #phi ",ptype[i].Data()),
3155                              nptbins,ptmin,ptmax,nphibins,phimin,phimax);
3156       fhMCPhi[i]->SetYTitle("#phi (rad)");
3157       fhMCPhi[i]->SetXTitle("#it{E} (GeV)");
3158       outputContainer->Add(fhMCPhi[i]) ;
3159
3160
3161       fhMCDeltaE[i]  = new TH2F (Form("hDeltaE_MC%s",pname[i].Data()),
3162                                  Form("MC - Reco E from %s",pname[i].Data()),
3163                                  nptbins,ptmin,ptmax, 200,-50,50);
3164       fhMCDeltaE[i]->SetYTitle("#Delta #it{E} (GeV)");
3165       fhMCDeltaE[i]->SetXTitle("#it{E} (GeV)");
3166       outputContainer->Add(fhMCDeltaE[i]);
3167
3168       fhMCDeltaPt[i]  = new TH2F (Form("hDeltaPt_MC%s",pname[i].Data()),
3169                                   Form("MC - Reco #it{p}_[T] from %s",pname[i].Data()),
3170                                   nptbins,ptmin,ptmax, 200,-50,50);
3171       fhMCDeltaPt[i]->SetXTitle("p_{T,rec} (GeV/#it{c})");
3172       fhMCDeltaPt[i]->SetYTitle("#Delta #it{p}_[T] (GeV/#it{c})");
3173       outputContainer->Add(fhMCDeltaPt[i]);
3174
3175       fhMC2E[i]  = new TH2F (Form("h2E_MC%s",pname[i].Data()),
3176                              Form("E distribution, reconstructed vs generated from %s",pname[i].Data()),
3177                              nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
3178       fhMC2E[i]->SetXTitle("#it{E}_{rec} (GeV)");
3179       fhMC2E[i]->SetYTitle("#it{E}_{gen} (GeV)");
3180       outputContainer->Add(fhMC2E[i]);
3181
3182       fhMC2Pt[i]  = new TH2F (Form("h2Pt_MC%s",pname[i].Data()),
3183                               Form("p_T distribution, reconstructed vs generated from %s",pname[i].Data()),
3184                               nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
3185       fhMC2Pt[i]->SetXTitle("p_{T,rec} (GeV/#it{c})");
3186       fhMC2Pt[i]->SetYTitle("p_{T,gen} (GeV/#it{c})");
3187       outputContainer->Add(fhMC2Pt[i]);
3188
3189
3190     }
3191
3192     TString pptype[] = { "#gamma", "#gamma_{#pi decay}","#gamma_{other decay}","hadron?",
3193       "#gamma_{prompt}","#gamma_{fragmentation}","#gamma_{ISR}"} ;
3194
3195     TString ppname[] = { "Photon","PhotonPi0Decay","PhotonOtherDecay","Hadron",
3196       "PhotonPrompt","PhotonFragmentation","PhotonISR"} ;
3197
3198     for(Int_t i = 0; i < fNPrimaryHistograms; i++)
3199     {
3200       fhEPrimMC[i]  = new TH1F(Form("hEPrim_MC%s",ppname[i].Data()),
3201                                Form("primary photon %s : E ",pptype[i].Data()),
3202                                nptbins,ptmin,ptmax);
3203       fhEPrimMC[i]->SetXTitle("#it{E} (GeV)");
3204       outputContainer->Add(fhEPrimMC[i]) ;
3205
3206       fhPtPrimMC[i]  = new TH1F(Form("hPtPrim_MC%s",ppname[i].Data()),
3207                                 Form("primary photon %s : #it{p}_[T] ",pptype[i].Data()),
3208                                 nptbins,ptmin,ptmax);
3209       fhPtPrimMC[i]->SetXTitle("#it{p}_[T] (GeV/#it{c})");
3210       outputContainer->Add(fhPtPrimMC[i]) ;
3211
3212       fhYPrimMC[i]  = new TH2F(Form("hYPrim_MC%s",ppname[i].Data()),
3213                                Form("primary photon %s : Rapidity ",pptype[i].Data()),
3214                                nptbins,ptmin,ptmax,200,-2,2);
3215       fhYPrimMC[i]->SetYTitle("Rapidity");
3216       fhYPrimMC[i]->SetXTitle("#it{E} (GeV)");
3217       outputContainer->Add(fhYPrimMC[i]) ;
3218
3219       fhEtaPrimMC[i]  = new TH2F(Form("hEtaPrim_MC%s",ppname[i].Data()),
3220                                Form("primary photon %s : #eta",pptype[i].Data()),
3221                                nptbins,ptmin,ptmax,200,-2,2);
3222       fhEtaPrimMC[i]->SetYTitle("#eta");
3223       fhEtaPrimMC[i]->SetXTitle("#it{E} (GeV)");
3224       outputContainer->Add(fhEtaPrimMC[i]) ;
3225
3226       fhPhiPrimMC[i]  = new TH2F(Form("hPhiPrim_MC%s",ppname[i].Data()),
3227                                  Form("primary photon %s : #phi ",pptype[i].Data()),
3228                                  nptbins,ptmin,ptmax,nphibins,0,TMath::TwoPi());
3229       fhPhiPrimMC[i]->SetYTitle("#phi (rad)");
3230       fhPhiPrimMC[i]->SetXTitle("#it{E} (GeV)");
3231       outputContainer->Add(fhPhiPrimMC[i]) ;
3232
3233
3234       fhEPrimMCAcc[i]  = new TH1F(Form("hEPrimAcc_MC%s",ppname[i].Data()),
3235                                   Form("primary photon %s in acceptance: E ",pptype[i].Data()),
3236                                   nptbins,ptmin,ptmax);
3237       fhEPrimMCAcc[i]->SetXTitle("#it{E} (GeV)");
3238       outputContainer->Add(fhEPrimMCAcc[i]) ;
3239
3240       fhPtPrimMCAcc[i]  = new TH1F(Form("hPtPrimAcc_MC%s",ppname[i].Data()),
3241                                    Form("primary photon %s in acceptance: #it{p}_[T] ",pptype[i].Data()),
3242                                    nptbins,ptmin,ptmax);
3243       fhPtPrimMCAcc[i]->SetXTitle("#it{p}_[T] (GeV/#it{c})");
3244       outputContainer->Add(fhPtPrimMCAcc[i]) ;
3245
3246       fhYPrimMCAcc[i]  = new TH2F(Form("hYPrimAcc_MC%s",ppname[i].Data()),
3247                                   Form("primary photon %s in acceptance: Rapidity ",pptype[i].Data()),
3248                                   nptbins,ptmin,ptmax,100,-1,1);
3249       fhYPrimMCAcc[i]->SetYTitle("Rapidity");
3250       fhYPrimMCAcc[i]->SetXTitle("#it{E} (GeV)");
3251       outputContainer->Add(fhYPrimMCAcc[i]) ;
3252
3253       fhEtaPrimMCAcc[i]  = new TH2F(Form("hEtaPrimAcc_MC%s",ppname[i].Data()),
3254                                   Form("primary photon %s in acceptance: #eta ",pptype[i].Data()),
3255                                   nptbins,ptmin,ptmax,netabins,etamin,etamax);
3256       fhEtaPrimMCAcc[i]->SetYTitle("#eta");
3257       fhEtaPrimMCAcc[i]->SetXTitle("#it{E} (GeV)");
3258       outputContainer->Add(fhEtaPrimMCAcc[i]) ;
3259
3260       fhPhiPrimMCAcc[i]  = new TH2F(Form("hPhiPrimAcc_MC%s",ppname[i].Data()),
3261                                     Form("primary photon %s in acceptance: #phi ",pptype[i].Data()),
3262                                     nptbins,ptmin,ptmax,nphibins,phimin,phimax);
3263       fhPhiPrimMCAcc[i]->SetYTitle("#phi (rad)");
3264       fhPhiPrimMCAcc[i]->SetXTitle("#it{E} (GeV)");
3265       outputContainer->Add(fhPhiPrimMCAcc[i]) ;
3266
3267     }
3268
3269     if(fFillSSHistograms)
3270     {
3271       TString ptypess[] = { "#gamma","hadron?","#pi^{0}","#eta","#gamma->e^{#pm}","e^{#pm}"} ;
3272
3273       TString pnamess[] = { "Photon","Hadron","Pi0","Eta","Conversion","Electron"} ;
3274
3275       for(Int_t i = 0; i < 6; i++)
3276       {
3277         fhMCELambda0[i]  = new TH2F(Form("hELambda0_MC%s",pnamess[i].Data()),
3278                                     Form("cluster from %s : E vs #lambda_{0}^{2}",ptypess[i].Data()),
3279                                     nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3280         fhMCELambda0[i]->SetYTitle("#lambda_{0}^{2}");
3281         fhMCELambda0[i]->SetXTitle("#it{E} (GeV)");
3282         outputContainer->Add(fhMCELambda0[i]) ;
3283
3284         fhMCELambda1[i]  = new TH2F(Form("hELambda1_MC%s",pnamess[i].Data()),
3285                                     Form("cluster from %s : E vs #lambda_{1}^{2}",ptypess[i].Data()),
3286                                     nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3287         fhMCELambda1[i]->SetYTitle("#lambda_{1}^{2}");
3288         fhMCELambda1[i]->SetXTitle("#it{E} (GeV)");
3289         outputContainer->Add(fhMCELambda1[i]) ;
3290
3291         fhMCEDispersion[i]  = new TH2F(Form("hEDispersion_MC%s",pnamess[i].Data()),
3292                                        Form("cluster from %s : E vs dispersion^{2}",ptypess[i].Data()),
3293                                        nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3294         fhMCEDispersion[i]->SetYTitle("D^{2}");
3295         fhMCEDispersion[i]->SetXTitle("#it{E} (GeV)");
3296         outputContainer->Add(fhMCEDispersion[i]) ;
3297
3298         fhMCNCellsE[i]  = new TH2F (Form("hNCellsE_MC%s",pnamess[i].Data()),
3299                                     Form("# of cells in cluster from %s vs E of clusters",ptypess[i].Data()),
3300                                     nptbins,ptmin,ptmax, nbins,nmin,nmax);
3301         fhMCNCellsE[i]->SetXTitle("#it{E} (GeV)");
3302         fhMCNCellsE[i]->SetYTitle("# of cells in cluster");
3303         outputContainer->Add(fhMCNCellsE[i]);
3304
3305         fhMCMaxCellDiffClusterE[i]  = new TH2F (Form("hMaxCellDiffClusterE_MC%s",pnamess[i].Data()),
3306                                                 Form("energy vs difference of cluster energy from %s - max cell energy / cluster energy, good clusters",ptypess[i].Data()),
3307                                                 nptbins,ptmin,ptmax, 500,0,1.);
3308         fhMCMaxCellDiffClusterE[i]->SetXTitle("#it{E}_{cluster} (GeV) ");
3309         fhMCMaxCellDiffClusterE[i]->SetYTitle("(#it{E}_{cluster} - #it{E}_{cell max})/ #it{E}_{cluster}");
3310         outputContainer->Add(fhMCMaxCellDiffClusterE[i]);
3311
3312         if(!fFillOnlySimpleSSHisto)
3313         {
3314           fhMCLambda0vsClusterMaxCellDiffE0[i]  = new TH2F(Form("hLambda0vsClusterMaxCellDiffE0_MC%s",pnamess[i].Data()),
3315                                                            Form("cluster from %s : #lambda^{2}_{0} vs fraction of energy carried by max cell, E < 2 GeV",ptypess[i].Data()),
3316                                                            ssbins,ssmin,ssmax,500,0,1.);
3317           fhMCLambda0vsClusterMaxCellDiffE0[i]->SetXTitle("#lambda_{0}^{2}");
3318           fhMCLambda0vsClusterMaxCellDiffE0[i]->SetYTitle("(#it{E}_{cluster} - #it{E}_{cell max})/ #it{E}_{cluster}");
3319           outputContainer->Add(fhMCLambda0vsClusterMaxCellDiffE0[i]) ;
3320
3321           fhMCLambda0vsClusterMaxCellDiffE2[i]  = new TH2F(Form("hLambda0vsClusterMaxCellDiffE2_MC%s",pnamess[i].Data()),
3322                                                            Form("cluster from %s : #lambda^{2}_{0} vs fraction of energy carried by max cell, 2< E < 6 GeV",ptypess[i].Data()),
3323                                                            ssbins,ssmin,ssmax,500,0,1.);
3324           fhMCLambda0vsClusterMaxCellDiffE2[i]->SetXTitle("#lambda_{0}^{2}");
3325           fhMCLambda0vsClusterMaxCellDiffE2[i]->SetYTitle("(#it{E}_{cluster} - #it{E}_{cell max})/ #it{E}_{cluster}");
3326           outputContainer->Add(fhMCLambda0vsClusterMaxCellDiffE2[i]) ;
3327
3328           fhMCLambda0vsClusterMaxCellDiffE6[i]  = new TH2F(Form("hLambda0vsClusterMaxCellDiffE6_MC%s",pnamess[i].Data()),
3329                                                            Form("cluster from %s : #lambda^{2}_{0} vs fraction of energy carried by max cell, #it{E} > 6 GeV",ptypess[i].Data()),
3330                                                            ssbins,ssmin,ssmax,500,0,1.);
3331           fhMCLambda0vsClusterMaxCellDiffE6[i]->SetXTitle("#lambda_{0}^{2}");
3332           fhMCLambda0vsClusterMaxCellDiffE6[i]->SetYTitle("(#it{E}_{cluster} - #it{E}_{cell max})/ #it{E}_{cluster}");
3333           outputContainer->Add(fhMCLambda0vsClusterMaxCellDiffE6[i]) ;
3334
3335           fhMCNCellsvsClusterMaxCellDiffE0[i]  = new TH2F(Form("hNCellsvsClusterMaxCellDiffE0_MC%s",pnamess[i].Data()),
3336                                                           Form("cluster from %s : N cells in cluster vs fraction of energy carried by max cell, E < 2 GeV",ptypess[i].Data()),
3337                                                           nbins/5,nmin,nmax/5,500,0,1.);
3338           fhMCNCellsvsClusterMaxCellDiffE0[i]->SetXTitle("N cells in cluster");
3339           fhMCNCellsvsClusterMaxCellDiffE0[i]->SetYTitle("(#it{E}_{cluster} - #it{E}_{cell max})/ #it{E}_{cluster}");
3340           outputContainer->Add(fhMCNCellsvsClusterMaxCellDiffE0[i]) ;
3341
3342           fhMCNCellsvsClusterMaxCellDiffE2[i]  = new TH2F(Form("hNCellsvsClusterMaxCellDiffE2_MC%s",pnamess[i].Data()),
3343                                                           Form("cluster from %s : N cells in cluster vs fraction of energy carried by max cell, 2< E < 6 GeV",ptypess[i].Data()),
3344                                                           nbins/5,nmin,nmax/5,500,0,1.);
3345           fhMCNCellsvsClusterMaxCellDiffE2[i]->SetXTitle("N cells in cluster");
3346           fhMCNCellsvsClusterMaxCellDiffE2[i]->SetYTitle("(#it{E}_{cluster} - #it{E}_{cell max})/ #it{E}_{cluster}");
3347           outputContainer->Add(fhMCNCellsvsClusterMaxCellDiffE2[i]) ;
3348
3349           fhMCNCellsvsClusterMaxCellDiffE6[i]  = new TH2F(Form("hNCellsvsClusterMaxCellDiffE6_MC%s",pnamess[i].Data()),
3350                                                           Form("cluster from %s : N cells in cluster vs fraction of energy carried by max cell, #it{E} > 6 GeV",ptypess[i].Data()),
3351                                                           nbins/5,nmin,nmax/5,500,0,1.);
3352           fhMCNCellsvsClusterMaxCellDiffE6[i]->SetXTitle("N cells in cluster");
3353           fhMCNCellsvsClusterMaxCellDiffE6[i]->SetYTitle("#it{E} (GeV)");
3354           outputContainer->Add(fhMCNCellsvsClusterMaxCellDiffE6[i]) ;
3355
3356           if(fCalorimeter=="EMCAL")
3357           {
3358             fhMCEDispEta[i]  = new TH2F (Form("hEDispEtaE_MC%s",pnamess[i].Data()),
3359                                          Form("cluster from %s : #sigma^{2}_{#eta #eta} = #Sigma w_{i}(#eta_{i} - <#eta>)^{2}/ #Sigma w_{i} vs E",ptypess[i].Data()),
3360                                          nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
3361             fhMCEDispEta[i]->SetXTitle("#it{E} (GeV)");
3362             fhMCEDispEta[i]->SetYTitle("#sigma^{2}_{#eta #eta}");
3363             outputContainer->Add(fhMCEDispEta[i]);
3364
3365             fhMCEDispPhi[i]  = new TH2F (Form("hEDispPhiE_MC%s",pnamess[i].Data()),
3366                                          Form("cluster from %s : #sigma^{2}_{#phi #phi} = #Sigma w_{i}(#phi_{i} - <#phi>)^{2} / #Sigma w_{i} vs E",ptypess[i].Data()),
3367                                          nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
3368             fhMCEDispPhi[i]->SetXTitle("#it{E} (GeV)");
3369             fhMCEDispPhi[i]->SetYTitle("#sigma^{2}_{#phi #phi}");
3370             outputContainer->Add(fhMCEDispPhi[i]);
3371
3372             fhMCESumEtaPhi[i]  = new TH2F (Form("hESumEtaPhiE_MC%s",pnamess[i].Data()),
3373                                            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()),
3374                                            nptbins,ptmin,ptmax, 2*ssbins,-ssmax,ssmax);
3375             fhMCESumEtaPhi[i]->SetXTitle("#it{E} (GeV)");
3376             fhMCESumEtaPhi[i]->SetYTitle("#delta^{2}_{#eta #phi}");
3377             outputContainer->Add(fhMCESumEtaPhi[i]);
3378
3379             fhMCEDispEtaPhiDiff[i]  = new TH2F (Form("hEDispEtaPhiDiffE_MC%s",pnamess[i].Data()),
3380                                                 Form("cluster from %s : #sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta} vs E",ptypess[i].Data()),
3381                                                 nptbins,ptmin,ptmax,200,-10,10);
3382             fhMCEDispEtaPhiDiff[i]->SetXTitle("#it{E} (GeV)");
3383             fhMCEDispEtaPhiDiff[i]->SetYTitle("#sigma^{2}_{#phi #phi}-#sigma^{2}_{#eta #eta}");
3384             outputContainer->Add(fhMCEDispEtaPhiDiff[i]);
3385
3386             fhMCESphericity[i]  = new TH2F (Form("hESphericity_MC%s",pnamess[i].Data()),
3387                                             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()),
3388                                             nptbins,ptmin,ptmax, 200,-1,1);
3389             fhMCESphericity[i]->SetXTitle("#it{E} (GeV)");
3390             fhMCESphericity[i]->SetYTitle("s = (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi})");
3391             outputContainer->Add(fhMCESphericity[i]);
3392
3393             for(Int_t ie = 0; ie < 7; ie++)
3394             {
3395               fhMCDispEtaDispPhi[ie][i] = new TH2F (Form("hMCDispEtaDispPhi_EBin%d_MC%s",ie,pnamess[i].Data()),
3396                                                     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]),
3397                                                     ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
3398               fhMCDispEtaDispPhi[ie][i]->SetXTitle("#sigma^{2}_{#eta #eta}");
3399               fhMCDispEtaDispPhi[ie][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
3400               outputContainer->Add(fhMCDispEtaDispPhi[ie][i]);
3401
3402               fhMCLambda0DispEta[ie][i] = new TH2F (Form("hMCLambda0DispEta_EBin%d_MC%s",ie,pnamess[i].Data()),
3403                                                     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]),
3404                                                     ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
3405               fhMCLambda0DispEta[ie][i]->SetXTitle("#lambda^{2}_{0}");
3406               fhMCLambda0DispEta[ie][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
3407               outputContainer->Add(fhMCLambda0DispEta[ie][i]);
3408
3409               fhMCLambda0DispPhi[ie][i] = new TH2F (Form("hMCLambda0DispPhi_EBin%d_MC%s",ie,pnamess[i].Data()),
3410                                                     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]),
3411                                                     ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
3412               fhMCLambda0DispPhi[ie][i]->SetXTitle("#lambda^{2}_{0}");
3413               fhMCLambda0DispPhi[ie][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
3414               outputContainer->Add(fhMCLambda0DispPhi[ie][i]);
3415             }
3416           }
3417         }
3418       }// loop
3419
3420       if(!GetReader()->IsEmbeddedClusterSelectionOn())
3421       {
3422         fhMCPhotonELambda0NoOverlap  = new TH2F("hELambda0_MCPhoton_NoOverlap",
3423                                                 "cluster from Photon : E vs #lambda_{0}^{2}",
3424                                                 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3425         fhMCPhotonELambda0NoOverlap->SetYTitle("#lambda_{0}^{2}");
3426         fhMCPhotonELambda0NoOverlap->SetXTitle("#it{E} (GeV)");
3427         outputContainer->Add(fhMCPhotonELambda0NoOverlap) ;
3428
3429         fhMCPhotonELambda0TwoOverlap  = new TH2F("hELambda0_MCPhoton_TwoOverlap",
3430                                                  "cluster from Photon : E vs #lambda_{0}^{2}",
3431                                                  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3432         fhMCPhotonELambda0TwoOverlap->SetYTitle("#lambda_{0}^{2}");
3433         fhMCPhotonELambda0TwoOverlap->SetXTitle("#it{E} (GeV)");
3434         outputContainer->Add(fhMCPhotonELambda0TwoOverlap) ;
3435
3436         fhMCPhotonELambda0NOverlap  = new TH2F("hELambda0_MCPhoton_NOverlap",
3437                                                "cluster from Photon : E vs #lambda_{0}^{2}",
3438                                                nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3439         fhMCPhotonELambda0NOverlap->SetYTitle("#lambda_{0}^{2}");
3440         fhMCPhotonELambda0NOverlap->SetXTitle("#it{E} (GeV)");
3441         outputContainer->Add(fhMCPhotonELambda0NOverlap) ;
3442
3443       } //No embedding
3444
3445       if(GetReader()->IsEmbeddedClusterSelectionOn())
3446       {
3447
3448         fhEmbeddedSignalFractionEnergy  = new TH2F("hEmbeddedSignal_FractionEnergy",
3449                                                    "Energy Fraction of embedded signal versus cluster energy",
3450                                                    nptbins,ptmin,ptmax,100,0.,1.);
3451         fhEmbeddedSignalFractionEnergy->SetYTitle("Fraction");
3452         fhEmbeddedSignalFractionEnergy->SetXTitle("#it{E} (GeV)");
3453         outputContainer->Add(fhEmbeddedSignalFractionEnergy) ;
3454
3455         fhEmbedPhotonELambda0FullSignal  = new TH2F("hELambda0_EmbedPhoton_FullSignal",
3456                                                     "cluster from Photon embedded with more than 90% energy in cluster : E vs #lambda_{0}^{2}",
3457                                                     nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3458         fhEmbedPhotonELambda0FullSignal->SetYTitle("#lambda_{0}^{2}");
3459         fhEmbedPhotonELambda0FullSignal->SetXTitle("#it{E} (GeV)");
3460         outputContainer->Add(fhEmbedPhotonELambda0FullSignal) ;
3461
3462         fhEmbedPhotonELambda0MostlySignal  = new TH2F("hELambda0_EmbedPhoton_MostlySignal",
3463                                                       "cluster from Photon embedded with 50% to 90% energy in cluster : E vs #lambda_{0}^{2}",
3464                                                       nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3465         fhEmbedPhotonELambda0MostlySignal->SetYTitle("#lambda_{0}^{2}");
3466         fhEmbedPhotonELambda0MostlySignal->SetXTitle("#it{E} (GeV)");
3467         outputContainer->Add(fhEmbedPhotonELambda0MostlySignal) ;
3468
3469         fhEmbedPhotonELambda0MostlyBkg  = new TH2F("hELambda0_EmbedPhoton_MostlyBkg",
3470                                                    "cluster from Photon embedded with 10% to 50% energy in cluster : E vs #lambda_{0}^{2}",
3471                                                    nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3472         fhEmbedPhotonELambda0MostlyBkg->SetYTitle("#lambda_{0}^{2}");
3473         fhEmbedPhotonELambda0MostlyBkg->SetXTitle("#it{E} (GeV)");
3474         outputContainer->Add(fhEmbedPhotonELambda0MostlyBkg) ;
3475
3476         fhEmbedPhotonELambda0FullBkg  = new TH2F("hELambda0_EmbedPhoton_FullBkg",
3477                                                  "cluster from Photonm embedded with 0% to 10% energy in cluster : E vs #lambda_{0}^{2}",
3478                                                  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3479         fhEmbedPhotonELambda0FullBkg->SetYTitle("#lambda_{0}^{2}");
3480         fhEmbedPhotonELambda0FullBkg->SetXTitle("#it{E} (GeV)");
3481         outputContainer->Add(fhEmbedPhotonELambda0FullBkg) ;
3482
3483         fhEmbedPi0ELambda0FullSignal  = new TH2F("hELambda0_EmbedPi0_FullSignal",
3484                                                  "cluster from Pi0 embedded with more than 90% energy in cluster : E vs #lambda_{0}^{2}",
3485                                                  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3486         fhEmbedPi0ELambda0FullSignal->SetYTitle("#lambda_{0}^{2}");
3487         fhEmbedPi0ELambda0FullSignal->SetXTitle("#it{E} (GeV)");
3488         outputContainer->Add(fhEmbedPi0ELambda0FullSignal) ;
3489
3490         fhEmbedPi0ELambda0MostlySignal  = new TH2F("hELambda0_EmbedPi0_MostlySignal",
3491                                                    "cluster from Pi0 embedded with 50% to 90% energy in cluster : E vs #lambda_{0}^{2}",
3492                                                    nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3493         fhEmbedPi0ELambda0MostlySignal->SetYTitle("#lambda_{0}^{2}");
3494         fhEmbedPi0ELambda0MostlySignal->SetXTitle("#it{E} (GeV)");
3495         outputContainer->Add(fhEmbedPi0ELambda0MostlySignal) ;
3496
3497         fhEmbedPi0ELambda0MostlyBkg  = new TH2F("hELambda0_EmbedPi0_MostlyBkg",
3498                                                 "cluster from Pi0 embedded with 10% to 50% energy in cluster : E vs #lambda_{0}^{2}",
3499                                                 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3500         fhEmbedPi0ELambda0MostlyBkg->SetYTitle("#lambda_{0}^{2}");
3501         fhEmbedPi0ELambda0MostlyBkg->SetXTitle("#it{E} (GeV)");
3502         outputContainer->Add(fhEmbedPi0ELambda0MostlyBkg) ;
3503
3504         fhEmbedPi0ELambda0FullBkg  = new TH2F("hELambda0_EmbedPi0_FullBkg",
3505                                               "cluster from Pi0 embedded with 0% to 10% energy in cluster : E vs #lambda_{0}^{2}",
3506                                               nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3507         fhEmbedPi0ELambda0FullBkg->SetYTitle("#lambda_{0}^{2}");
3508         fhEmbedPi0ELambda0FullBkg->SetXTitle("#it{E} (GeV)");
3509         outputContainer->Add(fhEmbedPi0ELambda0FullBkg) ;
3510
3511       }// embedded histograms
3512
3513
3514     }// Fill SS MC histograms
3515
3516   }//Histos with MC
3517
3518   return outputContainer ;
3519
3520 }
3521
3522 //_______________________
3523 void AliAnaPhoton::Init()
3524 {
3525
3526   //Init
3527   //Do some checks
3528   if(fCalorimeter == "PHOS" && !GetReader()->IsPHOSSwitchedOn() && NewOutputAOD())
3529   {
3530     printf("AliAnaPhoton::Init() - !!STOP: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!\n");
3531     abort();
3532   }
3533   else  if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn() && NewOutputAOD())
3534   {
3535     printf("AliAnaPhoton::Init() - !!STOP: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!\n");
3536     abort();
3537   }
3538
3539   if(GetReader()->GetDataType() == AliCaloTrackReader::kMC) GetCaloPID()->SwitchOnBayesian();
3540
3541 }
3542
3543 //____________________________________________________________________________
3544 void AliAnaPhoton::InitParameters()
3545 {
3546
3547   //Initialize the parameters of the analysis.
3548   AddToHistogramsName("AnaPhoton_");
3549
3550   fCalorimeter = "EMCAL" ;
3551   fMinDist     = 2.;
3552   fMinDist2    = 4.;
3553   fMinDist3    = 5.;
3554
3555   fTimeCutMin  =-1000000;
3556   fTimeCutMax  = 1000000;
3557   fNCellsCut   = 0;
3558
3559   fRejectTrackMatch       = kTRUE ;
3560
3561 }
3562
3563 //__________________________________________________________________
3564 void  AliAnaPhoton::MakeAnalysisFillAOD()
3565 {
3566   //Do photon analysis and fill aods
3567
3568   //Get the vertex
3569   Double_t v[3] = {0,0,0}; //vertex ;
3570   GetReader()->GetVertex(v);
3571
3572   //Select the Calorimeter of the photon
3573   TObjArray * pl = 0x0;
3574   AliVCaloCells* cells    = 0;
3575   if      (fCalorimeter == "PHOS" )
3576   {
3577     pl    = GetPHOSClusters();
3578     cells = GetPHOSCells();
3579   }
3580   else if (fCalorimeter == "EMCAL")
3581   {
3582     pl    = GetEMCALClusters();
3583     cells = GetEMCALCells();
3584   }
3585
3586   if(!pl)
3587   {
3588     Info("MakeAnalysisFillAOD","TObjArray with %s clusters is NULL!\n",fCalorimeter.Data());
3589     return;
3590   }
3591
3592   FillPileUpHistogramsPerEvent();
3593
3594   // Loop on raw clusters before filtering in the reader and fill control histogram
3595   if((GetReader()->GetEMCALClusterListName()=="" && fCalorimeter=="EMCAL") || fCalorimeter=="PHOS")
3596   {
3597     for(Int_t iclus = 0; iclus < GetReader()->GetInputEvent()->GetNumberOfCaloClusters(); iclus++ )
3598     {
3599       AliVCluster * clus = GetReader()->GetInputEvent()->GetCaloCluster(iclus);
3600       if     (fCalorimeter == "PHOS"  && clus->IsPHOS()  && clus->E() > GetReader()->GetPHOSPtMin() ) fhClusterCuts[0]->Fill(clus->E());
3601       else if(fCalorimeter == "EMCAL" && clus->IsEMCAL() && clus->E() > GetReader()->GetEMCALPtMin()) fhClusterCuts[0]->Fill(clus->E());
3602     }
3603   }
3604   else
3605   { // reclusterized
3606     TClonesArray * clusterList = 0;
3607
3608     if(GetReader()->GetInputEvent()->FindListObject(GetReader()->GetEMCALClusterListName()))
3609       clusterList = dynamic_cast<TClonesArray*> (GetReader()->GetInputEvent()->FindListObject(GetReader()->GetEMCALClusterListName()));
3610     else if(GetReader()->GetOutputEvent())
3611       clusterList = dynamic_cast<TClonesArray*> (GetReader()->GetOutputEvent()->FindListObject(GetReader()->GetEMCALClusterListName()));
3612
3613     if(clusterList)
3614     {
3615       Int_t nclusters = clusterList->GetEntriesFast();
3616       for (Int_t iclus =  0; iclus <  nclusters; iclus++)
3617       {
3618         AliVCluster * clus = dynamic_cast<AliVCluster*> (clusterList->At(iclus));
3619         if(clus)fhClusterCuts[0]->Fill(clus->E());
3620       }
3621     }
3622   }
3623
3624   // Fill some trigger related histograms
3625   Int_t  idTrig = GetReader()->GetTriggerClusterIndex();
3626   Bool_t exotic = GetReader()->IsExoticEvent();
3627   Bool_t bad    = GetReader()->IsBadCellTriggerEvent();
3628
3629   if( fFillEMCALBCHistograms && fCalorimeter=="EMCAL" &&
3630      ( bad || exotic )  && idTrig >= 0 && !GetReader()->AreBadTriggerEventsRemoved())
3631   {
3632     //    printf("Index %d, Id %d,  bad %d, exo %d\n",
3633     //           GetReader()->GetTriggerClusterIndex(),
3634     //           GetReader()->GetTriggerClusterId(),
3635     //           GetReader()->IsBadCellTriggerEvent(),
3636     //           GetReader()->IsExoticEvent() );
3637
3638     TClonesArray * clusterList = 0;
3639     TString  clusterListName   = GetReader()->GetEMCALClusterListName();
3640     if     (GetReader()->GetInputEvent()->FindListObject(clusterListName))
3641       clusterList = dynamic_cast<TClonesArray*> (GetReader()->GetInputEvent() ->FindListObject(clusterListName));
3642     else if(GetReader()->GetOutputEvent())
3643       clusterList = dynamic_cast<TClonesArray*> (GetReader()->GetOutputEvent()->FindListObject(clusterListName));
3644
3645     AliVCluster  *  badClusTrig = 0;
3646     if(clusterList) badClusTrig = (AliVCluster*) clusterList->At(idTrig);
3647     else            badClusTrig = GetReader()->GetInputEvent()->GetCaloCluster(idTrig);
3648
3649     if(!badClusTrig)
3650       printf("AliAnaPhoton::MakeAnalysisFillAOD() - No cluster (bad-exotic trigger) found with requested index %d \n",idTrig);
3651     else
3652     {
3653       TLorentzVector momBadClus;
3654
3655       badClusTrig->GetMomentum(momBadClus,GetVertex(0));
3656
3657       Float_t etaclusterBad = momBadClus.Eta();
3658       Float_t phiclusterBad = momBadClus.Phi();
3659       if( phiclusterBad < 0 ) phiclusterBad+=TMath::TwoPi();
3660       Float_t tofclusterBad = badClusTrig->GetTOF()*1.e9;
3661       Float_t eclusterBad   = badClusTrig->E();
3662
3663       if( bad && exotic )
3664       {
3665         if(GetReader()->IsTriggerMatched())
3666         {
3667           fhEtaPhiTriggerEMCALBCBadExoticCluster->Fill(etaclusterBad, phiclusterBad);
3668           fhTimeTriggerEMCALBCBadExoticCluster  ->Fill(eclusterBad,   tofclusterBad);
3669         }
3670         else
3671         {
3672           fhEtaPhiTriggerEMCALBCUMBadExoticCluster->Fill(etaclusterBad, phiclusterBad);
3673           fhTimeTriggerEMCALBCUMBadExoticCluster  ->Fill(eclusterBad,   tofclusterBad);
3674         }
3675       }
3676       else if( bad && !exotic )
3677       {
3678         if(GetReader()->IsTriggerMatched())
3679         {
3680           fhEtaPhiTriggerEMCALBCBadCluster->Fill(etaclusterBad, phiclusterBad);
3681           fhTimeTriggerEMCALBCBadCluster  ->Fill(eclusterBad,   tofclusterBad);
3682         }
3683         else
3684         {
3685           fhEtaPhiTriggerEMCALBCUMBadCluster->Fill(etaclusterBad, phiclusterBad);
3686           fhTimeTriggerEMCALBCUMBadCluster  ->Fill(eclusterBad,   tofclusterBad);
3687         }
3688       }// Bad cluster trigger
3689       else if( !bad && exotic )
3690       {
3691         if(GetReader()->IsTriggerMatched())
3692         {
3693           fhEtaPhiTriggerEMCALBCExoticCluster->Fill(etaclusterBad, phiclusterBad);
3694           fhTimeTriggerEMCALBCExoticCluster  ->Fill(eclusterBad, tofclusterBad);
3695         }
3696         else
3697         {
3698           fhEtaPhiTriggerEMCALBCUMExoticCluster->Fill(etaclusterBad, phiclusterBad);
3699           fhTimeTriggerEMCALBCUMExoticCluster  ->Fill(eclusterBad, tofclusterBad);
3700         }
3701       }
3702     }// cluster exists
3703   } // study bad/exotic trigger BC
3704
3705   //Init arrays, variables, get number of clusters
3706   TLorentzVector mom, mom2 ;
3707   Int_t nCaloClusters = pl->GetEntriesFast();
3708
3709   if(GetDebug() > 0) printf("AliAnaPhoton::MakeAnalysisFillAOD() - input %s cluster entries %d\n", fCalorimeter.Data(), nCaloClusters);
3710
3711   //----------------------------------------------------
3712   // Fill AOD with PHOS/EMCAL AliAODPWG4Particle objects
3713   //----------------------------------------------------
3714   // Loop on clusters
3715   for(Int_t icalo = 0; icalo < nCaloClusters; icalo++)
3716   {
3717           AliVCluster * calo =  (AliVCluster*) (pl->At(icalo));
3718     //printf("calo %d, %f\n",icalo,calo->E());
3719
3720     //Get the index where the cluster comes, to retrieve the corresponding vertex
3721     Int_t evtIndex = 0 ;
3722     if (GetMixedEvent())
3723     {
3724       evtIndex=GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ;
3725       //Get the vertex and check it is not too large in z
3726       if(TMath::Abs(GetVertex(evtIndex)[2])> GetZvertexCut()) continue;
3727     }
3728
3729     //Cluster selection, not charged, with photon id and in fiducial cut
3730     if(GetReader()->GetDataType() != AliCaloTrackReader::kMC)
3731     {
3732       calo->GetMomentum(mom,GetVertex(evtIndex)) ;
3733     }//Assume that come from vertex in straight line
3734     else
3735     {
3736       Double_t vertex[]={0,0,0};
3737       calo->GetMomentum(mom,vertex) ;
3738     }
3739
3740     //--------------------------------------
3741     // Cluster selection
3742     //--------------------------------------
3743     Int_t nMaxima = GetCaloUtils()->GetNumberOfLocalMaxima(calo, cells); // NLM
3744     if(!ClusterSelected(calo,mom,nMaxima)) continue;
3745
3746     //----------------------------
3747     //Create AOD for analysis
3748     //----------------------------
3749     AliAODPWG4Particle aodph = AliAODPWG4Particle(mom);
3750
3751     //...............................................
3752     //Set the indeces of the original caloclusters (MC, ID), and calorimeter
3753     Int_t label = calo->GetLabel();
3754     aodph.SetLabel(label);
3755     aodph.SetCaloLabel(calo->GetID(),-1);
3756     aodph.SetDetector(fCalorimeter);
3757     //printf("Index %d, Id %d, iaod %d\n",icalo, calo->GetID(),GetOutputAODBranch()->GetEntriesFast());
3758
3759     //...............................................
3760     //Set bad channel distance bit
3761     Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
3762     if     (distBad > fMinDist3) aodph.SetDistToBad(2) ;
3763     else if(distBad > fMinDist2) aodph.SetDistToBad(1) ;
3764     else                         aodph.SetDistToBad(0) ;
3765     //printf("DistBad %f Bit %d\n",distBad, aodph.DistToBad());
3766
3767     //--------------------------------------------------------------------------------------
3768     // Play with the MC stack if available
3769     //--------------------------------------------------------------------------------------
3770
3771     //Check origin of the candidates
3772     Int_t tag = -1;
3773
3774     if(IsDataMC())
3775     {
3776       tag = GetMCAnalysisUtils()->CheckOrigin(calo->GetLabels(),calo->GetNLabels(),GetReader());
3777       aodph.SetTag(tag);
3778
3779       if(GetDebug() > 0)
3780         printf("AliAnaPhoton::MakeAnalysisFillAOD() - Origin of candidate, bit map %d\n",aodph.GetTag());
3781     }//Work with stack also
3782
3783
3784     //--------------------------------------------------------------------------------------
3785     //Fill some shower shape histograms before PID is applied
3786     //--------------------------------------------------------------------------------------
3787
3788     FillShowerShapeHistograms(calo,tag);
3789
3790     //-------------------------------------
3791     //PID selection or bit setting
3792     //-------------------------------------
3793
3794     //...............................................
3795     // Data, PID check on
3796     if(IsCaloPIDOn())
3797     {
3798       // Get most probable PID, 2 options check bayesian PID weights or redo PID
3799       // By default, redo PID
3800
3801       aodph.SetIdentifiedParticleType(GetCaloPID()->GetIdentifiedParticleType(calo));
3802
3803       if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PDG of identified particle %d\n",aodph.GetIdentifiedParticleType());
3804
3805       //If cluster does not pass pid, not photon, skip it.
3806       if(aodph.GetIdentifiedParticleType() != AliCaloPID::kPhoton) continue ;
3807
3808     }
3809
3810     //...............................................
3811     // Data, PID check off
3812     else
3813     {
3814       //Set PID bits for later selection (AliAnaPi0 for example)
3815       //GetIdentifiedParticleType already called in SetPIDBits.
3816
3817       GetCaloPID()->SetPIDBits(calo,&aodph, GetCaloUtils(),GetReader()->GetInputEvent());
3818
3819       if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PID Bits set \n");
3820     }
3821
3822     if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - Photon selection cuts passed: pT %3.2f, pdg %d\n",
3823                               aodph.Pt(), aodph.GetIdentifiedParticleType());
3824
3825     fhClusterCuts[9]->Fill(calo->E());
3826
3827     Int_t   nSM  = GetModuleNumber(calo);
3828     if(nSM < GetCaloUtils()->GetNumberOfSuperModulesUsed() && nSM >=0)
3829     {
3830       fhEPhotonSM ->Fill(mom.E (),nSM);
3831       fhPtPhotonSM->Fill(mom.Pt(),nSM);
3832     }
3833
3834     fhNLocMax->Fill(calo->E(),nMaxima);
3835
3836     // Matching after cuts
3837     if(fFillTMHisto) FillTrackMatchingResidualHistograms(calo,1);
3838
3839     // Fill histograms to undertand pile-up before other cuts applied
3840     // Remember to relax time cuts in the reader
3841     FillPileUpHistograms(calo->E(),mom.Pt(),calo->GetTOF()*1e9);
3842
3843     // Add number of local maxima to AOD, method name in AOD to be FIXED
3844     aodph.SetFiducialArea(nMaxima);
3845
3846     if(fFillEMCALBCHistograms && fCalorimeter=="EMCAL" && idTrig >= 0)
3847     {
3848       Double_t calotof   = calo->GetTOF()*1e9;
3849       Float_t  calotofUS = TMath::Abs(calotof);
3850       Float_t phicluster = aodph.Phi();
3851       if(phicluster < 0) phicluster+=TMath::TwoPi();
3852
3853       if(calo->E() > 2)
3854       {
3855         if      (calotofUS < 25) fhEtaPhiPhotonEMCALBC0->Fill(aodph.Eta(), phicluster);
3856         else if (calotofUS < 75) fhEtaPhiPhotonEMCALBC1->Fill(aodph.Eta(), phicluster);
3857         else                     fhEtaPhiPhotonEMCALBCN->Fill(aodph.Eta(), phicluster);
3858       }
3859
3860       Int_t bc = GetReader()->GetTriggerClusterBC();
3861       Int_t histoBC = bc-5;
3862       if(GetReader()->AreBadTriggerEventsRemoved()) histoBC=0; // histograms created only for one BC since the others where rejected
3863       if(TMath::Abs(bc) < 6 && !GetReader()->IsBadCellTriggerEvent() && !GetReader()->IsExoticEvent())
3864       {
3865         if(GetReader()->IsTriggerMatched())
3866         {
3867           if(calo->E() > 2) fhEtaPhiPhotonTriggerEMCALBC[histoBC]->Fill(aodph.Eta(), phicluster);
3868           fhTimePhotonTriggerEMCALBC[histoBC]->Fill(calo->E(), calotof);
3869           if(GetReader()->IsPileUpFromSPD()) fhTimePhotonTriggerEMCALBCPileUpSPD[histoBC]->Fill(calo->E(), calotof);
3870         }
3871         else
3872         {
3873           if(calo->E() > 2) fhEtaPhiPhotonTriggerEMCALBCUM[histoBC]->Fill(aodph.Eta(), phicluster);
3874           fhTimePhotonTriggerEMCALBCUM[histoBC]->Fill(calo->E(), calotof);
3875
3876           if(bc==0)
3877           {
3878             if(GetReader()->IsTriggerMatchedOpenCuts(0)) fhTimePhotonTriggerEMCALBC0UMReMatchOpenTime   ->Fill(calo->E(), calotof);
3879             if(GetReader()->IsTriggerMatchedOpenCuts(1)) fhTimePhotonTriggerEMCALBC0UMReMatchCheckNeigh ->Fill(calo->E(), calotof);
3880             if(GetReader()->IsTriggerMatchedOpenCuts(2)) fhTimePhotonTriggerEMCALBC0UMReMatchBoth       ->Fill(calo->E(), calotof);
3881           }
3882         }
3883       }
3884       else if(TMath::Abs(bc) >= 6)
3885         printf("AliAnaPhoton::MakeAnalysisFillAOD() - Trigger BC not expected = %d\n",bc);
3886     }
3887
3888     //Add AOD with photon object to aod branch
3889     AddAODParticle(aodph);
3890
3891   }//loop
3892
3893   if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD()  End fill AODs, with %d entries \n",GetOutputAODBranch()->GetEntriesFast());
3894
3895 }
3896
3897 //__________________________________________________________________
3898 void  AliAnaPhoton::MakeAnalysisFillHistograms()
3899 {
3900   //Fill histograms
3901
3902   //In case of simulated data, fill acceptance histograms
3903   if(IsDataMC()) FillAcceptanceHistograms();
3904
3905   // Get vertex
3906   Double_t v[3] = {0,0,0}; //vertex ;
3907   GetReader()->GetVertex(v);
3908   //fhVertex->Fill(v[0],v[1],v[2]);
3909   if(TMath::Abs(v[2]) > GetZvertexCut()) return ; // done elsewhere for Single Event analysis, but there for mixed event
3910
3911   //----------------------------------
3912   //Loop on stored AOD photons
3913   Int_t naod = GetOutputAODBranch()->GetEntriesFast();
3914   if(GetDebug() > 0) printf("AliAnaPhoton::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod);
3915
3916   Float_t cen = GetEventCentrality();
3917   // printf("++++++++++ GetEventCentrality() %f\n",cen);
3918
3919   Float_t ep  = GetEventPlaneAngle();
3920
3921   for(Int_t iaod = 0; iaod < naod ; iaod++)
3922   {
3923     AliAODPWG4Particle* ph =  (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
3924     Int_t pdg = ph->GetIdentifiedParticleType();
3925
3926     if(GetDebug() > 3)
3927       printf("AliAnaPhoton::MakeAnalysisFillHistograms() - PDG %d, MC TAG %d, Calorimeter %s\n",
3928              ph->GetIdentifiedParticleType(),ph->GetTag(), (ph->GetDetector()).Data()) ;
3929
3930     //If PID used, fill histos with photons in Calorimeter fCalorimeter
3931     if(IsCaloPIDOn() && pdg != AliCaloPID::kPhoton) continue;
3932     if(ph->GetDetector() != fCalorimeter) continue;
3933
3934     if(GetDebug() > 2)
3935       printf("AliAnaPhoton::MakeAnalysisFillHistograms() - ID Photon: pt %f, phi %f, eta %f\n", ph->Pt(),ph->Phi(),ph->Eta()) ;
3936
3937     //................................
3938     //Fill photon histograms
3939     Float_t ptcluster  = ph->Pt();
3940     Float_t phicluster = ph->Phi();
3941     Float_t etacluster = ph->Eta();
3942     Float_t ecluster   = ph->E();
3943
3944     fhEPhoton   ->Fill(ecluster);
3945     fhPtPhoton  ->Fill(ptcluster);
3946     fhPhiPhoton ->Fill(ptcluster,phicluster);
3947     fhEtaPhoton ->Fill(ptcluster,etacluster);
3948     if     (ecluster   > 0.5) fhEtaPhiPhoton  ->Fill(etacluster, phicluster);
3949     else if(GetMinPt() < 0.5) fhEtaPhi05Photon->Fill(etacluster, phicluster);
3950
3951     fhPtCentralityPhoton ->Fill(ptcluster,cen) ;
3952     fhPtEventPlanePhoton ->Fill(ptcluster,ep ) ;
3953
3954     //Get original cluster, to recover some information
3955     AliVCaloCells* cells    = 0;
3956     TObjArray * clusters    = 0;
3957     if(fCalorimeter == "EMCAL")
3958     {
3959       cells    = GetEMCALCells();
3960       clusters = GetEMCALClusters();
3961     }
3962     else
3963     {
3964       cells    = GetPHOSCells();
3965       clusters = GetPHOSClusters();
3966     }
3967
3968     Int_t iclus = -1;
3969     AliVCluster *cluster = FindCluster(clusters,ph->GetCaloLabel(0),iclus);
3970     if(cluster)
3971     {
3972       Float_t maxCellFraction = 0;
3973       Int_t absID = GetCaloUtils()->GetMaxEnergyCell(cells, cluster, maxCellFraction);
3974       if( absID < 0 ) AliFatal("Wrong absID");
3975
3976       // Control histograms
3977       fhMaxCellDiffClusterE->Fill(ph->E() ,maxCellFraction);
3978       fhNCellsE            ->Fill(ph->E() ,cluster->GetNCells());
3979       fhTimePt             ->Fill(ph->Pt(),cluster->GetTOF()*1.e9);
3980
3981       if(cells)
3982       {
3983         for(Int_t icell = 0; icell <  cluster->GetNCells(); icell++)
3984           fhCellsE->Fill(ph->E(),cells->GetCellAmplitude(cluster->GetCellsAbsId()[icell]));
3985       }
3986     }
3987
3988     //.......................................
3989     //Play with the MC data if available
3990     if(IsDataMC())
3991     {
3992       if(GetDebug()>0)
3993       {
3994         if(GetReader()->ReadStack() && !GetMCStack())
3995         {
3996           printf("AliAnaPhoton::MakeAnalysisFillHistograms() - Stack not available, is the MC handler called?\n");
3997         }
3998         else if(GetReader()->ReadAODMCParticles() && !GetReader()->GetAODMCParticles())
3999         {
4000           printf("AliAnaPhoton::MakeAnalysisFillHistograms() -  Standard MCParticles not available!\n");
4001         }
4002       }
4003
4004       //....................................................................
4005       // Access MC information in stack if requested, check that it exists.
4006       Int_t label =ph->GetLabel();
4007
4008       if(label < 0)
4009       {
4010         if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** bad label ***:  label %d \n", label);
4011         continue;
4012       }
4013
4014       Float_t eprim   = 0;
4015       Float_t ptprim  = 0;
4016       Bool_t ok = kFALSE;
4017       TLorentzVector primary = GetMCAnalysisUtils()->GetMother(label,GetReader(),ok);
4018       if(ok)
4019       {
4020         eprim   = primary.Energy();
4021         ptprim  = primary.Pt();
4022       }
4023
4024       Int_t tag =ph->GetTag();
4025       Int_t mcParticleTag = -1;
4026       if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) && fhMCE[kmcPhoton])
4027       {
4028         fhMCE  [kmcPhoton] ->Fill(ecluster);
4029         fhMCPt [kmcPhoton] ->Fill(ptcluster);
4030         fhMCPhi[kmcPhoton] ->Fill(ecluster,phicluster);
4031         fhMCEta[kmcPhoton] ->Fill(ecluster,etacluster);
4032
4033         fhMC2E     [kmcPhoton] ->Fill(ecluster, eprim);
4034         fhMC2Pt    [kmcPhoton] ->Fill(ptcluster, ptprim);
4035         fhMCDeltaE [kmcPhoton] ->Fill(ecluster,eprim-ecluster);
4036         fhMCDeltaPt[kmcPhoton] ->Fill(ptcluster,ptprim-ptcluster);
4037
4038         if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion) &&
4039            GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton)     &&
4040            fhMCE[kmcConversion])
4041         {
4042           fhMCE  [kmcConversion] ->Fill(ecluster);
4043           fhMCPt [kmcConversion] ->Fill(ptcluster);
4044           fhMCPhi[kmcConversion] ->Fill(ecluster,phicluster);
4045           fhMCEta[kmcConversion] ->Fill(ecluster,etacluster);
4046
4047           fhMC2E     [kmcConversion] ->Fill(ecluster, eprim);
4048           fhMC2Pt    [kmcConversion] ->Fill(ptcluster, ptprim);
4049           fhMCDeltaE [kmcConversion] ->Fill(ecluster,eprim-ecluster);
4050           fhMCDeltaPt[kmcConversion] ->Fill(ptcluster,ptprim-ptcluster);
4051         }
4052
4053         if     (GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt))
4054         {
4055           mcParticleTag = kmcPrompt;
4056         }
4057         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation))
4058         {
4059           mcParticleTag = kmcFragmentation;
4060         }
4061         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCISR)&& fhMCE[kmcISR])
4062         {
4063           mcParticleTag = kmcISR;
4064         }
4065         else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay) &&
4066                 !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0))
4067         {
4068           mcParticleTag = kmcPi0Decay;
4069         }
4070         else if((( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay) &&
4071                   !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta)        ) ||
4072                  GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay) ))
4073         {
4074           mcParticleTag = kmcOtherDecay;
4075         }
4076         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0))
4077         {
4078           mcParticleTag = kmcPi0;
4079         }
4080         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta))
4081         {
4082           mcParticleTag = kmcEta;
4083         }
4084       }
4085       else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiNeutron))
4086       {
4087         mcParticleTag = kmcAntiNeutron;
4088       }
4089       else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiProton))
4090       {
4091         mcParticleTag = kmcAntiProton;
4092       }
4093       else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron))
4094       {
4095         mcParticleTag = kmcElectron;
4096       }
4097       else if( fhMCE[kmcOther])
4098       {
4099         mcParticleTag = kmcOther;
4100
4101         //               printf(" AliAnaPhoton::MakeAnalysisFillHistograms() - Label %d, pT %2.3f Unknown, bits set: ",
4102         //                                      ph->GetLabel(),ph->Pt());
4103         //                for(Int_t i = 0; i < 20; i++) {
4104         //                        if(GetMCAnalysisUtils()->CheckTagBit(tag,i)) printf(" %d, ",i);
4105         //                }
4106         //                printf("\n");
4107
4108       }
4109
4110       if(mcParticleTag >= 0 && fhMCE  [mcParticleTag])
4111       {
4112         fhMCE  [mcParticleTag]->Fill(ecluster);
4113         fhMCPt [mcParticleTag]->Fill(ptcluster);
4114         fhMCPhi[mcParticleTag]->Fill(ecluster,phicluster);
4115         fhMCEta[mcParticleTag]->Fill(ecluster,etacluster);
4116
4117         fhMC2E     [mcParticleTag]->Fill(ecluster, eprim);
4118         fhMC2Pt    [mcParticleTag]->Fill(ptcluster, ptprim);
4119         fhMCDeltaE [mcParticleTag]->Fill(ecluster,eprim-ecluster);
4120         fhMCDeltaPt[mcParticleTag]->Fill(ptcluster,ptprim-ptcluster);
4121       }
4122     }//Histograms with MC
4123
4124   }// aod loop
4125
4126 }
4127
4128
4129 //__________________________________________________________________
4130 void AliAnaPhoton::Print(const Option_t * opt) const
4131 {
4132   //Print some relevant parameters set for the analysis
4133
4134   if(! opt)
4135     return;
4136
4137   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
4138   AliAnaCaloTrackCorrBaseClass::Print(" ");
4139
4140   printf("Calorimeter            =     %s\n", fCalorimeter.Data()) ;
4141   printf("Min Distance to Bad Channel   = %2.1f\n",fMinDist);
4142   printf("Min Distance to Bad Channel 2 = %2.1f\n",fMinDist2);
4143   printf("Min Distance to Bad Channel 3 = %2.1f\n",fMinDist3);
4144   printf("Reject clusters with a track matched = %d\n",fRejectTrackMatch);
4145   printf("Time Cut: %3.1f < TOF  < %3.1f\n", fTimeCutMin, fTimeCutMax);
4146   printf("Number of cells in cluster is        > %d \n", fNCellsCut);
4147   printf("    \n") ;
4148
4149 }