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