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