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