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