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