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