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),
  66     // Histograms
  67     fhNCellsE(0),                 fhCellsE(0),   // Control histograms
  68     fhMaxCellDiffClusterE(0),     fhTimeE(0),    // Control histograms
  69     fhEPhoton(0),                 fhPtPhoton(0),
  70     fhPhiPhoton(0),               fhEtaPhoton(0),
  71     fhEtaPhiPhoton(0),            fhEtaPhi05Photon(0),
  72
  73     // Shower shape histograms
  74     fhNLocMax(0),
  75     fhDispE(0),                   fhLam0E(0),                   fhLam1E(0),
  76     fhDispETRD(0),                fhLam0ETRD(0),                fhLam1ETRD(0),
  77     fhDispETM(0),                 fhLam0ETM(0),                 fhLam1ETM(0),
  78     fhDispETMTRD(0),              fhLam0ETMTRD(0),              fhLam1ETMTRD(0),
  79
  80     fhNCellsLam0LowE(0),          fhNCellsLam1LowE(0),          fhNCellsDispLowE(0),
  81     fhNCellsLam0HighE(0),         fhNCellsLam1HighE(0),         fhNCellsDispHighE(0),
  82
  83     fhEtaLam0LowE(0),             fhPhiLam0LowE(0),
  84     fhEtaLam0HighE(0),            fhPhiLam0HighE(0),
  85     fhLam0DispLowE(0),            fhLam0DispHighE(0),
  86     fhLam1Lam0LowE(0),            fhLam1Lam0HighE(0),
  87     fhDispLam1LowE(0),            fhDispLam1HighE(0),
  88     fhDispEtaE(0),                fhDispPhiE(0),
  89     fhSumEtaE(0),                 fhSumPhiE(0),                 fhSumEtaPhiE(0),
  90     fhDispEtaPhiDiffE(0),         fhSphericityE(0),
  91     fhDispSumEtaDiffE(0),         fhDispSumPhiDiffE(0),
  92
  93     // MC histograms
  94     fhMCPhotonELambda0NoOverlap(0),       fhMCPhotonELambda0TwoOverlap(0),      fhMCPhotonELambda0NOverlap(0),
  95     // Embedding
  96     fhEmbeddedSignalFractionEnergy(0),
  97     fhEmbedPhotonELambda0FullSignal(0),   fhEmbedPhotonELambda0MostlySignal(0),
  98     fhEmbedPhotonELambda0MostlyBkg(0),    fhEmbedPhotonELambda0FullBkg(0),
  99     fhEmbedPi0ELambda0FullSignal(0),      fhEmbedPi0ELambda0MostlySignal(0),
 100     fhEmbedPi0ELambda0MostlyBkg(0),       fhEmbedPi0ELambda0FullBkg(0),
 101     // PileUp
 102     fhTimeENoCut(0),                      fhTimeESPD(0),        fhTimeESPDMulti(0),
 103     fhTimeNPileUpVertSPD(0),              fhTimeNPileUpVertTrack(0),
 104     fhTimeNPileUpVertContributors(0),
 105     fhTimePileUpMainVertexZDistance(0),   fhTimePileUpMainVertexZDiamond(0),
 106     fhClusterMultSPDPileUp(),             fhClusterMultNoPileUp(),
 107     fhEtaPhiBC0(0),  fhEtaPhiBCPlus(0),   fhEtaPhiBCMinus(0),
 108     fhEtaPhiBC0PileUpSPD(0),
 109     fhEtaPhiBCPlusPileUpSPD(0),
 110     fhEtaPhiBCMinusPileUpSPD(0)
 111  {
 112   //default ctor
 113
 114   for(Int_t i = 0; i < 14; i++)
 115   {
 116     fhMCPt     [i] = 0;
 117     fhMCE      [i] = 0;
 118     fhMCPhi    [i] = 0;
 119     fhMCEta    [i] = 0;
 120     fhMCDeltaE [i] = 0;
 121     fhMCDeltaPt[i] = 0;
 122     fhMC2E     [i] = 0;
 123     fhMC2Pt    [i] = 0;
 124   }
 125
 126   for(Int_t i = 0; i < 7; i++)
 127   {
 128     fhPtPrimMC [i] = 0;
 129     fhEPrimMC  [i] = 0;
 130     fhPhiPrimMC[i] = 0;
 131     fhYPrimMC  [i] = 0;
 132
 133     fhPtPrimMCAcc [i] = 0;
 134     fhEPrimMCAcc  [i] = 0;
 135     fhPhiPrimMCAcc[i] = 0;
 136     fhYPrimMCAcc  [i] = 0;
 137
 138     fhDispEtaDispPhi[i] = 0;
 139     fhLambda0DispPhi[i] = 0;
 140     fhLambda0DispEta[i] = 0;
 141
 142     fhPtPileUp       [i] = 0;
 143     fhPtChargedPileUp[i] = 0;
 144     fhPtPhotonPileUp [i] = 0;
 145
 146     fhLambda0PileUp       [i] = 0;
 147     fhLambda0ChargedPileUp[i] = 0;
 148
 149     fhClusterEFracLongTimePileUp  [i] = 0;
 150
 151     fhClusterTimeDiffPileUp       [i] = 0;
 152     fhClusterTimeDiffChargedPileUp[i] = 0;
 153     fhClusterTimeDiffPhotonPileUp [i] = 0;
 154
 155     for(Int_t j = 0; j < 6; j++)
 156     {
 157       fhMCDispEtaDispPhi[i][j] = 0;
 158       fhMCLambda0DispEta[i][j] = 0;
 159       fhMCLambda0DispPhi[i][j] = 0;
 160     }
 161   }
 162
 163   for(Int_t i = 0; i < 6; i++)
 164   {
 165     fhMCELambda0    [i]                  = 0;
 166     fhMCELambda1    [i]                  = 0;
 167     fhMCEDispersion [i]                  = 0;
 168     fhMCNCellsE     [i]                  = 0;
 169     fhMCMaxCellDiffClusterE[i]           = 0;
 170     fhLambda0DispEta[i]                  = 0;
 171     fhLambda0DispPhi[i]                  = 0;
 172
 173     fhMCLambda0vsClusterMaxCellDiffE0[i] = 0;
 174     fhMCLambda0vsClusterMaxCellDiffE2[i] = 0;
 175     fhMCLambda0vsClusterMaxCellDiffE6[i] = 0;
 176     fhMCNCellsvsClusterMaxCellDiffE0 [i] = 0;
 177     fhMCNCellsvsClusterMaxCellDiffE2 [i] = 0;
 178     fhMCNCellsvsClusterMaxCellDiffE6 [i] = 0;
 179
 180     fhMCEDispEta       [i]               = 0;
 181     fhMCEDispPhi       [i]               = 0;
 182     fhMCESumEtaPhi     [i]               = 0;
 183     fhMCEDispEtaPhiDiff[i]               = 0;
 184     fhMCESphericity    [i]               = 0;
 185   }
 186
 187    for(Int_t i = 0; i < 5; i++)
 188    {
 189      fhClusterCuts[i] = 0;
 190    }
 191
 192    // Track matching residuals
 193    for(Int_t i = 0; i < 2; i++)
 194    {
 195      fhTrackMatchedDEta[i] = 0;                fhTrackMatchedDPhi[i] = 0;         fhTrackMatchedDEtaDPhi[i] = 0;
 196      fhTrackMatchedDEtaTRD[i] = 0;             fhTrackMatchedDPhiTRD[i] = 0;
 197      fhTrackMatchedDEtaMCOverlap[i] = 0;       fhTrackMatchedDPhiMCOverlap[i] = 0;
 198      fhTrackMatchedDEtaMCNoOverlap[i] = 0;     fhTrackMatchedDPhiMCNoOverlap[i] = 0;
 199      fhTrackMatchedDEtaMCConversion[i] = 0;    fhTrackMatchedDPhiMCConversion[i] = 0;
 200      fhTrackMatchedMCParticle[i] = 0;          fhTrackMatchedMCParticle[i] = 0;
 201      fhdEdx[i] = 0;                            fhEOverP[i] = 0;
 202      fhEOverPTRD[i] = 0;
 203    }
 204
 205    for(Int_t i = 0; i < 4; i++)
 206    {
 207      fhClusterMultSPDPileUp[i] = 0;
 208      fhClusterMultNoPileUp [i] = 0;
 209    }
 210
 211   //Initialize parameters
 212   InitParameters();
 213
 214 }
 215
 216 //_____________________________________________________________________________________________________
 217 Bool_t  AliAnaPhoton::ClusterSelected(AliVCluster* calo, const TLorentzVector mom, const Int_t nMaxima)
 218 {
 219   //Select clusters if they pass different cuts
 220
 221   Float_t ptcluster = mom.Pt();
 222   Float_t ecluster  = mom.E();
 223   Float_t l0cluster = calo->GetM02();
 224
 225   if(GetDebug() > 2)
 226     printf("AliAnaPhoton::ClusterSelected() Current Event %d; Before selection : E %2.2f, pT %2.2f, phi %2.2f, eta %2.2f\n",
 227            GetReader()->GetEventNumber(),
 228            ecluster,ptcluster, mom.Phi()*TMath::RadToDeg(),mom.Eta());
 229
 230   fhClusterCuts[1]->Fill(ecluster);
 231
 232   //.......................................
 233   //If too small or big energy, skip it
 234   if(ecluster < GetMinEnergy() || ecluster > GetMaxEnergy() ) return kFALSE ;
 235
 236   if(GetDebug() > 2) printf("\t Cluster %d Pass E Cut \n",calo->GetID());
 237
 238   fhClusterCuts[2]->Fill(ecluster);
 239
 240   if(fFillPileUpHistograms)
 241   {
 242     // Get the fraction of the cluster energy that carries the cell with highest energy and its absId
 243     AliVCaloCells* cells = 0;
 244     if(fCalorimeter == "EMCAL") cells = GetEMCALCells();
 245     else                        cells = GetPHOSCells();
 246
 247     Float_t maxCellFraction = 0.;
 248     Int_t absIdMax = GetCaloUtils()->GetMaxEnergyCell(cells, calo,maxCellFraction);
 249
 250     Double_t tmax  = cells->GetCellTime(absIdMax);
 251     GetCaloUtils()->RecalibrateCellTime(tmax, fCalorimeter, absIdMax,GetReader()->GetInputEvent()->GetBunchCrossNumber());
 252     tmax*=1.e9;
 253
 254     Bool_t okPhoton = kFALSE;
 255     if( GetCaloPID()->GetIdentifiedParticleType(calo)== AliCaloPID::kPhoton) okPhoton = kTRUE;
 256
 257     Bool_t matched = IsTrackMatched(calo,GetReader()->GetInputEvent());
 258     Float_t clusterLongTimeE = 0;
 259     Float_t clusterOKTimeE   = 0;
 260     //Loop on cells inside cluster
 261     for (Int_t ipos = 0; ipos < calo->GetNCells(); ipos++)
 262     {
 263       Int_t absId  = calo->GetCellsAbsId()[ipos];
 264       //if(absId!=absIdMax && cells->GetCellAmplitude(absIdMax) > 0.01)
 265       if(cells->GetCellAmplitude(absIdMax) > 0.1)
 266       {
 267         Double_t time  = cells->GetCellTime(absId);
 268         Float_t  amp   = cells->GetCellAmplitude(absId);
 269         Int_t    bc    = GetReader()->GetInputEvent()->GetBunchCrossNumber();
 270         GetCaloUtils()->GetEMCALRecoUtils()->AcceptCalibrateCell(absId,bc,amp,time,cells);
 271         time*=1e9;
 272
 273         Float_t diff = (tmax-time);
 274
 275         if(GetReader()->IsInTimeWindow(time,amp)) clusterOKTimeE   += amp;
 276         else                                      clusterLongTimeE += amp;
 277
 278         if(GetReader()->IsPileUpFromSPD())
 279         {
 280           fhClusterTimeDiffPileUp[0]->Fill(ecluster, diff);
 281           if(!matched)
 282           {
 283             fhClusterTimeDiffChargedPileUp[0]->Fill(ecluster, diff);
 284             if(okPhoton)  fhClusterTimeDiffPhotonPileUp[0]->Fill(ecluster, diff);
 285           }
 286         }
 287
 288         if(GetReader()->IsPileUpFromEMCal())
 289         {
 290           fhClusterTimeDiffPileUp[1]->Fill(ecluster, diff);
 291           if(!matched)
 292           {
 293             fhClusterTimeDiffChargedPileUp[1]->Fill(ecluster, diff);
 294             if(okPhoton)  fhClusterTimeDiffPhotonPileUp[1]->Fill(ecluster, diff);
 295           }
 296         }
 297
 298         if(GetReader()->IsPileUpFromSPDOrEMCal())
 299         {
 300           fhClusterTimeDiffPileUp[2]->Fill(ecluster, diff);
 301           if(!matched)
 302           {
 303             fhClusterTimeDiffChargedPileUp[2]->Fill(ecluster, diff);
 304             if(okPhoton)  fhClusterTimeDiffPhotonPileUp[2]->Fill(ecluster, diff);
 305           }
 306         }
 307
 308         if(GetReader()->IsPileUpFromSPDAndEMCal())
 309         {
 310           fhClusterTimeDiffPileUp[3]->Fill(ecluster, diff);
 311           if(!matched)
 312           {
 313             fhClusterTimeDiffChargedPileUp[3]->Fill(ecluster, diff);
 314             if(okPhoton)  fhClusterTimeDiffPhotonPileUp[3]->Fill(ecluster, diff);
 315           }
 316         }
 317
 318         if(GetReader()->IsPileUpFromSPDAndNotEMCal())
 319         {
 320           fhClusterTimeDiffPileUp[4]->Fill(ecluster, diff);
 321           if(!matched)
 322           {
 323             fhClusterTimeDiffChargedPileUp[4]->Fill(ecluster, diff);
 324             if(okPhoton)  fhClusterTimeDiffPhotonPileUp[4]->Fill(ecluster, diff);
 325           }
 326         }
 327
 328         if(GetReader()->IsPileUpFromEMCalAndNotSPD())
 329         {
 330           fhClusterTimeDiffPileUp[5]->Fill(ecluster, diff);
 331           if(!matched)
 332           {
 333             fhClusterTimeDiffChargedPileUp[5]->Fill(ecluster, diff);
 334             if(okPhoton)  fhClusterTimeDiffPhotonPileUp[5]->Fill(ecluster, diff);
 335           }
 336         }
 337
 338         if(GetReader()->IsPileUpFromNotSPDAndNotEMCal())
 339         {
 340           fhClusterTimeDiffPileUp[6]->Fill(ecluster, diff);
 341           if(!matched)
 342           {
 343             fhClusterTimeDiffChargedPileUp[6]->Fill(ecluster, diff);
 344             if(okPhoton)  fhClusterTimeDiffPhotonPileUp[6]->Fill(ecluster, diff);
 345           }
 346         }
 347       }// Not max
 348     }//loop
 349
 350     Float_t frac = 0;
 351     if(clusterLongTimeE+clusterOKTimeE > 0.001)
 352       frac = clusterLongTimeE/(clusterLongTimeE+clusterOKTimeE);
 353     //printf("E long %f, E OK %f, Fraction large time %f, E %f\n",clusterLongTimeE,clusterOKTimeE,frac,ecluster);
 354
 355     if(GetReader()->IsPileUpFromSPD())               {fhPtPileUp[0]->Fill(ptcluster); fhLambda0PileUp[0]->Fill(ecluster,l0cluster); fhClusterEFracLongTimePileUp[0]->Fill(ecluster,frac);}
 356     if(GetReader()->IsPileUpFromEMCal())             {fhPtPileUp[1]->Fill(ptcluster); fhLambda0PileUp[1]->Fill(ecluster,l0cluster); fhClusterEFracLongTimePileUp[1]->Fill(ecluster,frac);}
 357     if(GetReader()->IsPileUpFromSPDOrEMCal())        {fhPtPileUp[2]->Fill(ptcluster); fhLambda0PileUp[2]->Fill(ecluster,l0cluster); fhClusterEFracLongTimePileUp[2]->Fill(ecluster,frac);}
 358     if(GetReader()->IsPileUpFromSPDAndEMCal())       {fhPtPileUp[3]->Fill(ptcluster); fhLambda0PileUp[3]->Fill(ecluster,l0cluster); fhClusterEFracLongTimePileUp[3]->Fill(ecluster,frac);}
 359     if(GetReader()->IsPileUpFromSPDAndNotEMCal())    {fhPtPileUp[4]->Fill(ptcluster); fhLambda0PileUp[4]->Fill(ecluster,l0cluster); fhClusterEFracLongTimePileUp[4]->Fill(ecluster,frac);}
 360     if(GetReader()->IsPileUpFromEMCalAndNotSPD())    {fhPtPileUp[5]->Fill(ptcluster); fhLambda0PileUp[5]->Fill(ecluster,l0cluster); fhClusterEFracLongTimePileUp[5]->Fill(ecluster,frac);}
 361     if(GetReader()->IsPileUpFromNotSPDAndNotEMCal()) {fhPtPileUp[6]->Fill(ptcluster); fhLambda0PileUp[6]->Fill(ecluster,l0cluster); fhClusterEFracLongTimePileUp[6]->Fill(ecluster,frac);}
 362
 363     if(tmax > -25 && tmax < 25) {fhEtaPhiBC0    ->Fill(mom.Eta(),mom.Phi()); if(GetReader()->IsPileUpFromSPD()) fhEtaPhiBC0PileUpSPD    ->Fill(mom.Eta(),mom.Phi()); }
 364     else if (tmax > 25)         {fhEtaPhiBCPlus ->Fill(mom.Eta(),mom.Phi()); if(GetReader()->IsPileUpFromSPD()) fhEtaPhiBCPlusPileUpSPD ->Fill(mom.Eta(),mom.Phi()); }
 365     else if (tmax <-25)         {fhEtaPhiBCMinus->Fill(mom.Eta(),mom.Phi()); if(GetReader()->IsPileUpFromSPD()) fhEtaPhiBCMinusPileUpSPD->Fill(mom.Eta(),mom.Phi()); }
 366   }
 367
 368   //.......................................
 369   // TOF cut, BE CAREFUL WITH THIS CUT
 370   Double_t tof = calo->GetTOF()*1e9;
 371   if(tof < fTimeCutMin || tof > fTimeCutMax) return kFALSE;
 372
 373   if(GetDebug() > 2)  printf("\t Cluster %d Pass Time Cut \n",calo->GetID());
 374
 375   fhClusterCuts[3]->Fill(ecluster);
 376
 377   //.......................................
 378   if(calo->GetNCells() <= fNCellsCut && GetReader()->GetDataType() != AliCaloTrackReader::kMC) return kFALSE;
 379
 380   if(GetDebug() > 2) printf("\t Cluster %d Pass NCell Cut \n",calo->GetID());
 381
 382   fhClusterCuts[4]->Fill(ecluster);
 383
 384   if(nMaxima < fNLMCutMin || nMaxima > fNLMCutMax) return kFALSE ;
 385   if(GetDebug() > 2) printf(" \t Cluster %d pass NLM %d of out of range \n",calo->GetID(), nMaxima);
 386
 387   fhClusterCuts[5]->Fill(ecluster);
 388
 389   //.......................................
 390   //Check acceptance selection
 391   if(IsFiducialCutOn())
 392   {
 393     Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter) ;
 394     if(! in ) return kFALSE ;
 395   }
 396
 397   if(GetDebug() > 2) printf("Fiducial cut passed \n");
 398
 399   fhClusterCuts[6]->Fill(ecluster);
 400
 401   //.......................................
 402   //Skip matched clusters with tracks
 403
 404   // Fill matching residual histograms before PID cuts
 405   if(fFillTMHisto) FillTrackMatchingResidualHistograms(calo,0);
 406
 407   if(fRejectTrackMatch)
 408   {
 409     if(IsTrackMatched(calo,GetReader()->GetInputEvent()))
 410     {
 411       if(GetDebug() > 2) printf("\t Reject track-matched clusters\n");
 412       return kFALSE ;
 413     }
 414     else
 415       if(GetDebug() > 2)  printf(" Track-matching cut passed \n");
 416   }// reject matched clusters
 417
 418   fhClusterCuts[7]->Fill(ecluster);
 419
 420   if(fFillPileUpHistograms)
 421   {
 422     if(GetReader()->IsPileUpFromSPD())               {fhPtChargedPileUp[0]->Fill(ptcluster); fhLambda0ChargedPileUp[0]->Fill(ecluster,l0cluster); }
 423     if(GetReader()->IsPileUpFromEMCal())             {fhPtChargedPileUp[1]->Fill(ptcluster); fhLambda0ChargedPileUp[1]->Fill(ecluster,l0cluster); }
 424     if(GetReader()->IsPileUpFromSPDOrEMCal())        {fhPtChargedPileUp[2]->Fill(ptcluster); fhLambda0ChargedPileUp[2]->Fill(ecluster,l0cluster); }
 425     if(GetReader()->IsPileUpFromSPDAndEMCal())       {fhPtChargedPileUp[3]->Fill(ptcluster); fhLambda0ChargedPileUp[3]->Fill(ecluster,l0cluster); }
 426     if(GetReader()->IsPileUpFromSPDAndNotEMCal())    {fhPtChargedPileUp[4]->Fill(ptcluster); fhLambda0ChargedPileUp[4]->Fill(ecluster,l0cluster); }
 427     if(GetReader()->IsPileUpFromEMCalAndNotSPD())    {fhPtChargedPileUp[5]->Fill(ptcluster); fhLambda0ChargedPileUp[5]->Fill(ecluster,l0cluster); }
 428     if(GetReader()->IsPileUpFromNotSPDAndNotEMCal()) {fhPtChargedPileUp[6]->Fill(ptcluster); fhLambda0ChargedPileUp[6]->Fill(ecluster,l0cluster); }
 429   }
 430
 431   //.......................................
 432   //Check Distance to Bad channel, set bit.
 433   Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
 434   if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;
 435   if(distBad < fMinDist)
 436   {//In bad channel (PHOS cristal size 2.2x2.2 cm), EMCAL ( cell units )
 437     return kFALSE ;
 438   }
 439   else if(GetDebug() > 2) printf("\t Bad channel cut passed %4.2f > %2.2f \n",distBad, fMinDist);
 440
 441   fhClusterCuts[8]->Fill(ecluster);
 442
 443   if(GetDebug() > 0)
 444     printf("AliAnaPhoton::ClusterSelected() Current Event %d; After  selection : E %2.2f, pT %2.2f, phi %2.2f, eta %2.2f\n",
 445            GetReader()->GetEventNumber(),
 446            ecluster, ptcluster,mom.Phi()*TMath::RadToDeg(),mom.Eta());
 447
 448   //All checks passed, cluster selected
 449   return kTRUE;
 450
 451 }
 452
 453 //___________________________________________
 454 void AliAnaPhoton::FillAcceptanceHistograms()
 455 {
 456   //Fill acceptance histograms if MC data is available
 457
 458   Double_t photonY   = -100 ;
 459   Double_t photonE   = -1 ;
 460   Double_t photonPt  = -1 ;
 461   Double_t photonPhi =  100 ;
 462   Double_t photonEta = -1 ;
 463
 464   Int_t    pdg       =  0 ;
 465   Int_t    tag       =  0 ;
 466   Int_t    mcIndex   =  0 ;
 467   Bool_t   inacceptance = kFALSE;
 468
 469   if(GetReader()->ReadStack())
 470   {
 471     AliStack * stack = GetMCStack();
 472     if(stack)
 473     {
 474       for(Int_t i=0 ; i<stack->GetNtrack(); i++)
 475       {
 476         TParticle * prim = stack->Particle(i) ;
 477         pdg = prim->GetPdgCode();
 478         //printf("i %d, %s %d  %s %d \n",i, stack->Particle(i)->GetName(), stack->Particle(i)->GetPdgCode(),
 479         //                             prim->GetName(), prim->GetPdgCode());
 480
 481         if(pdg == 22)
 482         {
 483           // Get tag of this particle photon from fragmentation, decay, prompt ...
 484           tag = GetMCAnalysisUtils()->CheckOrigin(i,GetReader(), 0);
 485           if(!GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton))
 486           {
 487             //A conversion photon from a hadron, skip this kind of photon
 488             // printf("AliAnaPhoton::FillAcceptanceHistograms() - not a photon, weird!\n ");
 489             // GetMCAnalysisUtils()->PrintMCTag(tag);
 490
 491             return;
 492           }
 493
 494           //Get photon kinematics
 495           if(prim->Energy() == TMath::Abs(prim->Pz()))  continue ; //Protection against floating point exception
 496
 497           photonY   = 0.5*TMath::Log((prim->Energy()-prim->Pz())/(prim->Energy()+prim->Pz())) ;
 498           photonE   = prim->Energy() ;
 499           photonPt  = prim->Pt() ;
 500           photonPhi = TMath::RadToDeg()*prim->Phi() ;
 501           if(photonPhi < 0) photonPhi+=TMath::TwoPi();
 502           photonEta = prim->Eta() ;
 503
 504           //Check if photons hit the Calorimeter
 505           TLorentzVector lv;
 506           prim->Momentum(lv);
 507           inacceptance = kFALSE;
 508           if     (fCalorimeter == "PHOS")
 509           {
 510             if(GetPHOSGeometry() && GetCaloUtils()->IsPHOSGeoMatrixSet())
 511             {
 512               Int_t mod ;
 513               Double_t x,z ;
 514               if(GetPHOSGeometry()->ImpactOnEmc(prim,mod,z,x))
 515                 inacceptance = kTRUE;
 516               if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
 517             }
 518             else
 519             {
 520               if(GetFiducialCut()->IsInFiducialCut(lv,fCalorimeter))
 521                 inacceptance = kTRUE ;
 522               if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
 523             }
 524           }
 525           else if(fCalorimeter == "EMCAL" && GetCaloUtils()->IsEMCALGeoMatrixSet())
 526           {
 527             if(GetEMCALGeometry())
 528             {
 529               Int_t absID=0;
 530
 531               GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(prim->Eta(),prim->Phi(),absID);
 532
 533               if( absID >= 0)
 534                 inacceptance = kTRUE;
 535
 536               //                  if(GetEMCALGeometry()->Impact(phot1) && GetEMCALGeometry()->Impact(phot2))
 537               //                    inacceptance = kTRUE;
 538               if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
 539             }
 540             else
 541             {
 542               if(GetFiducialCut()->IsInFiducialCut(lv,fCalorimeter))
 543                 inacceptance = kTRUE ;
 544               if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
 545             }
 546           }       //In EMCAL
 547
 548           //Fill histograms
 549           fhYPrimMC[kmcPPhoton]->Fill(photonPt, photonY) ;
 550           if(TMath::Abs(photonY) < 1.0)
 551           {
 552             fhEPrimMC  [kmcPPhoton]->Fill(photonE ) ;
 553             fhPtPrimMC [kmcPPhoton]->Fill(photonPt) ;
 554             fhPhiPrimMC[kmcPPhoton]->Fill(photonE , photonPhi) ;
 555             fhYPrimMC  [kmcPPhoton]->Fill(photonE , photonEta) ;
 556           }
 557           if(inacceptance)
 558           {
 559             fhEPrimMCAcc  [kmcPPhoton]->Fill(photonE ) ;
 560             fhPtPrimMCAcc [kmcPPhoton]->Fill(photonPt) ;
 561             fhPhiPrimMCAcc[kmcPPhoton]->Fill(photonE , photonPhi) ;
 562             fhYPrimMCAcc  [kmcPPhoton]->Fill(photonE , photonY) ;
 563           }//Accepted
 564
 565           //Origin of photon
 566           if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt) && fhEPrimMC[kmcPPrompt])
 567           {
 568             mcIndex = kmcPPrompt;
 569           }
 570           else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation) && fhEPrimMC[kmcPFragmentation])
 571           {
 572             mcIndex = kmcPFragmentation ;
 573           }
 574           else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCISR) && fhEPrimMC[kmcPISR])
 575           {
 576             mcIndex = kmcPISR;
 577           }
 578           else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay)&& fhEPrimMC[kmcPPi0Decay])
 579           {
 580             mcIndex = kmcPPi0Decay;
 581           }
 582           else if( (GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay) ||
 583                     GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay)) && fhEPrimMC[kmcPOtherDecay])
 584           {
 585             mcIndex = kmcPOtherDecay;
 586           }
 587           else if(fhEPrimMC[kmcPOther])
 588           {
 589             mcIndex = kmcPOther;
 590           }//Other origin
 591
 592           fhYPrimMC[mcIndex]->Fill(photonPt, photonY) ;
 593           if(TMath::Abs(photonY) < 1.0)
 594           {
 595             fhEPrimMC  [mcIndex]->Fill(photonE ) ;
 596             fhPtPrimMC [mcIndex]->Fill(photonPt) ;
 597             fhPhiPrimMC[mcIndex]->Fill(photonE , photonPhi) ;
 598             fhYPrimMC  [mcIndex]->Fill(photonE , photonEta) ;
 599           }
 600           if(inacceptance)
 601           {
 602             fhEPrimMCAcc  [mcIndex]->Fill(photonE ) ;
 603             fhPtPrimMCAcc [mcIndex]->Fill(photonPt) ;
 604             fhPhiPrimMCAcc[mcIndex]->Fill(photonE , photonPhi) ;
 605             fhYPrimMCAcc  [mcIndex]->Fill(photonE , photonY) ;
 606           }//Accepted
 607
 608         }// Primary photon
 609       }//loop on primaries
 610     }//stack exists and data is MC
 611   }//read stack
 612   else if(GetReader()->ReadAODMCParticles())
 613   {
 614     TClonesArray * mcparticles = GetReader()->GetAODMCParticles(0);
 615     if(mcparticles)
 616     {
 617       Int_t nprim = mcparticles->GetEntriesFast();
 618
 619       for(Int_t i=0; i < nprim; i++)
 620       {
 621         AliAODMCParticle * prim = (AliAODMCParticle *) mcparticles->At(i);
 622
 623         pdg = prim->GetPdgCode();
 624
 625         if(pdg == 22)
 626         {
 627           // Get tag of this particle photon from fragmentation, decay, prompt ...
 628           tag = GetMCAnalysisUtils()->CheckOrigin(i,GetReader(), 0);
 629           if(!GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton))
 630           {
 631             //A conversion photon from a hadron, skip this kind of photon
 632             //            printf("AliAnaPhoton::FillAcceptanceHistograms() - not a photon, weird!\n ");
 633             //            GetMCAnalysisUtils()->PrintMCTag(tag);
 634
 635             return;
 636           }
 637
 638           //Get photon kinematics
 639           if(prim->E() == TMath::Abs(prim->Pz()))  continue ; //Protection against floating point exception
 640
 641           photonY   = 0.5*TMath::Log((prim->E()-prim->Pz())/(prim->E()+prim->Pz())) ;
 642           photonE   = prim->E() ;
 643           photonPt  = prim->Pt() ;
 644           photonPhi = prim->Phi() ;
 645           if(photonPhi < 0) photonPhi+=TMath::TwoPi();
 646           photonEta = prim->Eta() ;
 647
 648           //Check if photons hit the Calorimeter
 649           TLorentzVector lv;
 650           lv.SetPxPyPzE(prim->Px(),prim->Py(),prim->Pz(),prim->E());
 651           inacceptance = kFALSE;
 652           if     (fCalorimeter == "PHOS")
 653           {
 654             if(GetPHOSGeometry() && GetCaloUtils()->IsPHOSGeoMatrixSet())
 655             {
 656               Int_t mod ;
 657               Double_t x,z ;
 658               Double_t vtx[]={prim->Xv(),prim->Yv(),prim->Zv()};
 659               if(GetPHOSGeometry()->ImpactOnEmc(vtx, prim->Theta(),prim->Phi(),mod,z,x))
 660                 inacceptance = kTRUE;
 661               if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
 662             }
 663             else
 664             {
 665               if(GetFiducialCut()->IsInFiducialCut(lv,fCalorimeter))
 666                 inacceptance = kTRUE ;
 667               if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
 668             }
 669           }
 670           else if(fCalorimeter == "EMCAL" && GetCaloUtils()->IsEMCALGeoMatrixSet())
 671           {
 672             if(GetEMCALGeometry())
 673             {
 674               Int_t absID=0;
 675
 676               GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(prim->Eta(),prim->Phi(),absID);
 677
 678               if( absID >= 0)
 679                 inacceptance = kTRUE;
 680
 681               if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
 682             }
 683             else
 684             {
 685               if(GetFiducialCut()->IsInFiducialCut(lv,fCalorimeter))
 686                 inacceptance = kTRUE ;
 687               if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
 688             }
 689           }       //In EMCAL
 690
 691           //Fill histograms
 692
 693           fhYPrimMC[kmcPPhoton]->Fill(photonPt, photonY) ;
 694           if(TMath::Abs(photonY) < 1.0)
 695           {
 696             fhEPrimMC  [kmcPPhoton]->Fill(photonE ) ;
 697             fhPtPrimMC [kmcPPhoton]->Fill(photonPt) ;
 698             fhPhiPrimMC[kmcPPhoton]->Fill(photonE , photonPhi) ;
 699             fhYPrimMC[kmcPPhoton]->Fill(photonE , photonEta) ;
 700           }
 701
 702           if(inacceptance)
 703           {
 704             fhEPrimMCAcc[kmcPPhoton]  ->Fill(photonE ) ;
 705             fhPtPrimMCAcc[kmcPPhoton] ->Fill(photonPt) ;
 706             fhPhiPrimMCAcc[kmcPPhoton]->Fill(photonE , photonPhi) ;
 707             fhYPrimMCAcc[kmcPPhoton]  ->Fill(photonE , photonY) ;
 708           }//Accepted
 709
 710           //Origin of photon
 711           if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt) && fhEPrimMC[kmcPPrompt])
 712           {
 713             mcIndex = kmcPPrompt;
 714           }
 715           else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation) && fhEPrimMC[kmcPFragmentation])
 716           {
 717             mcIndex = kmcPFragmentation ;
 718           }
 719           else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCISR) && fhEPrimMC[kmcPISR])
 720           {
 721             mcIndex = kmcPISR;
 722           }
 723           else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay)&& fhEPrimMC[kmcPPi0Decay])
 724           {
 725             mcIndex = kmcPPi0Decay;
 726           }
 727           else if( (GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay) ||
 728                     GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay)) && fhEPrimMC[kmcPOtherDecay])
 729           {
 730             mcIndex = kmcPOtherDecay;
 731           }
 732           else if(fhEPrimMC[kmcPOther])
 733           {
 734             mcIndex = kmcPOther;
 735           }//Other origin
 736
 737           fhYPrimMC[mcIndex]->Fill(photonPt, photonY) ;
 738           if(TMath::Abs(photonY) < 1.0)
 739           {
 740             fhEPrimMC  [mcIndex]->Fill(photonE ) ;
 741             fhPtPrimMC [mcIndex]->Fill(photonPt) ;
 742             fhPhiPrimMC[mcIndex]->Fill(photonE , photonPhi) ;
 743             fhYPrimMC  [mcIndex]->Fill(photonE , photonEta) ;
 744           }
 745           if(inacceptance)
 746           {
 747             fhEPrimMCAcc  [mcIndex]->Fill(photonE ) ;
 748             fhPtPrimMCAcc [mcIndex]->Fill(photonPt) ;
 749             fhPhiPrimMCAcc[mcIndex]->Fill(photonE , photonPhi) ;
 750             fhYPrimMCAcc  [mcIndex]->Fill(photonE , photonY) ;
 751           }//Accepted
 752
 753         }// Primary photon
 754       }//loop on primaries
 755
 756     }//kmc array exists and data is MC
 757   }     // read AOD MC
 758 }
 759
 760 //___________________________________________________________________
 761 void AliAnaPhoton::FillPileUpHistogramsPerEvent(TObjArray * clusters)
 762 {
 763   // Fill some histograms per event to understand pile-up
 764   // Open the time cut in the reader to be more meaningful
 765
 766   if(!fFillPileUpHistograms) return;
 767
 768   // Loop on clusters, get the maximum energy cluster as reference
 769   Int_t nclusters = clusters->GetEntriesFast();
 770   Int_t   idMax = 0;
 771   Float_t  eMax = 0;
 772   Float_t  tMax = 0;
 773   for(Int_t iclus = 0; iclus < nclusters ; iclus++)
 774   {
 775           AliVCluster * clus =  (AliVCluster*) (clusters->At(iclus));
 776     if(clus->E() > eMax && TMath::Abs(clus->GetTOF()*1e9) < 20)
 777     {
 778       eMax  = clus->E();
 779       tMax  = clus->GetTOF()*1e9;
 780       idMax = iclus;
 781     }
 782   }
 783
 784   if(eMax < 5) return;
 785
 786   // Loop again on clusters to compare this max cluster t and the rest of the clusters, if E > 0.3
 787   Int_t n20  = 0;
 788   Int_t n40  = 0;
 789   Int_t n    = 0;
 790   Int_t nOK  = 0;
 791
 792   for(Int_t iclus = 0; iclus < nclusters ; iclus++)
 793   {
 794           AliVCluster * clus =  (AliVCluster*) (clusters->At(iclus));
 795
 796     if(clus->E() < 0.3 || iclus==idMax) continue;
 797
 798     Float_t tdiff = TMath::Abs(tMax-clus->GetTOF()*1e9);
 799     n++;
 800     if(tdiff < 20) nOK++;
 801     else
 802     {
 803       n20++;
 804       if(tdiff > 40 ) n40++;
 805     }
 806   }
 807
 808   // Check pile-up and fill histograms depending on the different cluster multiplicities
 809   if(GetReader()->IsPileUpFromSPD())
 810   {
 811     fhClusterMultSPDPileUp[0]->Fill(eMax,n  );
 812     fhClusterMultSPDPileUp[1]->Fill(eMax,nOK);
 813     fhClusterMultSPDPileUp[2]->Fill(eMax,n20);
 814     fhClusterMultSPDPileUp[3]->Fill(eMax,n40);
 815   }
 816   else
 817   {
 818     fhClusterMultNoPileUp[0]->Fill(eMax,n  );
 819     fhClusterMultNoPileUp[1]->Fill(eMax,nOK);
 820     fhClusterMultNoPileUp[2]->Fill(eMax,n20);
 821     fhClusterMultNoPileUp[3]->Fill(eMax,n40);
 822   }
 823
 824 }
 825
 826
 827 //_________________________________________________________________________________________________
 828 void AliAnaPhoton::FillPileUpHistograms(const Float_t energy, const Float_t pt, const Float_t time)
 829 {
 830   // Fill some histograms to understand pile-up
 831   if(!fFillPileUpHistograms) return;
 832
 833   //printf("E %f, time %f\n",energy,time);
 834   AliVEvent * event = GetReader()->GetInputEvent();
 835
 836   if(GetReader()->IsPileUpFromSPD())               fhPtPhotonPileUp[0]->Fill(pt);
 837   if(GetReader()->IsPileUpFromEMCal())             fhPtPhotonPileUp[1]->Fill(pt);
 838   if(GetReader()->IsPileUpFromSPDOrEMCal())        fhPtPhotonPileUp[2]->Fill(pt);
 839   if(GetReader()->IsPileUpFromSPDAndEMCal())       fhPtPhotonPileUp[3]->Fill(pt);
 840   if(GetReader()->IsPileUpFromSPDAndNotEMCal())    fhPtPhotonPileUp[4]->Fill(pt);
 841   if(GetReader()->IsPileUpFromEMCalAndNotSPD())    fhPtPhotonPileUp[5]->Fill(pt);
 842   if(GetReader()->IsPileUpFromNotSPDAndNotEMCal()) fhPtPhotonPileUp[6]->Fill(pt);
 843
 844   fhTimeENoCut->Fill(energy,time);
 845   if(GetReader()->IsPileUpFromSPD())     fhTimeESPD     ->Fill(energy,time);
 846   if(event->IsPileupFromSPDInMultBins()) fhTimeESPDMulti->Fill(energy,time);
 847
 848   if(energy < 8) return; // Fill time figures for high energy clusters not too close to trigger threshold
 849
 850   AliESDEvent* esdEv = dynamic_cast<AliESDEvent*> (event);
 851   AliAODEvent* aodEv = dynamic_cast<AliAODEvent*> (event);
 852
 853   // N pile up vertices
 854   Int_t nVerticesSPD    = -1;
 855   Int_t nVerticesTracks = -1;
 856
 857   if      (esdEv)
 858   {
 859     nVerticesSPD    = esdEv->GetNumberOfPileupVerticesSPD();
 860     nVerticesTracks = esdEv->GetNumberOfPileupVerticesTracks();
 861
 862   }//ESD
 863   else if (aodEv)
 864   {
 865     nVerticesSPD    = aodEv->GetNumberOfPileupVerticesSPD();
 866     nVerticesTracks = aodEv->GetNumberOfPileupVerticesTracks();
 867   }//AOD
 868
 869   fhTimeNPileUpVertSPD  ->Fill(time,nVerticesSPD);
 870   fhTimeNPileUpVertTrack->Fill(time,nVerticesTracks);
 871
 872   //printf("Is SPD %d, Is SPD Multi %d, n spd %d, n track %d\n",
 873   //       GetReader()->IsPileUpFromSPD(),event->IsPileupFromSPDInMultBins(),nVerticesSPD,nVerticesTracks);
 874
 875   Int_t ncont = -1;
 876   Float_t z1 = -1, z2 = -1;
 877   Float_t diamZ = -1;
 878   for(Int_t iVert=0; iVert<nVerticesSPD;iVert++)
 879   {
 880     if      (esdEv)
 881     {
 882       const AliESDVertex* pv=esdEv->GetPileupVertexSPD(iVert);
 883       ncont=pv->GetNContributors();
 884       z1 = esdEv->GetPrimaryVertexSPD()->GetZ();
 885       z2 = pv->GetZ();
 886       diamZ = esdEv->GetDiamondZ();
 887     }//ESD
 888     else if (aodEv)
 889     {
 890       AliAODVertex *pv=aodEv->GetVertex(iVert);
 891       if(pv->GetType()!=AliAODVertex::kPileupSPD) continue;
 892       ncont=pv->GetNContributors();
 893       z1=aodEv->GetPrimaryVertexSPD()->GetZ();
 894       z2=pv->GetZ();
 895       diamZ = aodEv->GetDiamondZ();
 896     }// AOD
 897
 898     Double_t distZ  = TMath::Abs(z2-z1);
 899     diamZ  = TMath::Abs(z2-diamZ);
 900
 901     fhTimeNPileUpVertContributors  ->Fill(time,ncont);
 902     fhTimePileUpMainVertexZDistance->Fill(time,distZ);
 903     fhTimePileUpMainVertexZDiamond ->Fill(time,diamZ);
 904
 905   }// loop
 906 }
 907
 908 //____________________________________________________________________________________
 909 void  AliAnaPhoton::FillShowerShapeHistograms(AliVCluster* cluster, const Int_t mcTag)
 910 {
 911     //Fill cluster Shower Shape histograms
 912
 913   if(!fFillSSHistograms || GetMixedEvent()) return;
 914
 915   Float_t energy  = cluster->E();
 916   Int_t   ncells  = cluster->GetNCells();
 917   Float_t lambda0 = cluster->GetM02();
 918   Float_t lambda1 = cluster->GetM20();
 919   Float_t disp    = cluster->GetDispersion()*cluster->GetDispersion();
 920
 921   TLorentzVector mom;
 922   if(GetReader()->GetDataType() != AliCaloTrackReader::kMC)
 923   {
 924     cluster->GetMomentum(mom,GetVertex(0)) ;
 925   }//Assume that come from vertex in straight line
 926   else
 927   {
 928     Double_t vertex[]={0,0,0};
 929     cluster->GetMomentum(mom,vertex) ;
 930   }
 931
 932   Float_t eta = mom.Eta();
 933   Float_t phi = mom.Phi();
 934   if(phi < 0) phi+=TMath::TwoPi();
 935
 936   fhLam0E ->Fill(energy,lambda0);
 937   fhLam1E ->Fill(energy,lambda1);
 938   fhDispE ->Fill(energy,disp);
 939
 940   if(fCalorimeter == "EMCAL" && GetModuleNumber(cluster) > 5)
 941   {
 942     fhLam0ETRD->Fill(energy,lambda0);
 943     fhLam1ETRD->Fill(energy,lambda1);
 944     fhDispETRD->Fill(energy,disp);
 945   }
 946
 947   Float_t l0   = 0., l1   = 0.;
 948   Float_t dispp= 0., dEta = 0., dPhi    = 0.;
 949   Float_t sEta = 0., sPhi = 0., sEtaPhi = 0.;
 950   if(fCalorimeter == "EMCAL" && !fFillOnlySimpleSSHisto)
 951   {
 952     GetCaloUtils()->GetEMCALRecoUtils()->RecalculateClusterShowerShapeParameters(GetEMCALGeometry(), GetReader()->GetInputEvent()->GetEMCALCells(), cluster,
 953                                                                                  l0, l1, dispp, dEta, dPhi, sEta, sPhi, sEtaPhi);
 954     //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",
 955     //       l0, l1, dispp, dEta, dPhi, sEta, sPhi, sEtaPhi );
 956     //printf("AliAnaPhoton::FillShowerShapeHistogram - dispersion %f, dispersion eta+phi %f \n",
 957     //       disp, dPhi+dEta );
 958     fhDispEtaE        -> Fill(energy,dEta);
 959     fhDispPhiE        -> Fill(energy,dPhi);
 960     fhSumEtaE         -> Fill(energy,sEta);
 961     fhSumPhiE         -> Fill(energy,sPhi);
 962     fhSumEtaPhiE      -> Fill(energy,sEtaPhi);
 963     fhDispEtaPhiDiffE -> Fill(energy,dPhi-dEta);
 964     if(dEta+dPhi>0)fhSphericityE     -> Fill(energy,(dPhi-dEta)/(dEta+dPhi));
 965     if(dEta+sEta>0)fhDispSumEtaDiffE -> Fill(energy,(dEta-sEta)/((dEta+sEta)/2.));
 966     if(dPhi+sPhi>0)fhDispSumPhiDiffE -> Fill(energy,(dPhi-sPhi)/((dPhi+sPhi)/2.));
 967
 968     Int_t ebin = -1;
 969     if      (energy < 2 ) ebin = 0;
 970     else if (energy < 4 ) ebin = 1;
 971     else if (energy < 6 ) ebin = 2;
 972     else if (energy < 10) ebin = 3;
 973     else if (energy < 15) ebin = 4;
 974     else if (energy < 20) ebin = 5;
 975     else                  ebin = 6;
 976
 977     fhDispEtaDispPhi[ebin]->Fill(dEta   ,dPhi);
 978     fhLambda0DispEta[ebin]->Fill(lambda0,dEta);
 979     fhLambda0DispPhi[ebin]->Fill(lambda0,dPhi);
 980
 981   }
 982
 983   // if track-matching was of, check effect of track-matching residual cut
 984
 985   if(!fRejectTrackMatch)
 986   {
 987     Float_t dZ  = cluster->GetTrackDz();
 988     Float_t dR  = cluster->GetTrackDx();
 989     if(cluster->IsEMCAL() && GetCaloUtils()->IsRecalculationOfClusterTrackMatchingOn())
 990     {
 991       dR = 2000., dZ = 2000.;
 992       GetCaloUtils()->GetEMCALRecoUtils()->GetMatchedResiduals(cluster->GetID(),dZ,dR);
 993     }
 994
 995     if(TMath::Abs(dZ) < 0.05 && TMath::Abs(dR) < 0.05)
 996     {
 997       fhLam0ETM ->Fill(energy,lambda0);
 998       fhLam1ETM ->Fill(energy,lambda1);
 999       fhDispETM ->Fill(energy,disp);
1000
1001       if(fCalorimeter == "EMCAL" && GetModuleNumber(cluster) > 5)
1002       {
1003         fhLam0ETMTRD->Fill(energy,lambda0);
1004         fhLam1ETMTRD->Fill(energy,lambda1);
1005         fhDispETMTRD->Fill(energy,disp);
1006       }
1007     }
1008   }// if track-matching was of, check effect of matching residual cut
1009
1010
1011   if(!fFillOnlySimpleSSHisto){
1012     if(energy < 2)
1013     {
1014       fhNCellsLam0LowE ->Fill(ncells,lambda0);
1015       fhNCellsLam1LowE ->Fill(ncells,lambda1);
1016       fhNCellsDispLowE ->Fill(ncells,disp);
1017
1018       fhLam1Lam0LowE  ->Fill(lambda1,lambda0);
1019       fhLam0DispLowE  ->Fill(lambda0,disp);
1020       fhDispLam1LowE  ->Fill(disp,lambda1);
1021       fhEtaLam0LowE   ->Fill(eta,lambda0);
1022       fhPhiLam0LowE   ->Fill(phi,lambda0);
1023     }
1024     else
1025     {
1026       fhNCellsLam0HighE ->Fill(ncells,lambda0);
1027       fhNCellsLam1HighE ->Fill(ncells,lambda1);
1028       fhNCellsDispHighE ->Fill(ncells,disp);
1029
1030       fhLam1Lam0HighE  ->Fill(lambda1,lambda0);
1031       fhLam0DispHighE  ->Fill(lambda0,disp);
1032       fhDispLam1HighE  ->Fill(disp,lambda1);
1033       fhEtaLam0HighE   ->Fill(eta, lambda0);
1034       fhPhiLam0HighE   ->Fill(phi, lambda0);
1035     }
1036   }
1037
1038   if(IsDataMC())
1039   {
1040     AliVCaloCells* cells = 0;
1041     if(fCalorimeter == "EMCAL") cells = GetEMCALCells();
1042     else                        cells = GetPHOSCells();
1043
1044     //Fill histograms to check shape of embedded clusters
1045     Float_t fraction = 0;
1046     if(GetReader()->IsEmbeddedClusterSelectionOn())
1047     {//Only working for EMCAL
1048       Float_t clusterE = 0; // recalculate in case corrections applied.
1049       Float_t cellE    = 0;
1050       for(Int_t icell  = 0; icell < cluster->GetNCells(); icell++)
1051       {
1052         cellE    = cells->GetCellAmplitude(cluster->GetCellAbsId(icell));
1053         clusterE+=cellE;
1054         fraction+=cellE*cluster->GetCellAmplitudeFraction(icell);
1055       }
1056
1057       //Fraction of total energy due to the embedded signal
1058       fraction/=clusterE;
1059
1060       if(GetDebug() > 1 )
1061         printf("AliAnaPhoton::FillShowerShapeHistogram() - Energy fraction of embedded signal %2.3f, Energy %2.3f\n",fraction, clusterE);
1062
1063       fhEmbeddedSignalFractionEnergy->Fill(clusterE,fraction);
1064
1065     }  // embedded fraction
1066
1067     // Get the fraction of the cluster energy that carries the cell with highest energy
1068     Int_t   absID           =-1 ;
1069     Float_t maxCellFraction = 0.;
1070
1071     absID = GetCaloUtils()->GetMaxEnergyCell(cells, cluster,maxCellFraction);
1072
1073     // Check the origin and fill histograms
1074
1075     Int_t mcIndex = -1;
1076
1077     if( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton)     &&
1078        !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCConversion) &&
1079        !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0)        &&
1080        !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEta))
1081     {
1082       mcIndex = kmcssPhoton ;
1083
1084       if(!GetReader()->IsEmbeddedClusterSelectionOn())
1085       {
1086         //Check particle overlaps in cluster
1087
1088         // Compare the primary depositing more energy with the rest,
1089         // if no photon/electron as comon ancestor (conversions), count as other particle
1090         Int_t ancPDG = 0, ancStatus = -1;
1091         TLorentzVector momentum; TVector3 prodVertex;
1092         Int_t ancLabel = 0;
1093         Int_t noverlaps = 1;
1094         for (UInt_t ilab = 0; ilab < cluster->GetNLabels(); ilab++ )
1095         {
1096           ancLabel = GetMCAnalysisUtils()->CheckCommonAncestor(cluster->GetLabels()[0],cluster->GetLabels()[ilab],
1097                                                                GetReader(),ancPDG,ancStatus,momentum,prodVertex);
1098           if(ancPDG!=22 && TMath::Abs(ancPDG)!=11) noverlaps++;
1099         }
1100         //printf("N overlaps %d \n",noverlaps);
1101
1102         if(noverlaps == 1)
1103         {
1104           fhMCPhotonELambda0NoOverlap  ->Fill(energy, lambda0);
1105         }
1106         else if(noverlaps == 2)
1107         {
1108           fhMCPhotonELambda0TwoOverlap ->Fill(energy, lambda0);
1109         }
1110         else if(noverlaps > 2)
1111         {
1112           fhMCPhotonELambda0NOverlap   ->Fill(energy, lambda0);
1113         }
1114         else
1115         {
1116           printf("AliAnaPhoton::FillShowerShapeHistogram() - n overlaps = %d!!", noverlaps);
1117         }
1118       }//No embedding
1119
1120       //Fill histograms to check shape of embedded clusters
1121       if(GetReader()->IsEmbeddedClusterSelectionOn())
1122       {
1123         if     (fraction > 0.9)
1124         {
1125           fhEmbedPhotonELambda0FullSignal   ->Fill(energy, lambda0);
1126         }
1127         else if(fraction > 0.5)
1128         {
1129           fhEmbedPhotonELambda0MostlySignal ->Fill(energy, lambda0);
1130         }
1131         else if(fraction > 0.1)
1132         {
1133           fhEmbedPhotonELambda0MostlyBkg    ->Fill(energy, lambda0);
1134         }
1135         else
1136         {
1137           fhEmbedPhotonELambda0FullBkg      ->Fill(energy, lambda0);
1138         }
1139       } // embedded
1140
1141     }//photon   no conversion
1142     else if  ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCElectron))
1143     {
1144       mcIndex = kmcssElectron ;
1145     }//electron
1146     else if  ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton) &&
1147                GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCConversion) )
1148     {
1149       mcIndex = kmcssConversion ;
1150     }//conversion photon
1151     else if  ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0)  )
1152     {
1153       mcIndex = kmcssPi0 ;
1154
1155       //Fill histograms to check shape of embedded clusters
1156       if(GetReader()->IsEmbeddedClusterSelectionOn())
1157       {
1158         if     (fraction > 0.9)
1159         {
1160           fhEmbedPi0ELambda0FullSignal   ->Fill(energy, lambda0);
1161         }
1162         else if(fraction > 0.5)
1163         {
1164           fhEmbedPi0ELambda0MostlySignal ->Fill(energy, lambda0);
1165         }
1166         else if(fraction > 0.1)
1167         {
1168           fhEmbedPi0ELambda0MostlyBkg    ->Fill(energy, lambda0);
1169         }
1170         else
1171         {
1172           fhEmbedPi0ELambda0FullBkg      ->Fill(energy, lambda0);
1173         }
1174       } // embedded
1175
1176     }//pi0
1177     else if  ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEta)  )
1178     {
1179       mcIndex = kmcssEta ;
1180     }//eta
1181     else
1182     {
1183       mcIndex = kmcssOther ;
1184     }//other particles
1185
1186     fhMCELambda0           [mcIndex]->Fill(energy, lambda0);
1187     fhMCELambda1           [mcIndex]->Fill(energy, lambda1);
1188     fhMCEDispersion        [mcIndex]->Fill(energy, disp);
1189     fhMCNCellsE            [mcIndex]->Fill(energy, ncells);
1190     fhMCMaxCellDiffClusterE[mcIndex]->Fill(energy, maxCellFraction);
1191
1192     if(!fFillOnlySimpleSSHisto)
1193     {
1194       if     (energy < 2.)
1195       {
1196         fhMCLambda0vsClusterMaxCellDiffE0[mcIndex]->Fill(lambda0, maxCellFraction);
1197         fhMCNCellsvsClusterMaxCellDiffE0 [mcIndex]->Fill(ncells,  maxCellFraction);
1198       }
1199       else if(energy < 6.)
1200       {
1201         fhMCLambda0vsClusterMaxCellDiffE2[mcIndex]->Fill(lambda0, maxCellFraction);
1202         fhMCNCellsvsClusterMaxCellDiffE2 [mcIndex]->Fill(ncells,  maxCellFraction);
1203       }
1204       else
1205       {
1206         fhMCLambda0vsClusterMaxCellDiffE6[mcIndex]->Fill(lambda0, maxCellFraction);
1207         fhMCNCellsvsClusterMaxCellDiffE6 [mcIndex]->Fill(ncells,  maxCellFraction);
1208       }
1209
1210       if(fCalorimeter == "EMCAL")
1211       {
1212         fhMCEDispEta        [mcIndex]-> Fill(energy,dEta);
1213         fhMCEDispPhi        [mcIndex]-> Fill(energy,dPhi);
1214         fhMCESumEtaPhi      [mcIndex]-> Fill(energy,sEtaPhi);
1215         fhMCEDispEtaPhiDiff [mcIndex]-> Fill(energy,dPhi-dEta);
1216         if(dEta+dPhi>0)fhMCESphericity[mcIndex]-> Fill(energy,(dPhi-dEta)/(dEta+dPhi));
1217
1218         Int_t ebin = -1;
1219         if      (energy < 2 ) ebin = 0;
1220         else if (energy < 4 ) ebin = 1;
1221         else if (energy < 6 ) ebin = 2;
1222         else if (energy < 10) ebin = 3;
1223         else if (energy < 15) ebin = 4;
1224         else if (energy < 20) ebin = 5;
1225         else                  ebin = 6;
1226
1227         fhMCDispEtaDispPhi[ebin][mcIndex]->Fill(dEta   ,dPhi);
1228         fhMCLambda0DispEta[ebin][mcIndex]->Fill(lambda0,dEta);
1229         fhMCLambda0DispPhi[ebin][mcIndex]->Fill(lambda0,dPhi);
1230       }
1231     }
1232   }//MC data
1233
1234 }
1235
1236 //__________________________________________________________________________
1237 void AliAnaPhoton::FillTrackMatchingResidualHistograms(AliVCluster* cluster,
1238                                                        const Int_t cut)
1239 {
1240   // If selected, fill histograms with residuals of matched clusters, help to define track matching cut
1241   // Residual filled for different cuts 0 (No cut), after 1 PID cut
1242
1243   Float_t dZ = cluster->GetTrackDz();
1244   Float_t dR = cluster->GetTrackDx();
1245
1246   if(cluster->IsEMCAL() && GetCaloUtils()->IsRecalculationOfClusterTrackMatchingOn())
1247   {
1248     dR = 2000., dZ = 2000.;
1249     GetCaloUtils()->GetEMCALRecoUtils()->GetMatchedResiduals(cluster->GetID(),dZ,dR);
1250   }
1251
1252   if(fhTrackMatchedDEta[cut] && TMath::Abs(dR) < 999)
1253   {
1254     fhTrackMatchedDEta[cut]->Fill(cluster->E(),dZ);
1255     fhTrackMatchedDPhi[cut]->Fill(cluster->E(),dR);
1256
1257     if(cluster->E() > 0.5) fhTrackMatchedDEtaDPhi[cut]->Fill(dZ,dR);
1258
1259     Int_t nSMod = GetModuleNumber(cluster);
1260
1261     if(fCalorimeter=="EMCAL" &&  nSMod > 5)
1262     {
1263       fhTrackMatchedDEtaTRD[cut]->Fill(cluster->E(),dZ);
1264       fhTrackMatchedDPhiTRD[cut]->Fill(cluster->E(),dR);
1265     }
1266
1267     // Check dEdx and E/p of matched clusters
1268
1269     if(TMath::Abs(dZ) < 0.05 && TMath::Abs(dR) < 0.05)
1270     {
1271
1272       AliVTrack *track = GetCaloUtils()->GetMatchedTrack(cluster, GetReader()->GetInputEvent());
1273
1274       if(track)
1275       {
1276
1277         Float_t dEdx   = track->GetTPCsignal();
1278         Float_t eOverp = cluster->E()/track->P();
1279
1280         fhdEdx[cut]  ->Fill(cluster->E(), dEdx);
1281         fhEOverP[cut]->Fill(cluster->E(), eOverp);
1282
1283         if(fCalorimeter=="EMCAL" && nSMod > 5)
1284           fhEOverPTRD[cut]->Fill(cluster->E(), eOverp);
1285
1286
1287       }
1288       else
1289           printf("AliAnaPhoton::FillTrackMatchingResidualHistograms() - Residual OK but (dR, dZ)= (%2.4f,%2.4f) no track associated WHAT? \n", dR,dZ);
1290
1291
1292
1293       if(IsDataMC())
1294       {
1295
1296         Int_t tag = GetMCAnalysisUtils()->CheckOrigin(cluster->GetLabels(),cluster->GetNLabels(),GetReader(), 0);
1297
1298         if  ( !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion)  )
1299         {
1300           if       ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0)      ||
1301                      GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta)       ) fhTrackMatchedMCParticle[cut]->Fill(cluster->E(), 2.5 );
1302           else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton)    ) fhTrackMatchedMCParticle[cut]->Fill(cluster->E(), 0.5 );
1303           else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron)  ) fhTrackMatchedMCParticle[cut]->Fill(cluster->E(), 1.5 );
1304           else                                                                                 fhTrackMatchedMCParticle[cut]->Fill(cluster->E(), 3.5 );
1305
1306           // Check if several particles contributed to cluster and discard overlapped mesons
1307           if(!GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) ||
1308              !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta))
1309           {
1310             if(cluster->GetNLabels()==1)
1311             {
1312               fhTrackMatchedDEtaMCNoOverlap[cut]->Fill(cluster->E(),dZ);
1313               fhTrackMatchedDPhiMCNoOverlap[cut]->Fill(cluster->E(),dR);
1314             }
1315             else
1316             {
1317               fhTrackMatchedDEtaMCOverlap[cut]->Fill(cluster->E(),dZ);
1318               fhTrackMatchedDPhiMCOverlap[cut]->Fill(cluster->E(),dR);
1319             }
1320
1321           }// Check overlaps
1322
1323         }
1324         else
1325         {
1326           if       ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0)      ||
1327                      GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta)       ) fhTrackMatchedMCParticle[cut]->Fill(cluster->E(), 6.5 );
1328           else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton)    ) fhTrackMatchedMCParticle[cut]->Fill(cluster->E(), 4.5 );
1329           else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron)  ) fhTrackMatchedMCParticle[cut]->Fill(cluster->E(), 5.5 );
1330           else                                                                                 fhTrackMatchedMCParticle[cut]->Fill(cluster->E(), 7.5 );
1331
1332           // Check if several particles contributed to cluster
1333           if(!GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) ||
1334              !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta))
1335           {
1336             fhTrackMatchedDEtaMCConversion[cut]->Fill(cluster->E(),dZ);
1337             fhTrackMatchedDPhiMCConversion[cut]->Fill(cluster->E(),dR);
1338
1339           }// Check overlaps
1340
1341         }
1342
1343       } // MC
1344
1345     } // residuals window
1346
1347   } // Small residual
1348
1349 }
1350
1351 //___________________________________________
1352 TObjString *  AliAnaPhoton::GetAnalysisCuts()
1353 {
1354   //Save parameters used for analysis
1355   TString parList ; //this will be list of parameters used for this analysis.
1356   const Int_t buffersize = 255;
1357   char onePar[buffersize] ;
1358
1359   snprintf(onePar,buffersize,"--- AliAnaPhoton ---\n") ;
1360   parList+=onePar ;
1361   snprintf(onePar,buffersize,"Calorimeter: %s\n",fCalorimeter.Data()) ;
1362   parList+=onePar ;
1363   snprintf(onePar,buffersize,"fMinDist =%2.2f (Minimal distance to bad channel to accept cluster) \n",fMinDist) ;
1364   parList+=onePar ;
1365   snprintf(onePar,buffersize,"fMinDist2=%2.2f (Cuts on Minimal distance to study acceptance evaluation) \n",fMinDist2) ;
1366   parList+=onePar ;
1367   snprintf(onePar,buffersize,"fMinDist3=%2.2f (One more cut on distance used for acceptance-efficiency study) \n",fMinDist3) ;
1368   parList+=onePar ;
1369   snprintf(onePar,buffersize,"fRejectTrackMatch: %d\n",fRejectTrackMatch) ;
1370   parList+=onePar ;
1371
1372   //Get parameters set in base class.
1373   parList += GetBaseParametersList() ;
1374
1375   //Get parameters set in PID class.
1376   parList += GetCaloPID()->GetPIDParametersList() ;
1377
1378   //Get parameters set in FiducialCut class (not available yet)
1379   //parlist += GetFidCut()->GetFidCutParametersList()
1380
1381   return new TObjString(parList) ;
1382 }
1383
1384 //________________________________________________________________________
1385 TList *  AliAnaPhoton::GetCreateOutputObjects()
1386 {
1387   // Create histograms to be saved in output file and
1388   // store them in outputContainer
1389   TList * outputContainer = new TList() ;
1390   outputContainer->SetName("PhotonHistos") ;
1391
1392   Int_t nptbins  = GetHistogramRanges()->GetHistoPtBins();  Float_t ptmax  = GetHistogramRanges()->GetHistoPtMax();  Float_t ptmin  = GetHistogramRanges()->GetHistoPtMin();
1393   Int_t nphibins = GetHistogramRanges()->GetHistoPhiBins(); Float_t phimax = GetHistogramRanges()->GetHistoPhiMax(); Float_t phimin = GetHistogramRanges()->GetHistoPhiMin();
1394   Int_t netabins = GetHistogramRanges()->GetHistoEtaBins(); Float_t etamax = GetHistogramRanges()->GetHistoEtaMax(); Float_t etamin = GetHistogramRanges()->GetHistoEtaMin();
1395   Int_t ssbins   = GetHistogramRanges()->GetHistoShowerShapeBins();  Float_t ssmax   = GetHistogramRanges()->GetHistoShowerShapeMax();  Float_t ssmin   = GetHistogramRanges()->GetHistoShowerShapeMin();
1396   Int_t nbins    = GetHistogramRanges()->GetHistoNClusterCellBins(); Int_t   nmax    = GetHistogramRanges()->GetHistoNClusterCellMax(); Int_t   nmin    = GetHistogramRanges()->GetHistoNClusterCellMin();
1397   Int_t ntimebins= GetHistogramRanges()->GetHistoTimeBins();         Float_t timemax = GetHistogramRanges()->GetHistoTimeMax();         Float_t timemin = GetHistogramRanges()->GetHistoTimeMin();
1398
1399   Int_t   nresetabins = GetHistogramRanges()->GetHistoTrackResidualEtaBins();
1400   Float_t resetamax   = GetHistogramRanges()->GetHistoTrackResidualEtaMax();
1401   Float_t resetamin   = GetHistogramRanges()->GetHistoTrackResidualEtaMin();
1402   Int_t   nresphibins = GetHistogramRanges()->GetHistoTrackResidualPhiBins();
1403   Float_t resphimax   = GetHistogramRanges()->GetHistoTrackResidualPhiMax();
1404   Float_t resphimin   = GetHistogramRanges()->GetHistoTrackResidualPhiMin();
1405
1406   Int_t   ndedxbins   = GetHistogramRanges()->GetHistodEdxBins();
1407   Float_t dedxmax     = GetHistogramRanges()->GetHistodEdxMax();
1408   Float_t dedxmin     = GetHistogramRanges()->GetHistodEdxMin();
1409   Int_t   nPoverEbins = GetHistogramRanges()->GetHistoPOverEBins();
1410   Float_t pOverEmax   = GetHistogramRanges()->GetHistoPOverEMax();
1411   Float_t pOverEmin   = GetHistogramRanges()->GetHistoPOverEMin();
1412
1413   Int_t bin[] = {0,2,4,6,10,15,20,100}; // energy bins for SS studies
1414
1415   TString cut[] = {"Open","Reader","E","Time","NCells","NLM","Fidutial","Matching","Bad","PID"};
1416   for (Int_t i = 0; i < 10 ;  i++)
1417   {
1418     fhClusterCuts[i] = new TH1F(Form("hCut_%d_%s", i, cut[i].Data()),
1419                                 Form("Number of clusters that pass cuts <= %d, %s", i, cut[i].Data()),
1420                                 nptbins,ptmin,ptmax);
1421     fhClusterCuts[i]->SetYTitle("dN/dE ");
1422     fhClusterCuts[i]->SetXTitle("E (GeV)");
1423     outputContainer->Add(fhClusterCuts[i]) ;
1424   }
1425
1426   fhNCellsE  = new TH2F ("hNCellsE","# of cells in cluster vs E of clusters", nptbins,ptmin,ptmax, nbins,nmin,nmax);
1427   fhNCellsE->SetXTitle("E (GeV)");
1428   fhNCellsE->SetYTitle("# of cells in cluster");
1429   outputContainer->Add(fhNCellsE);
1430
1431   fhCellsE  = new TH2F ("hCellsE","energy of cells in cluster vs E of clusters", nptbins,ptmin,ptmax, nptbins*2,ptmin,ptmax);
1432   fhCellsE->SetXTitle("E_{cluster} (GeV)");
1433   fhCellsE->SetYTitle("E_{cell} (GeV)");
1434   outputContainer->Add(fhCellsE);
1435
1436   fhTimeE  = new TH2F ("hTimeE","time of cluster vs E of clusters", nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
1437   fhTimeE->SetXTitle("E (GeV)");
1438   fhTimeE->SetYTitle("time (ns)");
1439   outputContainer->Add(fhTimeE);
1440
1441   fhMaxCellDiffClusterE  = new TH2F ("hMaxCellDiffClusterE","energy vs difference of cluster energy - max cell energy / cluster energy, good clusters",
1442                                     nptbins,ptmin,ptmax, 500,0,1.);
1443   fhMaxCellDiffClusterE->SetXTitle("E_{cluster} (GeV) ");
1444   fhMaxCellDiffClusterE->SetYTitle("(E_{cluster} - E_{cell max})/ E_{cluster}");
1445   outputContainer->Add(fhMaxCellDiffClusterE);
1446
1447   fhEPhoton  = new TH1F("hEPhoton","Number of #gamma over calorimeter vs energy",nptbins,ptmin,ptmax);
1448   fhEPhoton->SetYTitle("N");
1449   fhEPhoton->SetXTitle("E_{#gamma}(GeV)");
1450   outputContainer->Add(fhEPhoton) ;
1451
1452   fhPtPhoton  = new TH1F("hPtPhoton","Number of #gamma over calorimeter vs p_{T}",nptbins,ptmin,ptmax);
1453   fhPtPhoton->SetYTitle("N");
1454   fhPtPhoton->SetXTitle("p_{T #gamma}(GeV/c)");
1455   outputContainer->Add(fhPtPhoton) ;
1456
1457   fhPhiPhoton  = new TH2F
1458     ("hPhiPhoton","#phi_{#gamma} vs p_{T}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
1459   fhPhiPhoton->SetYTitle("#phi (rad)");
1460   fhPhiPhoton->SetXTitle("p_{T #gamma} (GeV/c)");
1461   outputContainer->Add(fhPhiPhoton) ;
1462
1463   fhEtaPhoton  = new TH2F
1464     ("hEtaPhoton","#eta_{#gamma} vs p_{T}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
1465   fhEtaPhoton->SetYTitle("#eta");
1466   fhEtaPhoton->SetXTitle("p_{T #gamma} (GeV/c)");
1467   outputContainer->Add(fhEtaPhoton) ;
1468
1469   fhEtaPhiPhoton  = new TH2F
1470   ("hEtaPhiPhoton","#eta vs #phi",netabins,etamin,etamax,nphibins,phimin,phimax);
1471   fhEtaPhiPhoton->SetYTitle("#phi (rad)");
1472   fhEtaPhiPhoton->SetXTitle("#eta");
1473   outputContainer->Add(fhEtaPhiPhoton) ;
1474   if(GetMinPt() < 0.5)
1475   {
1476     fhEtaPhi05Photon  = new TH2F
1477     ("hEtaPhi05Photon","#eta vs #phi, E > 0.5",netabins,etamin,etamax,nphibins,phimin,phimax);
1478     fhEtaPhi05Photon->SetYTitle("#phi (rad)");
1479     fhEtaPhi05Photon->SetXTitle("#eta");
1480     outputContainer->Add(fhEtaPhi05Photon) ;
1481   }
1482
1483
1484   fhNLocMax = new TH2F("hNLocMax","Number of local maxima in cluster",
1485                        nptbins,ptmin,ptmax,10,0,10);
1486   fhNLocMax ->SetYTitle("N maxima");
1487   fhNLocMax ->SetXTitle("E (GeV)");
1488   outputContainer->Add(fhNLocMax) ;
1489
1490   //Shower shape
1491   if(fFillSSHistograms)
1492   {
1493     fhLam0E  = new TH2F ("hLam0E","#lambda_{0}^{2} vs E", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1494     fhLam0E->SetYTitle("#lambda_{0}^{2}");
1495     fhLam0E->SetXTitle("E (GeV)");
1496     outputContainer->Add(fhLam0E);
1497
1498     fhLam1E  = new TH2F ("hLam1E","#lambda_{1}^{2} vs E", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1499     fhLam1E->SetYTitle("#lambda_{1}^{2}");
1500     fhLam1E->SetXTitle("E (GeV)");
1501     outputContainer->Add(fhLam1E);
1502
1503     fhDispE  = new TH2F ("hDispE"," dispersion^{2} vs E", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1504     fhDispE->SetYTitle("D^{2}");
1505     fhDispE->SetXTitle("E (GeV) ");
1506     outputContainer->Add(fhDispE);
1507
1508     if(!fRejectTrackMatch)
1509     {
1510       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);
1511       fhLam0ETM->SetYTitle("#lambda_{0}^{2}");
1512       fhLam0ETM->SetXTitle("E (GeV)");
1513       outputContainer->Add(fhLam0ETM);
1514
1515       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);
1516       fhLam1ETM->SetYTitle("#lambda_{1}^{2}");
1517       fhLam1ETM->SetXTitle("E (GeV)");
1518       outputContainer->Add(fhLam1ETM);
1519
1520       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);
1521       fhDispETM->SetYTitle("D^{2}");
1522       fhDispETM->SetXTitle("E (GeV) ");
1523       outputContainer->Add(fhDispETM);
1524     }
1525
1526     if(fCalorimeter == "EMCAL")
1527     {
1528       fhLam0ETRD  = new TH2F ("hLam0ETRD","#lambda_{0}^{2} vs E, EMCAL SM covered by TRD", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1529       fhLam0ETRD->SetYTitle("#lambda_{0}^{2}");
1530       fhLam0ETRD->SetXTitle("E (GeV)");
1531       outputContainer->Add(fhLam0ETRD);
1532
1533       fhLam1ETRD  = new TH2F ("hLam1ETRD","#lambda_{1}^{2} vs E, EMCAL SM covered by TRD", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1534       fhLam1ETRD->SetYTitle("#lambda_{1}^{2}");
1535       fhLam1ETRD->SetXTitle("E (GeV)");
1536       outputContainer->Add(fhLam1ETRD);
1537
1538       fhDispETRD  = new TH2F ("hDispETRD"," dispersion^{2} vs E, EMCAL SM covered by TRD", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1539       fhDispETRD->SetYTitle("Dispersion^{2}");
1540       fhDispETRD->SetXTitle("E (GeV) ");
1541       outputContainer->Add(fhDispETRD);
1542
1543       if(!fRejectTrackMatch)
1544       {
1545         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);
1546         fhLam0ETMTRD->SetYTitle("#lambda_{0}^{2}");
1547         fhLam0ETMTRD->SetXTitle("E (GeV)");
1548         outputContainer->Add(fhLam0ETMTRD);
1549
1550         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);
1551         fhLam1ETMTRD->SetYTitle("#lambda_{1}^{2}");
1552         fhLam1ETMTRD->SetXTitle("E (GeV)");
1553         outputContainer->Add(fhLam1ETMTRD);
1554
1555         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);
1556         fhDispETMTRD->SetYTitle("Dispersion^{2}");
1557         fhDispETMTRD->SetXTitle("E (GeV) ");
1558         outputContainer->Add(fhDispETMTRD);
1559       }
1560     }
1561
1562     if(!fFillOnlySimpleSSHisto)
1563     {
1564       fhNCellsLam0LowE  = new TH2F ("hNCellsLam0LowE","N_{cells} in cluster vs #lambda_{0}^{2}, E < 2 GeV", nbins,nmin, nmax, ssbins,ssmin,ssmax);
1565       fhNCellsLam0LowE->SetXTitle("N_{cells}");
1566       fhNCellsLam0LowE->SetYTitle("#lambda_{0}^{2}");
1567       outputContainer->Add(fhNCellsLam0LowE);
1568
1569       fhNCellsLam0HighE  = new TH2F ("hNCellsLam0HighE","N_{cells} in cluster vs #lambda_{0}^{2}, E > 2 GeV", nbins,nmin, nmax, ssbins,ssmin,ssmax);
1570       fhNCellsLam0HighE->SetXTitle("N_{cells}");
1571       fhNCellsLam0HighE->SetYTitle("#lambda_{0}^{2}");
1572       outputContainer->Add(fhNCellsLam0HighE);
1573
1574       fhNCellsLam1LowE  = new TH2F ("hNCellsLam1LowE","N_{cells} in cluster vs #lambda_{1}^{2}, E < 2 GeV", nbins,nmin, nmax, ssbins,ssmin,ssmax);
1575       fhNCellsLam1LowE->SetXTitle("N_{cells}");
1576       fhNCellsLam1LowE->SetYTitle("#lambda_{0}^{2}");
1577       outputContainer->Add(fhNCellsLam1LowE);
1578
1579       fhNCellsLam1HighE  = new TH2F ("hNCellsLam1HighE","N_{cells} in cluster vs #lambda_{1}^{2}, E > 2 GeV", nbins,nmin, nmax, ssbins,ssmin,ssmax);
1580       fhNCellsLam1HighE->SetXTitle("N_{cells}");
1581       fhNCellsLam1HighE->SetYTitle("#lambda_{0}^{2}");
1582       outputContainer->Add(fhNCellsLam1HighE);
1583
1584       fhNCellsDispLowE  = new TH2F ("hNCellsDispLowE","N_{cells} in cluster vs dispersion^{2}, E < 2 GeV", nbins,nmin, nmax, ssbins,ssmin,ssmax);
1585       fhNCellsDispLowE->SetXTitle("N_{cells}");
1586       fhNCellsDispLowE->SetYTitle("D^{2}");
1587       outputContainer->Add(fhNCellsDispLowE);
1588
1589       fhNCellsDispHighE  = new TH2F ("hNCellsDispHighE","N_{cells} in cluster vs dispersion^{2}, E < 2 GeV", nbins,nmin, nmax, ssbins,ssmin,ssmax);
1590       fhNCellsDispHighE->SetXTitle("N_{cells}");
1591       fhNCellsDispHighE->SetYTitle("D^{2}");
1592       outputContainer->Add(fhNCellsDispHighE);
1593
1594       fhEtaLam0LowE  = new TH2F ("hEtaLam0LowE","#eta vs #lambda_{0}^{2}, E < 2 GeV", netabins,etamin,etamax, ssbins,ssmin,ssmax);
1595       fhEtaLam0LowE->SetYTitle("#lambda_{0}^{2}");
1596       fhEtaLam0LowE->SetXTitle("#eta");
1597       outputContainer->Add(fhEtaLam0LowE);
1598
1599       fhPhiLam0LowE  = new TH2F ("hPhiLam0LowE","#phi vs #lambda_{0}^{2}, E < 2 GeV", nphibins,phimin,phimax, ssbins,ssmin,ssmax);
1600       fhPhiLam0LowE->SetYTitle("#lambda_{0}^{2}");
1601       fhPhiLam0LowE->SetXTitle("#phi");
1602       outputContainer->Add(fhPhiLam0LowE);
1603
1604       fhEtaLam0HighE  = new TH2F ("hEtaLam0HighE","#eta vs #lambda_{0}^{2}, E > 2 GeV", netabins,etamin,etamax, ssbins,ssmin,ssmax);
1605       fhEtaLam0HighE->SetYTitle("#lambda_{0}^{2}");
1606       fhEtaLam0HighE->SetXTitle("#eta");
1607       outputContainer->Add(fhEtaLam0HighE);
1608
1609       fhPhiLam0HighE  = new TH2F ("hPhiLam0HighE","#phi vs #lambda_{0}^{2}, E > 2 GeV", nphibins,phimin,phimax, ssbins,ssmin,ssmax);
1610       fhPhiLam0HighE->SetYTitle("#lambda_{0}^{2}");
1611       fhPhiLam0HighE->SetXTitle("#phi");
1612       outputContainer->Add(fhPhiLam0HighE);
1613
1614       fhLam1Lam0LowE  = new TH2F ("hLam1Lam0LowE","#lambda_{0}^{2} vs #lambda_{1}^{2} in cluster of E < 2 GeV",  ssbins,ssmin,ssmax, ssbins,ssmin,ssmax);
1615       fhLam1Lam0LowE->SetYTitle("#lambda_{0}^{2}");
1616       fhLam1Lam0LowE->SetXTitle("#lambda_{1}^{2}");
1617       outputContainer->Add(fhLam1Lam0LowE);
1618
1619       fhLam1Lam0HighE  = new TH2F ("hLam1Lam0HighE","#lambda_{0}^{2} vs #lambda_{1}^{2} in cluster of E > 2 GeV",  ssbins,ssmin,ssmax, ssbins,ssmin,ssmax);
1620       fhLam1Lam0HighE->SetYTitle("#lambda_{0}^{2}");
1621       fhLam1Lam0HighE->SetXTitle("#lambda_{1}^{2}");
1622       outputContainer->Add(fhLam1Lam0HighE);
1623
1624       fhLam0DispLowE  = new TH2F ("hLam0DispLowE","#lambda_{0}^{2} vs dispersion^{2} in cluster of E < 2 GeV",  ssbins,ssmin,ssmax, ssbins,ssmin,ssmax);
1625       fhLam0DispLowE->SetXTitle("#lambda_{0}^{2}");
1626       fhLam0DispLowE->SetYTitle("D^{2}");
1627       outputContainer->Add(fhLam0DispLowE);
1628
1629       fhLam0DispHighE  = new TH2F ("hLam0DispHighE","#lambda_{0}^{2} vs dispersion^{2} in cluster of E > 2 GeV",  ssbins,ssmin,ssmax, ssbins,ssmin,ssmax);
1630       fhLam0DispHighE->SetXTitle("#lambda_{0}^{2}");
1631       fhLam0DispHighE->SetYTitle("D^{2}");
1632       outputContainer->Add(fhLam0DispHighE);
1633
1634       fhDispLam1LowE  = new TH2F ("hDispLam1LowE","Dispersion^{2} vs #lambda_{1}^{2} in cluster of E < 2 GeV",  ssbins,ssmin,ssmax, ssbins,ssmin,ssmax);
1635       fhDispLam1LowE->SetXTitle("D^{2}");
1636       fhDispLam1LowE->SetYTitle("#lambda_{1}^{2}");
1637       outputContainer->Add(fhDispLam1LowE);
1638
1639       fhDispLam1HighE  = new TH2F ("hDispLam1HighE","Dispersion^{2} vs #lambda_{1^{2}} in cluster of E > 2 GeV",  ssbins,ssmin,ssmax, ssbins,ssmin,ssmax);
1640       fhDispLam1HighE->SetXTitle("D^{2}");
1641       fhDispLam1HighE->SetYTitle("#lambda_{1}^{2}");
1642       outputContainer->Add(fhDispLam1HighE);
1643
1644       if(fCalorimeter == "EMCAL")
1645       {
1646         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);
1647         fhDispEtaE->SetXTitle("E (GeV)");
1648         fhDispEtaE->SetYTitle("#sigma^{2}_{#eta #eta}");
1649         outputContainer->Add(fhDispEtaE);
1650
1651         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);
1652         fhDispPhiE->SetXTitle("E (GeV)");
1653         fhDispPhiE->SetYTitle("#sigma^{2}_{#phi #phi}");
1654         outputContainer->Add(fhDispPhiE);
1655
1656         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);
1657         fhSumEtaE->SetXTitle("E (GeV)");
1658         fhSumEtaE->SetYTitle("#delta^{2}_{#eta #eta}");
1659         outputContainer->Add(fhSumEtaE);
1660
1661         fhSumPhiE  = new TH2F ("hSumPhiE","#delta^{2}_{#phi #phi} = #Sigma w_{i}(#phi_{i})^{2}/ #Sigma w_{i} - <#phi>^{2} vs E",
1662                                nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
1663         fhSumPhiE->SetXTitle("E (GeV)");
1664         fhSumPhiE->SetYTitle("#delta^{2}_{#phi #phi}");
1665         outputContainer->Add(fhSumPhiE);
1666
1667         fhSumEtaPhiE  = new TH2F ("hSumEtaPhiE","#delta^{2}_{#eta #phi} = #Sigma w_{i}(#phi_{i} #eta_{i} ) / #Sigma w_{i} - <#phi><#eta> vs E",
1668                                   nptbins,ptmin,ptmax, 2*ssbins,-ssmax,ssmax);
1669         fhSumEtaPhiE->SetXTitle("E (GeV)");
1670         fhSumEtaPhiE->SetYTitle("#delta^{2}_{#eta #phi}");
1671         outputContainer->Add(fhSumEtaPhiE);
1672
1673         fhDispEtaPhiDiffE  = new TH2F ("hDispEtaPhiDiffE","#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta} vs E",
1674                                        nptbins,ptmin,ptmax,200, -10,10);
1675         fhDispEtaPhiDiffE->SetXTitle("E (GeV)");
1676         fhDispEtaPhiDiffE->SetYTitle("#sigma^{2}_{#phi #phi}-#sigma^{2}_{#eta #eta}");
1677         outputContainer->Add(fhDispEtaPhiDiffE);
1678
1679         fhSphericityE  = new TH2F ("hSphericityE","(#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi}) vs E",
1680                                    nptbins,ptmin,ptmax, 200, -1,1);
1681         fhSphericityE->SetXTitle("E (GeV)");
1682         fhSphericityE->SetYTitle("s = (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi})");
1683         outputContainer->Add(fhSphericityE);
1684
1685         fhDispSumEtaDiffE  = new TH2F ("hDispSumEtaDiffE","#sigma^{2}_{#eta #eta} - #delta^{2}_{#eta #eta} / average vs E",  nptbins,ptmin,ptmax, 200,-0.01,0.01);
1686         fhDispSumEtaDiffE->SetXTitle("E (GeV)");
1687         fhDispSumEtaDiffE->SetYTitle("#sigma^{2}_{#eta #eta} - #delta^{2}_{#eta #eta} / average");
1688         outputContainer->Add(fhDispSumEtaDiffE);
1689
1690         fhDispSumPhiDiffE  = new TH2F ("hDispSumPhiDiffE","#sigma^{2}_{#phi #phi} - #delta^{2}_{#phi #phi} / average vs E",  nptbins,ptmin,ptmax, 200,-0.01,0.01);
1691         fhDispSumPhiDiffE->SetXTitle("E (GeV)");
1692         fhDispSumPhiDiffE->SetYTitle("#sigma^{2}_{#phi #phi} - #delta^{2}_{#phi #phi} / average");
1693         outputContainer->Add(fhDispSumPhiDiffE);
1694
1695         for(Int_t i = 0; i < 7; i++)
1696         {
1697           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]),
1698                                           ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
1699           fhDispEtaDispPhi[i]->SetXTitle("#sigma^{2}_{#eta #eta}");
1700           fhDispEtaDispPhi[i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1701           outputContainer->Add(fhDispEtaDispPhi[i]);
1702
1703           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]),
1704                                           ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
1705           fhLambda0DispEta[i]->SetXTitle("#lambda^{2}_{0}");
1706           fhLambda0DispEta[i]->SetYTitle("#sigma^{2}_{#eta #eta}");
1707           outputContainer->Add(fhLambda0DispEta[i]);
1708
1709           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]),
1710                                           ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
1711           fhLambda0DispPhi[i]->SetXTitle("#lambda^{2}_{0}");
1712           fhLambda0DispPhi[i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1713           outputContainer->Add(fhLambda0DispPhi[i]);
1714         }
1715       }
1716     }
1717   } // Shower shape
1718
1719   // Track Matching
1720
1721   if(fFillTMHisto)
1722   {
1723     fhTrackMatchedDEta[0]  = new TH2F
1724     ("hTrackMatchedDEtaNoCut",
1725      "d#eta of cluster-track vs cluster energy, no photon cuts",
1726      nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
1727     fhTrackMatchedDEta[0]->SetYTitle("d#eta");
1728     fhTrackMatchedDEta[0]->SetXTitle("E_{cluster} (GeV)");
1729
1730     fhTrackMatchedDPhi[0]  = new TH2F
1731     ("hTrackMatchedDPhiNoCut",
1732      "d#phi of cluster-track vs cluster energy, no photon cuts",
1733      nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
1734     fhTrackMatchedDPhi[0]->SetYTitle("d#phi (rad)");
1735     fhTrackMatchedDPhi[0]->SetXTitle("E_{cluster} (GeV)");
1736
1737     fhTrackMatchedDEtaDPhi[0]  = new TH2F
1738     ("hTrackMatchedDEtaDPhiNoCut",
1739      "d#eta vs d#phi of cluster-track vs cluster energy, no photon cuts",
1740      nresetabins,resetamin,resetamax,nresphibins,resphimin,resphimax);
1741     fhTrackMatchedDEtaDPhi[0]->SetYTitle("d#phi (rad)");
1742     fhTrackMatchedDEtaDPhi[0]->SetXTitle("d#eta");
1743
1744     fhdEdx[0]  = new TH2F ("hdEdxNoCut","matched track <dE/dx> vs cluster E, no photon cuts ",
1745                            nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
1746     fhdEdx[0]->SetXTitle("E (GeV)");
1747     fhdEdx[0]->SetYTitle("<dE/dx>");
1748
1749     fhEOverP[0]  = new TH2F ("hEOverPNoCut","matched track E/p vs cluster E, no photon cuts ",
1750                              nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
1751     fhEOverP[0]->SetXTitle("E (GeV)");
1752     fhEOverP[0]->SetYTitle("E/p");
1753
1754     outputContainer->Add(fhTrackMatchedDEta[0]) ;
1755     outputContainer->Add(fhTrackMatchedDPhi[0]) ;
1756     outputContainer->Add(fhTrackMatchedDEtaDPhi[0]) ;
1757     outputContainer->Add(fhdEdx[0]);
1758     outputContainer->Add(fhEOverP[0]);
1759
1760     fhTrackMatchedDEta[1]  = new TH2F
1761     ("hTrackMatchedDEta",
1762      "d#eta of cluster-track vs cluster energy, no photon cuts",
1763      nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
1764     fhTrackMatchedDEta[1]->SetYTitle("d#eta");
1765     fhTrackMatchedDEta[1]->SetXTitle("E_{cluster} (GeV)");
1766
1767     fhTrackMatchedDPhi[1]  = new TH2F
1768     ("hTrackMatchedDPhi",
1769      "d#phi of cluster-track vs cluster energy, no photon cuts",
1770      nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
1771     fhTrackMatchedDPhi[1]->SetYTitle("d#phi (rad)");
1772     fhTrackMatchedDPhi[1]->SetXTitle("E_{cluster} (GeV)");
1773
1774     fhTrackMatchedDEtaDPhi[1]  = new TH2F
1775     ("hTrackMatchedDEtaDPhi",
1776      "d#eta vs d#phi of cluster-track vs cluster energy, no photon cuts",
1777      nresetabins,resetamin,resetamax,nresphibins,resphimin,resphimax);
1778     fhTrackMatchedDEtaDPhi[1]->SetYTitle("d#phi (rad)");
1779     fhTrackMatchedDEtaDPhi[1]->SetXTitle("d#eta");
1780
1781     fhdEdx[1]  = new TH2F ("hdEdx","matched track <dE/dx> vs cluster E ",
1782                            nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
1783     fhdEdx[1]->SetXTitle("E (GeV)");
1784     fhdEdx[1]->SetYTitle("<dE/dx>");
1785
1786     fhEOverP[1]  = new TH2F ("hEOverP","matched track E/p vs cluster E ",
1787                              nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
1788     fhEOverP[1]->SetXTitle("E (GeV)");
1789     fhEOverP[1]->SetYTitle("E/p");
1790
1791     outputContainer->Add(fhTrackMatchedDEta[1]) ;
1792     outputContainer->Add(fhTrackMatchedDPhi[1]) ;
1793     outputContainer->Add(fhTrackMatchedDEtaDPhi[1]) ;
1794     outputContainer->Add(fhdEdx[1]);
1795     outputContainer->Add(fhEOverP[1]);
1796
1797     if(fCalorimeter=="EMCAL")
1798     {
1799       fhTrackMatchedDEtaTRD[0]  = new TH2F
1800       ("hTrackMatchedDEtaTRDNoCut",
1801        "d#eta of cluster-track vs cluster energy, SM behind TRD, no photon cuts",
1802        nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
1803       fhTrackMatchedDEtaTRD[0]->SetYTitle("d#eta");
1804       fhTrackMatchedDEtaTRD[0]->SetXTitle("E_{cluster} (GeV)");
1805
1806       fhTrackMatchedDPhiTRD[0]  = new TH2F
1807       ("hTrackMatchedDPhiTRDNoCut",
1808        "d#phi of cluster-track vs cluster energy, SM behing TRD, no photon cuts",
1809        nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
1810       fhTrackMatchedDPhiTRD[0]->SetYTitle("d#phi (rad)");
1811       fhTrackMatchedDPhiTRD[0]->SetXTitle("E_{cluster} (GeV)");
1812
1813       fhEOverPTRD[0]  = new TH2F ("hEOverPTRDNoCut","matched track E/p vs cluster E, behind TRD, no photon cuts ",
1814                                   nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
1815       fhEOverPTRD[0]->SetXTitle("E (GeV)");
1816       fhEOverPTRD[0]->SetYTitle("E/p");
1817
1818       outputContainer->Add(fhTrackMatchedDEtaTRD[0]) ;
1819       outputContainer->Add(fhTrackMatchedDPhiTRD[0]) ;
1820       outputContainer->Add(fhEOverPTRD[0]);
1821
1822       fhTrackMatchedDEtaTRD[1]  = new TH2F
1823       ("hTrackMatchedDEtaTRD",
1824        "d#eta of cluster-track vs cluster energy, SM behind TRD",
1825        nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
1826       fhTrackMatchedDEtaTRD[1]->SetYTitle("d#eta");
1827       fhTrackMatchedDEtaTRD[1]->SetXTitle("E_{cluster} (GeV)");
1828
1829       fhTrackMatchedDPhiTRD[1]  = new TH2F
1830       ("hTrackMatchedDPhiTRD",
1831        "d#phi of cluster-track vs cluster energy, SM behing TRD",
1832        nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
1833       fhTrackMatchedDPhiTRD[1]->SetYTitle("d#phi (rad)");
1834       fhTrackMatchedDPhiTRD[1]->SetXTitle("E_{cluster} (GeV)");
1835
1836       fhEOverPTRD[1]  = new TH2F ("hEOverPTRD","matched track E/p vs cluster E, behind TRD ",
1837                                   nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
1838       fhEOverPTRD[1]->SetXTitle("E (GeV)");
1839       fhEOverPTRD[1]->SetYTitle("E/p");
1840
1841       outputContainer->Add(fhTrackMatchedDEtaTRD[1]) ;
1842       outputContainer->Add(fhTrackMatchedDPhiTRD[1]) ;
1843       outputContainer->Add(fhEOverPTRD[1]);
1844
1845     }
1846
1847     if(IsDataMC())
1848     {
1849       fhTrackMatchedDEtaMCNoOverlap[0]  = new TH2F
1850       ("hTrackMatchedDEtaMCNoOverlapNoCut",
1851        "d#eta of cluster-track vs cluster energy, no other MC particles overlap",
1852        nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
1853       fhTrackMatchedDEtaMCNoOverlap[0]->SetYTitle("d#eta");
1854       fhTrackMatchedDEtaMCNoOverlap[0]->SetXTitle("E_{cluster} (GeV)");
1855
1856       fhTrackMatchedDPhiMCNoOverlap[0]  = new TH2F
1857       ("hTrackMatchedDPhiMCNoOverlapNoCut",
1858        "d#phi of cluster-track vs cluster energy, no other MC particles overlap",
1859        nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
1860       fhTrackMatchedDPhiMCNoOverlap[0]->SetYTitle("d#phi (rad)");
1861       fhTrackMatchedDPhiMCNoOverlap[0]->SetXTitle("E_{cluster} (GeV)");
1862
1863       outputContainer->Add(fhTrackMatchedDEtaMCNoOverlap[0]) ;
1864       outputContainer->Add(fhTrackMatchedDPhiMCNoOverlap[0]) ;
1865
1866       fhTrackMatchedDEtaMCNoOverlap[1]  = new TH2F
1867       ("hTrackMatchedDEtaMCNoOverlap",
1868        "d#eta of cluster-track vs cluster energy, no other MC particles overlap",
1869        nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
1870       fhTrackMatchedDEtaMCNoOverlap[1]->SetYTitle("d#eta");
1871       fhTrackMatchedDEtaMCNoOverlap[1]->SetXTitle("E_{cluster} (GeV)");
1872
1873       fhTrackMatchedDPhiMCNoOverlap[1]  = new TH2F
1874       ("hTrackMatchedDPhiMCNoOverlap",
1875        "d#phi of cluster-track vs cluster energy, no other MC particles overlap",
1876        nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
1877       fhTrackMatchedDPhiMCNoOverlap[1]->SetYTitle("d#phi (rad)");
1878       fhTrackMatchedDPhiMCNoOverlap[1]->SetXTitle("E_{cluster} (GeV)");
1879
1880       outputContainer->Add(fhTrackMatchedDEtaMCNoOverlap[1]) ;
1881       outputContainer->Add(fhTrackMatchedDPhiMCNoOverlap[1]) ;
1882
1883       fhTrackMatchedDEtaMCOverlap[0]  = new TH2F
1884       ("hTrackMatchedDEtaMCOverlapNoCut",
1885        "d#eta of cluster-track vs cluster energy, several MC particles overlap",
1886        nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
1887       fhTrackMatchedDEtaMCOverlap[0]->SetYTitle("d#eta");
1888       fhTrackMatchedDEtaMCOverlap[0]->SetXTitle("E_{cluster} (GeV)");
1889
1890       fhTrackMatchedDPhiMCOverlap[0]  = new TH2F
1891       ("hTrackMatchedDPhiMCOverlapNoCut",
1892        "d#phi of cluster-track vs cluster energy, several MC particles overlap",
1893        nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
1894       fhTrackMatchedDPhiMCOverlap[0]->SetYTitle("d#phi (rad)");
1895       fhTrackMatchedDPhiMCOverlap[0]->SetXTitle("E_{cluster} (GeV)");
1896
1897       outputContainer->Add(fhTrackMatchedDEtaMCOverlap[0]) ;
1898       outputContainer->Add(fhTrackMatchedDPhiMCOverlap[0]) ;
1899
1900       fhTrackMatchedDEtaMCOverlap[1]  = new TH2F
1901       ("hTrackMatchedDEtaMCOverlap",
1902        "d#eta of cluster-track vs cluster energy, several MC particles overlap",
1903        nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
1904       fhTrackMatchedDEtaMCOverlap[1]->SetYTitle("d#eta");
1905       fhTrackMatchedDEtaMCOverlap[1]->SetXTitle("E_{cluster} (GeV)");
1906
1907       fhTrackMatchedDPhiMCOverlap[1]  = new TH2F
1908       ("hTrackMatchedDPhiMCOverlap",
1909        "d#phi of cluster-track vs cluster energy, several MC particles overlap",
1910        nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
1911       fhTrackMatchedDPhiMCOverlap[1]->SetYTitle("d#phi (rad)");
1912       fhTrackMatchedDPhiMCOverlap[1]->SetXTitle("E_{cluster} (GeV)");
1913
1914       outputContainer->Add(fhTrackMatchedDEtaMCOverlap[1]) ;
1915       outputContainer->Add(fhTrackMatchedDPhiMCOverlap[1]) ;
1916
1917       fhTrackMatchedDEtaMCConversion[0]  = new TH2F
1918       ("hTrackMatchedDEtaMCConversionNoCut",
1919        "d#eta of cluster-track vs cluster energy, no other MC particles overlap appart from conversions",
1920        nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
1921       fhTrackMatchedDEtaMCConversion[0]->SetYTitle("d#eta");
1922       fhTrackMatchedDEtaMCConversion[0]->SetXTitle("E_{cluster} (GeV)");
1923
1924       fhTrackMatchedDPhiMCConversion[0]  = new TH2F
1925       ("hTrackMatchedDPhiMCConversionNoCut",
1926        "d#phi of cluster-track vs cluster energy, no other MC particles overlap appart from conversions",
1927        nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
1928       fhTrackMatchedDPhiMCConversion[0]->SetYTitle("d#phi (rad)");
1929       fhTrackMatchedDPhiMCConversion[0]->SetXTitle("E_{cluster} (GeV)");
1930
1931       outputContainer->Add(fhTrackMatchedDEtaMCConversion[0]) ;
1932       outputContainer->Add(fhTrackMatchedDPhiMCConversion[0]) ;
1933
1934
1935       fhTrackMatchedDEtaMCConversion[1]  = new TH2F
1936       ("hTrackMatchedDEtaMCConversion",
1937        "d#eta of cluster-track vs cluster energy, no other MC particles overlap appart from conversions",
1938        nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
1939       fhTrackMatchedDEtaMCConversion[1]->SetYTitle("d#eta");
1940       fhTrackMatchedDEtaMCConversion[1]->SetXTitle("E_{cluster} (GeV)");
1941
1942       fhTrackMatchedDPhiMCConversion[1]  = new TH2F
1943       ("hTrackMatchedDPhiMCConversion",
1944        "d#phi of cluster-track vs cluster energy, no other MC particles overlap appart from conversions",
1945        nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
1946       fhTrackMatchedDPhiMCConversion[1]->SetYTitle("d#phi (rad)");
1947       fhTrackMatchedDPhiMCConversion[1]->SetXTitle("E_{cluster} (GeV)");
1948
1949       outputContainer->Add(fhTrackMatchedDEtaMCConversion[1]) ;
1950       outputContainer->Add(fhTrackMatchedDPhiMCConversion[1]) ;
1951
1952
1953       fhTrackMatchedMCParticle[0]  = new TH2F
1954       ("hTrackMatchedMCParticleNoCut",
1955        "Origin of particle vs energy",
1956        nptbins,ptmin,ptmax,8,0,8);
1957       fhTrackMatchedMCParticle[0]->SetXTitle("E (GeV)");
1958       //fhTrackMatchedMCParticle[0]->SetYTitle("Particle type");
1959
1960       fhTrackMatchedMCParticle[0]->GetYaxis()->SetBinLabel(1 ,"Photon");
1961       fhTrackMatchedMCParticle[0]->GetYaxis()->SetBinLabel(2 ,"Electron");
1962       fhTrackMatchedMCParticle[0]->GetYaxis()->SetBinLabel(3 ,"Meson Merged");
1963       fhTrackMatchedMCParticle[0]->GetYaxis()->SetBinLabel(4 ,"Rest");
1964       fhTrackMatchedMCParticle[0]->GetYaxis()->SetBinLabel(5 ,"Conv. Photon");
1965       fhTrackMatchedMCParticle[0]->GetYaxis()->SetBinLabel(6 ,"Conv. Electron");
1966       fhTrackMatchedMCParticle[0]->GetYaxis()->SetBinLabel(7 ,"Conv. Merged");
1967       fhTrackMatchedMCParticle[0]->GetYaxis()->SetBinLabel(8 ,"Conv. Rest");
1968
1969       fhTrackMatchedMCParticle[1]  = new TH2F
1970       ("hTrackMatchedMCParticle",
1971        "Origin of particle vs energy",
1972        nptbins,ptmin,ptmax,8,0,8);
1973       fhTrackMatchedMCParticle[1]->SetXTitle("E (GeV)");
1974       //fhTrackMatchedMCParticle[1]->SetYTitle("Particle type");
1975
1976       fhTrackMatchedMCParticle[1]->GetYaxis()->SetBinLabel(1 ,"Photon");
1977       fhTrackMatchedMCParticle[1]->GetYaxis()->SetBinLabel(2 ,"Electron");
1978       fhTrackMatchedMCParticle[1]->GetYaxis()->SetBinLabel(3 ,"Meson Merged");
1979       fhTrackMatchedMCParticle[1]->GetYaxis()->SetBinLabel(4 ,"Rest");
1980       fhTrackMatchedMCParticle[1]->GetYaxis()->SetBinLabel(5 ,"Conv. Photon");
1981       fhTrackMatchedMCParticle[1]->GetYaxis()->SetBinLabel(6 ,"Conv. Electron");
1982       fhTrackMatchedMCParticle[1]->GetYaxis()->SetBinLabel(7 ,"Conv. Merged");
1983       fhTrackMatchedMCParticle[1]->GetYaxis()->SetBinLabel(8 ,"Conv. Rest");
1984
1985       outputContainer->Add(fhTrackMatchedMCParticle[0]);
1986       outputContainer->Add(fhTrackMatchedMCParticle[1]);
1987
1988     }
1989   }
1990
1991   if(fFillPileUpHistograms)
1992   {
1993
1994     TString pileUpName[] = {"SPD","EMCAL","SPDOrEMCAL","SPDAndEMCAL","SPDAndNotEMCAL","EMCALAndNotSPD","NotSPDAndNotEMCAL"} ;
1995
1996     for(Int_t i = 0 ; i < 7 ; i++)
1997     {
1998       fhPtPileUp[i]  = new TH1F(Form("hPtPileUp%s",pileUpName[i].Data()),
1999                                        Form("Cluster  p_{T} distribution, %s Pile-Up event",pileUpName[i].Data()), nptbins,ptmin,ptmax);
2000       fhPtPileUp[i]->SetXTitle("p_{T} (GeV/c)");
2001       outputContainer->Add(fhPtPileUp[i]);
2002
2003       fhPtChargedPileUp[i]  = new TH1F(Form("hPtChargedPileUp%s",pileUpName[i].Data()),
2004                                       Form("Charged clusters p_{T} distribution, %s Pile-Up event",pileUpName[i].Data()), nptbins,ptmin,ptmax);
2005       fhPtChargedPileUp[i]->SetXTitle("p_{T} (GeV/c)");
2006       outputContainer->Add(fhPtChargedPileUp[i]);
2007
2008       fhPtPhotonPileUp[i]  = new TH1F(Form("hPtPhotonPileUp%s",pileUpName[i].Data()),
2009                                       Form("Selected photon p_{T} distribution, %s Pile-Up event",pileUpName[i].Data()), nptbins,ptmin,ptmax);
2010       fhPtPhotonPileUp[i]->SetXTitle("p_{T} (GeV/c)");
2011       outputContainer->Add(fhPtPhotonPileUp[i]);
2012
2013
2014       fhClusterEFracLongTimePileUp[i]  = new TH2F(Form("hClusterEFracLongTimePileUp%s",pileUpName[i].Data()),
2015                                              Form("Cluster E vs fraction of cluster energy from large T cells, %s Pile-Up event",pileUpName[i].Data()),
2016                                              nptbins,ptmin,ptmax,200,0,1);
2017       fhClusterEFracLongTimePileUp[i]->SetXTitle("E (GeV)");
2018       fhClusterEFracLongTimePileUp[i]->SetYTitle("E(large time) / E");
2019       outputContainer->Add(fhClusterEFracLongTimePileUp[i]);
2020
2021       fhClusterTimeDiffPileUp[i]  = new TH2F(Form("hClusterTimeDiffPileUp%s",pileUpName[i].Data()),
2022                                 Form("Cluster E vs t_{max}-t_{cell} in cluster, %s Pile-Up event",pileUpName[i].Data()),
2023                                              nptbins,ptmin,ptmax,200,-100,100);
2024       fhClusterTimeDiffPileUp[i]->SetXTitle("E (GeV)");
2025       fhClusterTimeDiffPileUp[i]->SetYTitle("t_{max}-t_{cell} (ns)");
2026       outputContainer->Add(fhClusterTimeDiffPileUp[i]);
2027
2028       fhClusterTimeDiffChargedPileUp[i]  = new TH2F(Form("hClusterTimeDiffChargedPileUp%s",pileUpName[i].Data()),
2029                                        Form("Charged clusters E vs t_{max}-t_{cell} in cluster, %s Pile-Up event",pileUpName[i].Data()),
2030                                                     nptbins,ptmin,ptmax,200,-100,100);
2031       fhClusterTimeDiffChargedPileUp[i]->SetXTitle("E (GeV)");
2032       fhClusterTimeDiffChargedPileUp[i]->SetYTitle("t_{max}-t_{cell} (ns)");
2033       outputContainer->Add(fhClusterTimeDiffChargedPileUp[i]);
2034
2035       fhClusterTimeDiffPhotonPileUp[i]  = new TH2F(Form("hClusterTimeDiffPhotonPileUp%s",pileUpName[i].Data()),
2036                                       Form("Selected photon E vs t_{max}-t_{cell} in cluster, %s Pile-Up event",pileUpName[i].Data()),
2037                                                    nptbins,ptmin,ptmax,200,-100,100);
2038       fhClusterTimeDiffPhotonPileUp[i]->SetXTitle("E (GeV)");
2039       fhClusterTimeDiffPhotonPileUp[i]->SetYTitle("t_{max}-t_{cell} (ns)");
2040       outputContainer->Add(fhClusterTimeDiffPhotonPileUp[i]);
2041
2042       fhLambda0PileUp[i]  = new TH2F(Form("hLambda0PileUp%s",pileUpName[i].Data()),
2043                                      Form("Cluster E vs #lambda^{2}_{0} in cluster, %s Pile-Up event",pileUpName[i].Data()),
2044                                      nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2045       fhLambda0PileUp[i]->SetXTitle("E (GeV)");
2046       fhLambda0PileUp[i]->SetYTitle("#lambda^{2}_{0}");
2047       outputContainer->Add(fhLambda0PileUp[i]);
2048
2049       fhLambda0ChargedPileUp[i]  = new TH2F(Form("hLambda0ChargedPileUp%s",pileUpName[i].Data()),
2050                                                     Form("Charged clusters E vs #lambda^{2}_{0}in cluster, %s Pile-Up event",pileUpName[i].Data()), nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2051       fhLambda0ChargedPileUp[i]->SetXTitle("E (GeV)");
2052       fhLambda0ChargedPileUp[i]->SetYTitle("#lambda^{2}_{0}");
2053       outputContainer->Add(fhLambda0ChargedPileUp[i]);
2054
2055     }
2056
2057     fhEtaPhiBC0  = new TH2F ("hEtaPhiBC0","eta-phi for clusters tof corresponding to BC=0",netabins,etamin,etamax, nphibins,phimin,phimax);
2058     fhEtaPhiBC0->SetXTitle("#eta ");
2059     fhEtaPhiBC0->SetYTitle("#phi (rad)");
2060     outputContainer->Add(fhEtaPhiBC0);
2061
2062     fhEtaPhiBCPlus  = new TH2F ("hEtaPhiBCPlus","eta-phi for clusters tof corresponding to BC>0",netabins,etamin,etamax, nphibins,phimin,phimax);
2063     fhEtaPhiBCPlus->SetXTitle("#eta ");
2064     fhEtaPhiBCPlus->SetYTitle("#phi (rad)");
2065     outputContainer->Add(fhEtaPhiBCPlus);
2066
2067     fhEtaPhiBCMinus  = new TH2F ("hEtaPhiBCMinus","eta-phi for clusters tof corresponding to BC<0",netabins,etamin,etamax, nphibins,phimin,phimax);
2068     fhEtaPhiBCMinus->SetXTitle("#eta ");
2069     fhEtaPhiBCMinus->SetYTitle("#phi (rad)");
2070     outputContainer->Add(fhEtaPhiBCMinus);
2071
2072     fhEtaPhiBC0PileUpSPD  = new TH2F ("hEtaPhiBC0PileUpSPD","eta-phi for clusters tof corresponding to BC=0, SPD pile-up",netabins,etamin,etamax, nphibins,phimin,phimax);
2073     fhEtaPhiBC0PileUpSPD->SetXTitle("#eta ");
2074     fhEtaPhiBC0PileUpSPD->SetYTitle("#phi (rad)");
2075     outputContainer->Add(fhEtaPhiBC0PileUpSPD);
2076
2077     fhEtaPhiBCPlusPileUpSPD  = new TH2F ("hEtaPhiBCPlusPileUpSPD","eta-phi for clusters tof corresponding to BC>0, SPD pile-up",netabins,etamin,etamax, nphibins,phimin,phimax);
2078     fhEtaPhiBCPlusPileUpSPD->SetXTitle("#eta ");
2079     fhEtaPhiBCPlusPileUpSPD->SetYTitle("#phi (rad)");
2080     outputContainer->Add(fhEtaPhiBCPlusPileUpSPD);
2081
2082     fhEtaPhiBCMinusPileUpSPD  = new TH2F ("hEtaPhiBCMinusPileUpSPD","eta-phi for clusters tof corresponding to BC<0, SPD pile-up",netabins,etamin,etamax, nphibins,phimin,phimax);
2083     fhEtaPhiBCMinusPileUpSPD->SetXTitle("#eta ");
2084     fhEtaPhiBCMinusPileUpSPD->SetYTitle("#phi (rad)");
2085     outputContainer->Add(fhEtaPhiBCMinusPileUpSPD);
2086
2087     fhTimeENoCut  = new TH2F ("hTimeE_NoCut","time of cluster vs E of clusters, no cut", nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
2088     fhTimeENoCut->SetXTitle("E (GeV)");
2089     fhTimeENoCut->SetYTitle("time (ns)");
2090     outputContainer->Add(fhTimeENoCut);
2091
2092     fhTimeESPD  = new TH2F ("hTimeE_SPD","time of cluster vs E of clusters, SPD cut", nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
2093     fhTimeESPD->SetXTitle("E (GeV)");
2094     fhTimeESPD->SetYTitle("time (ns)");
2095     outputContainer->Add(fhTimeESPD);
2096
2097     fhTimeESPDMulti  = new TH2F ("hTimeE_SPDMulti","time of cluster vs E of clusters, SPD multi cut", nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
2098     fhTimeESPDMulti->SetXTitle("E (GeV)");
2099     fhTimeESPDMulti->SetYTitle("time (ns)");
2100     outputContainer->Add(fhTimeESPDMulti);
2101
2102     fhTimeNPileUpVertSPD  = new TH2F ("hTime_NPileUpVertSPD","time of cluster vs N pile-up SPD vertex", ntimebins,timemin,timemax,50,0,50);
2103     fhTimeNPileUpVertSPD->SetYTitle("# vertex ");
2104     fhTimeNPileUpVertSPD->SetXTitle("time (ns)");
2105     outputContainer->Add(fhTimeNPileUpVertSPD);
2106
2107     fhTimeNPileUpVertTrack  = new TH2F ("hTime_NPileUpVertTracks","time of cluster vs N pile-up Tracks vertex", ntimebins,timemin,timemax, 50,0,50 );
2108     fhTimeNPileUpVertTrack->SetYTitle("# vertex ");
2109     fhTimeNPileUpVertTrack->SetXTitle("time (ns)");
2110     outputContainer->Add(fhTimeNPileUpVertTrack);
2111
2112     fhTimeNPileUpVertContributors  = new TH2F ("hTime_NPileUpVertContributors","time of cluster vs N constributors to pile-up SPD vertex", ntimebins,timemin,timemax,50,0,50);
2113     fhTimeNPileUpVertContributors->SetYTitle("# vertex ");
2114     fhTimeNPileUpVertContributors->SetXTitle("time (ns)");
2115     outputContainer->Add(fhTimeNPileUpVertContributors);
2116
2117     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);
2118     fhTimePileUpMainVertexZDistance->SetYTitle("distance Z (cm) ");
2119     fhTimePileUpMainVertexZDistance->SetXTitle("time (ns)");
2120     outputContainer->Add(fhTimePileUpMainVertexZDistance);
2121
2122     fhTimePileUpMainVertexZDiamond  = new TH2F ("hTime_PileUpMainVertexZDiamond","time of cluster vs distance in Z pile-up SPD vertex - z diamond",ntimebins,timemin,timemax,100,0,50);
2123     fhTimePileUpMainVertexZDiamond->SetYTitle("diamond distance Z (cm) ");
2124     fhTimePileUpMainVertexZDiamond->SetXTitle("time (ns)");
2125     outputContainer->Add(fhTimePileUpMainVertexZDiamond);
2126
2127     TString title[] = {"no |t diff| cut","|t diff|<20 ns","|t diff|>20 ns","|t diff|>40 ns"};
2128     TString name [] = {"TDiffNoCut","TDiffSmaller20ns","TDiffLarger20ns","TDiffLarger40ns"};
2129     for(Int_t i = 0; i < 4; i++)
2130     {
2131       fhClusterMultSPDPileUp[i] = new TH2F(Form("fhClusterMultSPDPileUp_%s", name[i].Data()),
2132                                            Form("Number of clusters per pile up event with E > 0.5 and %s respect cluster max vs cluster max E ",title[i].Data()),
2133                                            nptbins,ptmin,ptmax,100,0,100);
2134       fhClusterMultSPDPileUp[i]->SetYTitle("n clusters ");
2135       fhClusterMultSPDPileUp[i]->SetXTitle("E_{cluster max} (GeV)");
2136       outputContainer->Add(fhClusterMultSPDPileUp[i]) ;
2137
2138       fhClusterMultNoPileUp[i] = new TH2F(Form("fhClusterMultNoPileUp_%s", name[i].Data()),
2139                                           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()),
2140                                           nptbins,ptmin,ptmax,100,0,100);
2141       fhClusterMultNoPileUp[i]->SetYTitle("n clusters ");
2142       fhClusterMultNoPileUp[i]->SetXTitle("E_{cluster max} (GeV)");
2143       outputContainer->Add(fhClusterMultNoPileUp[i]) ;
2144     }
2145
2146   }
2147
2148   if(IsDataMC())
2149   {
2150     TString ptype[] = { "#gamma", "#gamma_{#pi decay}","#gamma_{other decay}", "#pi^{0}","#eta",
2151                         "e^{#pm}","#gamma->e^{#pm}","hadron?","Anti-N","Anti-P",
2152                         "#gamma_{prompt}","#gamma_{fragmentation}","#gamma_{ISR}","String"                                } ;
2153
2154     TString pname[] = { "Photon","PhotonPi0Decay","PhotonOtherDecay","Pi0","Eta","Electron",
2155                         "Conversion", "Hadron", "AntiNeutron","AntiProton",
2156                         "PhotonPrompt","PhotonFragmentation","PhotonISR","String" } ;
2157
2158     for(Int_t i = 0; i < fNOriginHistograms; i++)
2159     {
2160       fhMCE[i]  = new TH1F(Form("hE_MC%s",pname[i].Data()),
2161                                 Form("cluster from %s : E ",ptype[i].Data()),
2162                                 nptbins,ptmin,ptmax);
2163       fhMCE[i]->SetXTitle("E (GeV)");
2164       outputContainer->Add(fhMCE[i]) ;
2165
2166       fhMCPt[i]  = new TH1F(Form("hPt_MC%s",pname[i].Data()),
2167                            Form("cluster from %s : p_{T} ",ptype[i].Data()),
2168                            nptbins,ptmin,ptmax);
2169       fhMCPt[i]->SetXTitle("p_{T} (GeV/c)");
2170       outputContainer->Add(fhMCPt[i]) ;
2171
2172       fhMCEta[i]  = new TH2F(Form("hEta_MC%s",pname[i].Data()),
2173                            Form("cluster from %s : #eta ",ptype[i].Data()),
2174                            nptbins,ptmin,ptmax,netabins,etamin,etamax);
2175       fhMCEta[i]->SetYTitle("#eta");
2176       fhMCEta[i]->SetXTitle("E (GeV)");
2177       outputContainer->Add(fhMCEta[i]) ;
2178
2179       fhMCPhi[i]  = new TH2F(Form("hPhi_MC%s",pname[i].Data()),
2180                            Form("cluster from %s : #phi ",ptype[i].Data()),
2181                            nptbins,ptmin,ptmax,nphibins,phimin,phimax);
2182       fhMCPhi[i]->SetYTitle("#phi (rad)");
2183       fhMCPhi[i]->SetXTitle("E (GeV)");
2184       outputContainer->Add(fhMCPhi[i]) ;
2185
2186
2187       fhMCDeltaE[i]  = new TH2F (Form("hDeltaE_MC%s",pname[i].Data()),
2188                                  Form("MC - Reco E from %s",pname[i].Data()),
2189                                  nptbins,ptmin,ptmax, 200,-50,50);
2190       fhMCDeltaE[i]->SetXTitle("#Delta E (GeV)");
2191       outputContainer->Add(fhMCDeltaE[i]);
2192
2193       fhMCDeltaPt[i]  = new TH2F (Form("hDeltaPt_MC%s",pname[i].Data()),
2194                                   Form("MC - Reco p_{T} from %s",pname[i].Data()),
2195                                   nptbins,ptmin,ptmax, 200,-50,50);
2196       fhMCDeltaPt[i]->SetXTitle("#Delta p_{T} (GeV/c)");
2197       outputContainer->Add(fhMCDeltaPt[i]);
2198
2199       fhMC2E[i]  = new TH2F (Form("h2E_MC%s",pname[i].Data()),
2200                              Form("E distribution, reconstructed vs generated from %s",pname[i].Data()),
2201                              nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2202       fhMC2E[i]->SetXTitle("E_{rec} (GeV)");
2203       fhMC2E[i]->SetYTitle("E_{gen} (GeV)");
2204       outputContainer->Add(fhMC2E[i]);
2205
2206       fhMC2Pt[i]  = new TH2F (Form("h2Pt_MC%s",pname[i].Data()),
2207                               Form("p_T distribution, reconstructed vs generated from %s",pname[i].Data()),
2208                               nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2209       fhMC2Pt[i]->SetXTitle("p_{T,rec} (GeV/c)");
2210       fhMC2Pt[i]->SetYTitle("p_{T,gen} (GeV/c)");
2211       outputContainer->Add(fhMC2Pt[i]);
2212
2213
2214     }
2215
2216     TString pptype[] = { "#gamma", "#gamma_{#pi decay}","#gamma_{other decay}","hadron?",
2217                          "#gamma_{prompt}","#gamma_{fragmentation}","#gamma_{ISR}"} ;
2218
2219     TString ppname[] = { "Photon","PhotonPi0Decay","PhotonOtherDecay","Hadron",
2220                          "PhotonPrompt","PhotonFragmentation","PhotonISR"} ;
2221
2222     for(Int_t i = 0; i < fNPrimaryHistograms; i++)
2223     {
2224       fhEPrimMC[i]  = new TH1F(Form("hEPrim_MC%s",ppname[i].Data()),
2225                            Form("primary photon %s : E ",pptype[i].Data()),
2226                            nptbins,ptmin,ptmax);
2227       fhEPrimMC[i]->SetXTitle("E (GeV)");
2228       outputContainer->Add(fhEPrimMC[i]) ;
2229
2230       fhPtPrimMC[i]  = new TH1F(Form("hPtPrim_MC%s",ppname[i].Data()),
2231                             Form("primary photon %s : p_{T} ",pptype[i].Data()),
2232                             nptbins,ptmin,ptmax);
2233       fhPtPrimMC[i]->SetXTitle("p_{T} (GeV/c)");
2234       outputContainer->Add(fhPtPrimMC[i]) ;
2235
2236       fhYPrimMC[i]  = new TH2F(Form("hYPrim_MC%s",ppname[i].Data()),
2237                              Form("primary photon %s : Rapidity ",pptype[i].Data()),
2238                              nptbins,ptmin,ptmax,800,-8,8);
2239       fhYPrimMC[i]->SetYTitle("Rapidity");
2240       fhYPrimMC[i]->SetXTitle("E (GeV)");
2241       outputContainer->Add(fhYPrimMC[i]) ;
2242
2243       fhPhiPrimMC[i]  = new TH2F(Form("hPhiPrim_MC%s",ppname[i].Data()),
2244                              Form("primary photon %s : #phi ",pptype[i].Data()),
2245                              nptbins,ptmin,ptmax,nphibins,phimin,phimax);
2246       fhPhiPrimMC[i]->SetYTitle("#phi (rad)");
2247       fhPhiPrimMC[i]->SetXTitle("E (GeV)");
2248       outputContainer->Add(fhPhiPrimMC[i]) ;
2249
2250
2251       fhEPrimMCAcc[i]  = new TH1F(Form("hEPrimAcc_MC%s",ppname[i].Data()),
2252                                Form("primary photon %s in acceptance: E ",pptype[i].Data()),
2253                                nptbins,ptmin,ptmax);
2254       fhEPrimMCAcc[i]->SetXTitle("E (GeV)");
2255       outputContainer->Add(fhEPrimMCAcc[i]) ;
2256
2257       fhPtPrimMCAcc[i]  = new TH1F(Form("hPtPrimAcc_MC%s",ppname[i].Data()),
2258                                 Form("primary photon %s in acceptance: p_{T} ",pptype[i].Data()),
2259                                 nptbins,ptmin,ptmax);
2260       fhPtPrimMCAcc[i]->SetXTitle("p_{T} (GeV/c)");
2261       outputContainer->Add(fhPtPrimMCAcc[i]) ;
2262
2263       fhYPrimMCAcc[i]  = new TH2F(Form("hYPrimAcc_MC%s",ppname[i].Data()),
2264                                  Form("primary photon %s in acceptance: Rapidity ",pptype[i].Data()),
2265                                  nptbins,ptmin,ptmax,100,-1,1);
2266       fhYPrimMCAcc[i]->SetYTitle("Rapidity");
2267       fhYPrimMCAcc[i]->SetXTitle("E (GeV)");
2268       outputContainer->Add(fhYPrimMCAcc[i]) ;
2269
2270       fhPhiPrimMCAcc[i]  = new TH2F(Form("hPhiPrimAcc_MC%s",ppname[i].Data()),
2271                                  Form("primary photon %s in acceptance: #phi ",pptype[i].Data()),
2272                                  nptbins,ptmin,ptmax,nphibins,phimin,phimax);
2273       fhPhiPrimMCAcc[i]->SetYTitle("#phi (rad)");
2274       fhPhiPrimMCAcc[i]->SetXTitle("E (GeV)");
2275       outputContainer->Add(fhPhiPrimMCAcc[i]) ;
2276
2277     }
2278
2279     if(fFillSSHistograms)
2280     {
2281       TString ptypess[] = { "#gamma","hadron?","#pi^{0}","#eta","#gamma->e^{#pm}","e^{#pm}"} ;
2282
2283       TString pnamess[] = { "Photon","Hadron","Pi0","Eta","Conversion","Electron"} ;
2284
2285       for(Int_t i = 0; i < 6; i++)
2286       {
2287         fhMCELambda0[i]  = new TH2F(Form("hELambda0_MC%s",pnamess[i].Data()),
2288                                     Form("cluster from %s : E vs #lambda_{0}^{2}",ptypess[i].Data()),
2289                                     nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2290         fhMCELambda0[i]->SetYTitle("#lambda_{0}^{2}");
2291         fhMCELambda0[i]->SetXTitle("E (GeV)");
2292         outputContainer->Add(fhMCELambda0[i]) ;
2293
2294         fhMCELambda1[i]  = new TH2F(Form("hELambda1_MC%s",pnamess[i].Data()),
2295                                     Form("cluster from %s : E vs #lambda_{1}^{2}",ptypess[i].Data()),
2296                                     nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2297         fhMCELambda1[i]->SetYTitle("#lambda_{1}^{2}");
2298         fhMCELambda1[i]->SetXTitle("E (GeV)");
2299         outputContainer->Add(fhMCELambda1[i]) ;
2300
2301         fhMCEDispersion[i]  = new TH2F(Form("hEDispersion_MC%s",pnamess[i].Data()),
2302                                        Form("cluster from %s : E vs dispersion^{2}",ptypess[i].Data()),
2303                                        nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2304         fhMCEDispersion[i]->SetYTitle("D^{2}");
2305         fhMCEDispersion[i]->SetXTitle("E (GeV)");
2306         outputContainer->Add(fhMCEDispersion[i]) ;
2307
2308         fhMCNCellsE[i]  = new TH2F (Form("hNCellsE_MC%s",pnamess[i].Data()),
2309                                     Form("# of cells in cluster from %s vs E of clusters",ptypess[i].Data()),
2310                                     nptbins,ptmin,ptmax, nbins,nmin,nmax);
2311         fhMCNCellsE[i]->SetXTitle("E (GeV)");
2312         fhMCNCellsE[i]->SetYTitle("# of cells in cluster");
2313         outputContainer->Add(fhMCNCellsE[i]);
2314
2315         fhMCMaxCellDiffClusterE[i]  = new TH2F (Form("hMaxCellDiffClusterE_MC%s",pnamess[i].Data()),
2316                                                 Form("energy vs difference of cluster energy from %s - max cell energy / cluster energy, good clusters",ptypess[i].Data()),
2317                                                 nptbins,ptmin,ptmax, 500,0,1.);
2318         fhMCMaxCellDiffClusterE[i]->SetXTitle("E_{cluster} (GeV) ");
2319         fhMCMaxCellDiffClusterE[i]->SetYTitle("(E_{cluster} - E_{cell max})/ E_{cluster}");
2320         outputContainer->Add(fhMCMaxCellDiffClusterE[i]);
2321
2322         if(!fFillOnlySimpleSSHisto)
2323         {
2324           fhMCLambda0vsClusterMaxCellDiffE0[i]  = new TH2F(Form("hLambda0vsClusterMaxCellDiffE0_MC%s",pnamess[i].Data()),
2325                                                            Form("cluster from %s : #lambda^{2}_{0} vs fraction of energy carried by max cell, E < 2 GeV",ptypess[i].Data()),
2326                                                            ssbins,ssmin,ssmax,500,0,1.);
2327           fhMCLambda0vsClusterMaxCellDiffE0[i]->SetXTitle("#lambda_{0}^{2}");
2328           fhMCLambda0vsClusterMaxCellDiffE0[i]->SetYTitle("(E_{cluster} - E_{cell max})/ E_{cluster}");
2329           outputContainer->Add(fhMCLambda0vsClusterMaxCellDiffE0[i]) ;
2330
2331           fhMCLambda0vsClusterMaxCellDiffE2[i]  = new TH2F(Form("hLambda0vsClusterMaxCellDiffE2_MC%s",pnamess[i].Data()),
2332                                                            Form("cluster from %s : #lambda^{2}_{0} vs fraction of energy carried by max cell, 2< E < 6 GeV",ptypess[i].Data()),
2333                                                            ssbins,ssmin,ssmax,500,0,1.);
2334           fhMCLambda0vsClusterMaxCellDiffE2[i]->SetXTitle("#lambda_{0}^{2}");
2335           fhMCLambda0vsClusterMaxCellDiffE2[i]->SetYTitle("(E_{cluster} - E_{cell max})/ E_{cluster}");
2336           outputContainer->Add(fhMCLambda0vsClusterMaxCellDiffE2[i]) ;
2337
2338           fhMCLambda0vsClusterMaxCellDiffE6[i]  = new TH2F(Form("hLambda0vsClusterMaxCellDiffE6_MC%s",pnamess[i].Data()),
2339                                                            Form("cluster from %s : #lambda^{2}_{0} vs fraction of energy carried by max cell, E > 6 GeV",ptypess[i].Data()),
2340                                                            ssbins,ssmin,ssmax,500,0,1.);
2341           fhMCLambda0vsClusterMaxCellDiffE6[i]->SetXTitle("#lambda_{0}^{2}");
2342           fhMCLambda0vsClusterMaxCellDiffE6[i]->SetYTitle("(E_{cluster} - E_{cell max})/ E_{cluster}");
2343           outputContainer->Add(fhMCLambda0vsClusterMaxCellDiffE6[i]) ;
2344
2345           fhMCNCellsvsClusterMaxCellDiffE0[i]  = new TH2F(Form("hNCellsvsClusterMaxCellDiffE0_MC%s",pnamess[i].Data()),
2346                                                           Form("cluster from %s : N cells in cluster vs fraction of energy carried by max cell, E < 2 GeV",ptypess[i].Data()),
2347                                                           nbins/5,nmin,nmax/5,500,0,1.);
2348           fhMCNCellsvsClusterMaxCellDiffE0[i]->SetXTitle("N cells in cluster");
2349           fhMCNCellsvsClusterMaxCellDiffE0[i]->SetYTitle("(E_{cluster} - E_{cell max})/ E_{cluster}");
2350           outputContainer->Add(fhMCNCellsvsClusterMaxCellDiffE0[i]) ;
2351
2352           fhMCNCellsvsClusterMaxCellDiffE2[i]  = new TH2F(Form("hNCellsvsClusterMaxCellDiffE2_MC%s",pnamess[i].Data()),
2353                                                           Form("cluster from %s : N cells in cluster vs fraction of energy carried by max cell, 2< E < 6 GeV",ptypess[i].Data()),
2354                                                           nbins/5,nmin,nmax/5,500,0,1.);
2355           fhMCNCellsvsClusterMaxCellDiffE2[i]->SetXTitle("N cells in cluster");
2356           fhMCNCellsvsClusterMaxCellDiffE2[i]->SetYTitle("(E_{cluster} - E_{cell max})/ E_{cluster}");
2357           outputContainer->Add(fhMCNCellsvsClusterMaxCellDiffE2[i]) ;
2358
2359           fhMCNCellsvsClusterMaxCellDiffE6[i]  = new TH2F(Form("hNCellsvsClusterMaxCellDiffE6_MC%s",pnamess[i].Data()),
2360                                                           Form("cluster from %s : N cells in cluster vs fraction of energy carried by max cell, E > 6 GeV",ptypess[i].Data()),
2361                                                           nbins/5,nmin,nmax/5,500,0,1.);
2362           fhMCNCellsvsClusterMaxCellDiffE6[i]->SetXTitle("N cells in cluster");
2363           fhMCNCellsvsClusterMaxCellDiffE6[i]->SetYTitle("E (GeV)");
2364           outputContainer->Add(fhMCNCellsvsClusterMaxCellDiffE6[i]) ;
2365
2366           if(fCalorimeter=="EMCAL")
2367           {
2368             fhMCEDispEta[i]  = new TH2F (Form("hEDispEtaE_MC%s",pnamess[i].Data()),
2369                                          Form("cluster from %s : #sigma^{2}_{#eta #eta} = #Sigma w_{i}(#eta_{i} - <#eta>)^{2}/ #Sigma w_{i} vs E",ptypess[i].Data()),
2370                                          nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
2371             fhMCEDispEta[i]->SetXTitle("E (GeV)");
2372             fhMCEDispEta[i]->SetYTitle("#sigma^{2}_{#eta #eta}");
2373             outputContainer->Add(fhMCEDispEta[i]);
2374
2375             fhMCEDispPhi[i]  = new TH2F (Form("hEDispPhiE_MC%s",pnamess[i].Data()),
2376                                          Form("cluster from %s : #sigma^{2}_{#phi #phi} = #Sigma w_{i}(#phi_{i} - <#phi>)^{2} / #Sigma w_{i} vs E",ptypess[i].Data()),
2377                                          nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
2378             fhMCEDispPhi[i]->SetXTitle("E (GeV)");
2379             fhMCEDispPhi[i]->SetYTitle("#sigma^{2}_{#phi #phi}");
2380             outputContainer->Add(fhMCEDispPhi[i]);
2381
2382             fhMCESumEtaPhi[i]  = new TH2F (Form("hESumEtaPhiE_MC%s",pnamess[i].Data()),
2383                                            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()),
2384                                            nptbins,ptmin,ptmax, 2*ssbins,-ssmax,ssmax);
2385             fhMCESumEtaPhi[i]->SetXTitle("E (GeV)");
2386             fhMCESumEtaPhi[i]->SetYTitle("#delta^{2}_{#eta #phi}");
2387             outputContainer->Add(fhMCESumEtaPhi[i]);
2388
2389             fhMCEDispEtaPhiDiff[i]  = new TH2F (Form("hEDispEtaPhiDiffE_MC%s",pnamess[i].Data()),
2390                                                 Form("cluster from %s : #sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta} vs E",ptypess[i].Data()),
2391                                                 nptbins,ptmin,ptmax,200,-10,10);
2392             fhMCEDispEtaPhiDiff[i]->SetXTitle("E (GeV)");
2393             fhMCEDispEtaPhiDiff[i]->SetYTitle("#sigma^{2}_{#phi #phi}-#sigma^{2}_{#eta #eta}");
2394             outputContainer->Add(fhMCEDispEtaPhiDiff[i]);
2395
2396             fhMCESphericity[i]  = new TH2F (Form("hESphericity_MC%s",pnamess[i].Data()),
2397                                             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()),
2398                                             nptbins,ptmin,ptmax, 200,-1,1);
2399             fhMCESphericity[i]->SetXTitle("E (GeV)");
2400             fhMCESphericity[i]->SetYTitle("s = (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi})");
2401             outputContainer->Add(fhMCESphericity[i]);
2402
2403             for(Int_t ie = 0; ie < 7; ie++)
2404             {
2405               fhMCDispEtaDispPhi[ie][i] = new TH2F (Form("hMCDispEtaDispPhi_EBin%d_MC%s",ie,pnamess[i].Data()),
2406                                                     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]),
2407                                                     ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
2408               fhMCDispEtaDispPhi[ie][i]->SetXTitle("#sigma^{2}_{#eta #eta}");
2409               fhMCDispEtaDispPhi[ie][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
2410               outputContainer->Add(fhMCDispEtaDispPhi[ie][i]);
2411
2412               fhMCLambda0DispEta[ie][i] = new TH2F (Form("hMCLambda0DispEta_EBin%d_MC%s",ie,pnamess[i].Data()),
2413                                                     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]),
2414                                                     ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
2415               fhMCLambda0DispEta[ie][i]->SetXTitle("#lambda^{2}_{0}");
2416               fhMCLambda0DispEta[ie][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
2417               outputContainer->Add(fhMCLambda0DispEta[ie][i]);
2418
2419               fhMCLambda0DispPhi[ie][i] = new TH2F (Form("hMCLambda0DispPhi_EBin%d_MC%s",ie,pnamess[i].Data()),
2420                                                     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]),
2421                                                     ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
2422               fhMCLambda0DispPhi[ie][i]->SetXTitle("#lambda^{2}_{0}");
2423               fhMCLambda0DispPhi[ie][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
2424               outputContainer->Add(fhMCLambda0DispPhi[ie][i]);
2425             }
2426           }
2427         }
2428       }// loop
2429
2430       if(!GetReader()->IsEmbeddedClusterSelectionOn())
2431       {
2432         fhMCPhotonELambda0NoOverlap  = new TH2F("hELambda0_MCPhoton_NoOverlap",
2433                                                 "cluster from Photon : E vs #lambda_{0}^{2}",
2434                                                 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2435         fhMCPhotonELambda0NoOverlap->SetYTitle("#lambda_{0}^{2}");
2436         fhMCPhotonELambda0NoOverlap->SetXTitle("E (GeV)");
2437         outputContainer->Add(fhMCPhotonELambda0NoOverlap) ;
2438
2439         fhMCPhotonELambda0TwoOverlap  = new TH2F("hELambda0_MCPhoton_TwoOverlap",
2440                                                  "cluster from Photon : E vs #lambda_{0}^{2}",
2441                                                  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2442         fhMCPhotonELambda0TwoOverlap->SetYTitle("#lambda_{0}^{2}");
2443         fhMCPhotonELambda0TwoOverlap->SetXTitle("E (GeV)");
2444         outputContainer->Add(fhMCPhotonELambda0TwoOverlap) ;
2445
2446         fhMCPhotonELambda0NOverlap  = new TH2F("hELambda0_MCPhoton_NOverlap",
2447                                                "cluster from Photon : E vs #lambda_{0}^{2}",
2448                                                nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2449         fhMCPhotonELambda0NOverlap->SetYTitle("#lambda_{0}^{2}");
2450         fhMCPhotonELambda0NOverlap->SetXTitle("E (GeV)");
2451         outputContainer->Add(fhMCPhotonELambda0NOverlap) ;
2452
2453       } //No embedding
2454
2455       //Fill histograms to check shape of embedded clusters
2456       if(GetReader()->IsEmbeddedClusterSelectionOn())
2457       {
2458
2459         fhEmbeddedSignalFractionEnergy  = new TH2F("hEmbeddedSignal_FractionEnergy",
2460                                                    "Energy Fraction of embedded signal versus cluster energy",
2461                                                    nptbins,ptmin,ptmax,100,0.,1.);
2462         fhEmbeddedSignalFractionEnergy->SetYTitle("Fraction");
2463         fhEmbeddedSignalFractionEnergy->SetXTitle("E (GeV)");
2464         outputContainer->Add(fhEmbeddedSignalFractionEnergy) ;
2465
2466         fhEmbedPhotonELambda0FullSignal  = new TH2F("hELambda0_EmbedPhoton_FullSignal",
2467                                                     "cluster from Photon embedded with more than 90% energy in cluster : E vs #lambda_{0}^{2}",
2468                                                     nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2469         fhEmbedPhotonELambda0FullSignal->SetYTitle("#lambda_{0}^{2}");
2470         fhEmbedPhotonELambda0FullSignal->SetXTitle("E (GeV)");
2471         outputContainer->Add(fhEmbedPhotonELambda0FullSignal) ;
2472
2473         fhEmbedPhotonELambda0MostlySignal  = new TH2F("hELambda0_EmbedPhoton_MostlySignal",
2474                                                       "cluster from Photon embedded with 50% to 90% energy in cluster : E vs #lambda_{0}^{2}",
2475                                                       nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2476         fhEmbedPhotonELambda0MostlySignal->SetYTitle("#lambda_{0}^{2}");
2477         fhEmbedPhotonELambda0MostlySignal->SetXTitle("E (GeV)");
2478         outputContainer->Add(fhEmbedPhotonELambda0MostlySignal) ;
2479
2480         fhEmbedPhotonELambda0MostlyBkg  = new TH2F("hELambda0_EmbedPhoton_MostlyBkg",
2481                                                    "cluster from Photon embedded with 10% to 50% energy in cluster : E vs #lambda_{0}^{2}",
2482                                                    nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2483         fhEmbedPhotonELambda0MostlyBkg->SetYTitle("#lambda_{0}^{2}");
2484         fhEmbedPhotonELambda0MostlyBkg->SetXTitle("E (GeV)");
2485         outputContainer->Add(fhEmbedPhotonELambda0MostlyBkg) ;
2486
2487         fhEmbedPhotonELambda0FullBkg  = new TH2F("hELambda0_EmbedPhoton_FullBkg",
2488                                                  "cluster from Photonm embedded with 0% to 10% energy in cluster : E vs #lambda_{0}^{2}",
2489                                                  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2490         fhEmbedPhotonELambda0FullBkg->SetYTitle("#lambda_{0}^{2}");
2491         fhEmbedPhotonELambda0FullBkg->SetXTitle("E (GeV)");
2492         outputContainer->Add(fhEmbedPhotonELambda0FullBkg) ;
2493
2494         fhEmbedPi0ELambda0FullSignal  = new TH2F("hELambda0_EmbedPi0_FullSignal",
2495                                                  "cluster from Pi0 embedded with more than 90% energy in cluster : E vs #lambda_{0}^{2}",
2496                                                  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2497         fhEmbedPi0ELambda0FullSignal->SetYTitle("#lambda_{0}^{2}");
2498         fhEmbedPi0ELambda0FullSignal->SetXTitle("E (GeV)");
2499         outputContainer->Add(fhEmbedPi0ELambda0FullSignal) ;
2500
2501         fhEmbedPi0ELambda0MostlySignal  = new TH2F("hELambda0_EmbedPi0_MostlySignal",
2502                                                    "cluster from Pi0 embedded with 50% to 90% energy in cluster : E vs #lambda_{0}^{2}",
2503                                                    nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2504         fhEmbedPi0ELambda0MostlySignal->SetYTitle("#lambda_{0}^{2}");
2505         fhEmbedPi0ELambda0MostlySignal->SetXTitle("E (GeV)");
2506         outputContainer->Add(fhEmbedPi0ELambda0MostlySignal) ;
2507
2508         fhEmbedPi0ELambda0MostlyBkg  = new TH2F("hELambda0_EmbedPi0_MostlyBkg",
2509                                                 "cluster from Pi0 embedded with 10% to 50% energy in cluster : E vs #lambda_{0}^{2}",
2510                                                 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2511         fhEmbedPi0ELambda0MostlyBkg->SetYTitle("#lambda_{0}^{2}");
2512         fhEmbedPi0ELambda0MostlyBkg->SetXTitle("E (GeV)");
2513         outputContainer->Add(fhEmbedPi0ELambda0MostlyBkg) ;
2514
2515         fhEmbedPi0ELambda0FullBkg  = new TH2F("hELambda0_EmbedPi0_FullBkg",
2516                                               "cluster from Pi0 embedded with 0% to 10% energy in cluster : E vs #lambda_{0}^{2}",
2517                                               nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2518         fhEmbedPi0ELambda0FullBkg->SetYTitle("#lambda_{0}^{2}");
2519         fhEmbedPi0ELambda0FullBkg->SetXTitle("E (GeV)");
2520         outputContainer->Add(fhEmbedPi0ELambda0FullBkg) ;
2521
2522       }// embedded histograms
2523
2524
2525     }// Fill SS MC histograms
2526
2527   }//Histos with MC
2528
2529   return outputContainer ;
2530
2531 }
2532
2533 //_______________________
2534 void AliAnaPhoton::Init()
2535 {
2536
2537   //Init
2538   //Do some checks
2539   if(fCalorimeter == "PHOS" && !GetReader()->IsPHOSSwitchedOn() && NewOutputAOD())
2540   {
2541     printf("AliAnaPhoton::Init() - !!STOP: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!\n");
2542     abort();
2543   }
2544   else  if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn() && NewOutputAOD())
2545   {
2546     printf("AliAnaPhoton::Init() - !!STOP: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!\n");
2547     abort();
2548   }
2549
2550   if(GetReader()->GetDataType() == AliCaloTrackReader::kMC) GetCaloPID()->SwitchOnBayesian();
2551
2552 }
2553
2554 //____________________________________________________________________________
2555 void AliAnaPhoton::InitParameters()
2556 {
2557
2558   //Initialize the parameters of the analysis.
2559   AddToHistogramsName("AnaPhoton_");
2560
2561   fCalorimeter = "EMCAL" ;
2562   fMinDist     = 2.;
2563   fMinDist2    = 4.;
2564   fMinDist3    = 5.;
2565
2566   fTimeCutMin  =-1000000;
2567   fTimeCutMax  = 1000000;
2568   fNCellsCut   = 0;
2569
2570   fRejectTrackMatch       = kTRUE ;
2571
2572 }
2573
2574 //__________________________________________________________________
2575 void  AliAnaPhoton::MakeAnalysisFillAOD()
2576 {
2577   //Do photon analysis and fill aods
2578
2579   //Get the vertex
2580   Double_t v[3] = {0,0,0}; //vertex ;
2581   GetReader()->GetVertex(v);
2582
2583   //Select the Calorimeter of the photon
2584   TObjArray * pl = 0x0;
2585   AliVCaloCells* cells    = 0;
2586   if      (fCalorimeter == "PHOS" )
2587   {
2588     pl    = GetPHOSClusters();
2589     cells = GetPHOSCells();
2590   }
2591   else if (fCalorimeter == "EMCAL")
2592   {
2593     pl    = GetEMCALClusters();
2594     cells = GetEMCALCells();
2595   }
2596
2597   if(!pl)
2598   {
2599     Info("MakeAnalysisFillAOD","TObjArray with %s clusters is NULL!\n",fCalorimeter.Data());
2600     return;
2601   }
2602
2603   FillPileUpHistogramsPerEvent(pl);
2604
2605   // Loop on raw clusters before filtering in the reader and fill control histogram
2606   if((GetReader()->GetEMCALClusterListName()=="" && fCalorimeter=="EMCAL") || fCalorimeter=="PHOS")
2607   {
2608     for(Int_t iclus = 0; iclus < GetReader()->GetInputEvent()->GetNumberOfCaloClusters(); iclus++ )
2609     {
2610       AliVCluster * clus = GetReader()->GetInputEvent()->GetCaloCluster(iclus);
2611       if     (fCalorimeter == "PHOS"  && clus->IsPHOS()  && clus->E() > GetReader()->GetPHOSPtMin() ) fhClusterCuts[0]->Fill(clus->E());
2612       else if(fCalorimeter == "EMCAL" && clus->IsEMCAL() && clus->E() > GetReader()->GetEMCALPtMin()) fhClusterCuts[0]->Fill(clus->E());
2613     }
2614   }
2615   else
2616   { // reclusterized
2617     TClonesArray * clusterList = 0;
2618
2619     if(GetReader()->GetInputEvent()->FindListObject(GetReader()->GetEMCALClusterListName()))
2620       clusterList = dynamic_cast<TClonesArray*> (GetReader()->GetInputEvent()->FindListObject(GetReader()->GetEMCALClusterListName()));
2621     else if(GetReader()->GetOutputEvent())
2622       clusterList = dynamic_cast<TClonesArray*> (GetReader()->GetOutputEvent()->FindListObject(GetReader()->GetEMCALClusterListName()));
2623
2624     if(clusterList)
2625     {
2626       Int_t nclusters = clusterList->GetEntriesFast();
2627       for (Int_t iclus =  0; iclus <  nclusters; iclus++)
2628       {
2629         AliVCluster * clus = dynamic_cast<AliVCluster*> (clusterList->At(iclus));
2630         if(clus)fhClusterCuts[0]->Fill(clus->E());
2631       }
2632     }
2633   }
2634
2635   //Init arrays, variables, get number of clusters
2636   TLorentzVector mom, mom2 ;
2637   Int_t nCaloClusters = pl->GetEntriesFast();
2638
2639   if(GetDebug() > 0) printf("AliAnaPhoton::MakeAnalysisFillAOD() - input %s cluster entries %d\n", fCalorimeter.Data(), nCaloClusters);
2640
2641   //----------------------------------------------------
2642   // Fill AOD with PHOS/EMCAL AliAODPWG4Particle objects
2643   //----------------------------------------------------
2644   // Loop on clusters
2645   for(Int_t icalo = 0; icalo < nCaloClusters; icalo++)
2646   {
2647           AliVCluster * calo =  (AliVCluster*) (pl->At(icalo));
2648     //printf("calo %d, %f\n",icalo,calo->E());
2649
2650     //Get the index where the cluster comes, to retrieve the corresponding vertex
2651     Int_t evtIndex = 0 ;
2652     if (GetMixedEvent())
2653     {
2654       evtIndex=GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ;
2655       //Get the vertex and check it is not too large in z
2656       if(TMath::Abs(GetVertex(evtIndex)[2])> GetZvertexCut()) continue;
2657     }
2658
2659     //Cluster selection, not charged, with photon id and in fiducial cut
2660     if(GetReader()->GetDataType() != AliCaloTrackReader::kMC)
2661     {
2662       calo->GetMomentum(mom,GetVertex(evtIndex)) ;}//Assume that come from vertex in straight line
2663     else
2664     {
2665       Double_t vertex[]={0,0,0};
2666       calo->GetMomentum(mom,vertex) ;
2667     }
2668
2669     //--------------------------------------
2670     // Cluster selection
2671     //--------------------------------------
2672     Int_t nMaxima = GetCaloUtils()->GetNumberOfLocalMaxima(calo, cells); // NLM
2673     if(!ClusterSelected(calo,mom,nMaxima)) continue;
2674
2675     //----------------------------
2676     //Create AOD for analysis
2677     //----------------------------
2678     AliAODPWG4Particle aodph = AliAODPWG4Particle(mom);
2679
2680     //...............................................
2681     //Set the indeces of the original caloclusters (MC, ID), and calorimeter
2682     Int_t label = calo->GetLabel();
2683     aodph.SetLabel(label);
2684     aodph.SetCaloLabel(calo->GetID(),-1);
2685     aodph.SetDetector(fCalorimeter);
2686     //printf("Index %d, Id %d, iaod %d\n",icalo, calo->GetID(),GetOutputAODBranch()->GetEntriesFast());
2687
2688     //...............................................
2689     //Set bad channel distance bit
2690     Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
2691     if     (distBad > fMinDist3) aodph.SetDistToBad(2) ;
2692     else if(distBad > fMinDist2) aodph.SetDistToBad(1) ;
2693     else                         aodph.SetDistToBad(0) ;
2694     //printf("DistBad %f Bit %d\n",distBad, aodph.DistToBad());
2695
2696     //--------------------------------------------------------------------------------------
2697     // Play with the MC stack if available
2698     //--------------------------------------------------------------------------------------
2699
2700     //Check origin of the candidates
2701     Int_t tag = -1;
2702
2703     if(IsDataMC())
2704     {
2705       tag = GetMCAnalysisUtils()->CheckOrigin(calo->GetLabels(),calo->GetNLabels(),GetReader(), aodph.GetInputFileIndex());
2706       aodph.SetTag(tag);
2707
2708       if(GetDebug() > 0)
2709         printf("AliAnaPhoton::MakeAnalysisFillAOD() - Origin of candidate, bit map %d\n",aodph.GetTag());
2710     }//Work with stack also
2711
2712
2713     //--------------------------------------------------------------------------------------
2714     //Fill some shower shape histograms before PID is applied
2715     //--------------------------------------------------------------------------------------
2716
2717     FillShowerShapeHistograms(calo,tag);
2718
2719     //-------------------------------------
2720     //PID selection or bit setting
2721     //-------------------------------------
2722
2723     //...............................................
2724     // Data, PID check on
2725     if(IsCaloPIDOn())
2726     {
2727       // Get most probable PID, 2 options check bayesian PID weights or redo PID
2728       // By default, redo PID
2729
2730       aodph.SetIdentifiedParticleType(GetCaloPID()->GetIdentifiedParticleType(calo));
2731
2732       if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PDG of identified particle %d\n",aodph.GetIdentifiedParticleType());
2733
2734       //If cluster does not pass pid, not photon, skip it.
2735       if(aodph.GetIdentifiedParticleType() != AliCaloPID::kPhoton) continue ;
2736
2737     }
2738
2739     //...............................................
2740     // Data, PID check off
2741     else
2742     {
2743       //Set PID bits for later selection (AliAnaPi0 for example)
2744       //GetIdentifiedParticleType already called in SetPIDBits.
2745
2746       GetCaloPID()->SetPIDBits(calo,&aodph, GetCaloUtils(),GetReader()->GetInputEvent());
2747
2748       if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PID Bits set \n");
2749     }
2750
2751     if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - Photon selection cuts passed: pT %3.2f, pdg %d\n",
2752                               aodph.Pt(), aodph.GetIdentifiedParticleType());
2753
2754     fhClusterCuts[9]->Fill(calo->E());
2755
2756     fhNLocMax->Fill(calo->E(),nMaxima);
2757
2758     // Matching after cuts
2759     if(fFillTMHisto) FillTrackMatchingResidualHistograms(calo,1);
2760
2761     // Fill histograms to undertand pile-up before other cuts applied
2762     // Remember to relax time cuts in the reader
2763     FillPileUpHistograms(calo->E(),mom.Pt(),calo->GetTOF()*1e9);
2764
2765     // Add number of local maxima to AOD, method name in AOD to be FIXED
2766     aodph.SetFiducialArea(nMaxima);
2767
2768
2769     //Add AOD with photon object to aod branch
2770     AddAODParticle(aodph);
2771
2772   }//loop
2773
2774   if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD()  End fill AODs, with %d entries \n",GetOutputAODBranch()->GetEntriesFast());
2775
2776 }
2777
2778 //__________________________________________________________________
2779 void  AliAnaPhoton::MakeAnalysisFillHistograms()
2780 {
2781   //Fill histograms
2782
2783   // Get vertex
2784   Double_t v[3] = {0,0,0}; //vertex ;
2785   GetReader()->GetVertex(v);
2786   //fhVertex->Fill(v[0],v[1],v[2]);
2787   if(TMath::Abs(v[2]) > GetZvertexCut()) return ; // done elsewhere for Single Event analysis, but there for mixed event
2788
2789   //----------------------------------
2790   //Loop on stored AOD photons
2791   Int_t naod = GetOutputAODBranch()->GetEntriesFast();
2792   if(GetDebug() > 0) printf("AliAnaPhoton::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod);
2793
2794   for(Int_t iaod = 0; iaod < naod ; iaod++)
2795   {
2796     AliAODPWG4Particle* ph =  (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
2797     Int_t pdg = ph->GetIdentifiedParticleType();
2798
2799     if(GetDebug() > 3)
2800       printf("AliAnaPhoton::MakeAnalysisFillHistograms() - PDG %d, MC TAG %d, Calorimeter %s\n",
2801              ph->GetIdentifiedParticleType(),ph->GetTag(), (ph->GetDetector()).Data()) ;
2802
2803     //If PID used, fill histos with photons in Calorimeter fCalorimeter
2804     if(IsCaloPIDOn() && pdg != AliCaloPID::kPhoton) continue;
2805     if(ph->GetDetector() != fCalorimeter) continue;
2806
2807     if(GetDebug() > 2)
2808       printf("AliAnaPhoton::MakeAnalysisFillHistograms() - ID Photon: pt %f, phi %f, eta %f\n", ph->Pt(),ph->Phi(),ph->Eta()) ;
2809
2810     //................................
2811     //Fill photon histograms
2812     Float_t ptcluster  = ph->Pt();
2813     Float_t phicluster = ph->Phi();
2814     Float_t etacluster = ph->Eta();
2815     Float_t ecluster   = ph->E();
2816
2817     fhEPhoton   ->Fill(ecluster);
2818     fhPtPhoton  ->Fill(ptcluster);
2819     fhPhiPhoton ->Fill(ptcluster,phicluster);
2820     fhEtaPhoton ->Fill(ptcluster,etacluster);
2821     if     (ecluster   > 0.5) fhEtaPhiPhoton  ->Fill(etacluster, phicluster);
2822     else if(GetMinPt() < 0.5) fhEtaPhi05Photon->Fill(etacluster, phicluster);
2823
2824     //Get original cluster, to recover some information
2825     Int_t absID             = 0;
2826     Float_t maxCellFraction = 0;
2827     AliVCaloCells* cells    = 0;
2828     TObjArray * clusters    = 0;
2829     if(fCalorimeter == "EMCAL")
2830     {
2831       cells    = GetEMCALCells();
2832       clusters = GetEMCALClusters();
2833     }
2834     else
2835     {
2836       cells    = GetPHOSCells();
2837       clusters = GetPHOSClusters();
2838     }
2839
2840     Int_t iclus = -1;
2841     AliVCluster *cluster = FindCluster(clusters,ph->GetCaloLabel(0),iclus);
2842     if(cluster)
2843     {
2844       absID = GetCaloUtils()->GetMaxEnergyCell(cells, cluster,maxCellFraction);
2845
2846       // Control histograms
2847       fhMaxCellDiffClusterE->Fill(ph->E(),maxCellFraction);
2848       fhNCellsE            ->Fill(ph->E(),cluster->GetNCells());
2849       fhTimeE              ->Fill(ph->E(),cluster->GetTOF()*1.e9);
2850       if(cells)
2851       {
2852       for(Int_t icell = 0; icell <  cluster->GetNCells(); icell++)
2853         fhCellsE->Fill(ph->E(),cells->GetCellAmplitude(cluster->GetCellsAbsId()[icell]));
2854       }
2855     }
2856
2857     //.......................................
2858     //Play with the MC data if available
2859     if(IsDataMC())
2860     {
2861       if(GetDebug()>0)
2862       {
2863         if(GetReader()->ReadStack() && !GetMCStack())
2864         {
2865           printf("AliAnaPhoton::MakeAnalysisFillHistograms() - Stack not available, is the MC handler called?\n");
2866         }
2867         else if(GetReader()->ReadAODMCParticles() && !GetReader()->GetAODMCParticles(0))
2868         {
2869           printf("AliAnaPhoton::MakeAnalysisFillHistograms() -  Standard MCParticles not available!\n");
2870         }
2871       }
2872
2873       FillAcceptanceHistograms();
2874
2875       //....................................................................
2876       // Access MC information in stack if requested, check that it exists.
2877       Int_t label =ph->GetLabel();
2878
2879       if(label < 0)
2880       {
2881         if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** bad label ***:  label %d \n", label);
2882         continue;
2883       }
2884
2885       Float_t eprim   = 0;
2886       Float_t ptprim  = 0;
2887       Bool_t ok = kFALSE;
2888       TLorentzVector primary = GetMCAnalysisUtils()->GetMother(label,GetReader(),ok);
2889       if(ok)
2890       {
2891         eprim   = primary.Energy();
2892         ptprim  = primary.Pt();
2893       }
2894
2895       Int_t tag =ph->GetTag();
2896       Int_t mcParticleTag = -1;
2897       if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) && fhMCE[kmcPhoton])
2898       {
2899         fhMCE  [kmcPhoton] ->Fill(ecluster);
2900         fhMCPt [kmcPhoton] ->Fill(ptcluster);
2901         fhMCPhi[kmcPhoton] ->Fill(ecluster,phicluster);
2902         fhMCEta[kmcPhoton] ->Fill(ecluster,etacluster);
2903
2904         fhMC2E     [kmcPhoton] ->Fill(ecluster, eprim);
2905         fhMC2Pt    [kmcPhoton] ->Fill(ptcluster, ptprim);
2906         fhMCDeltaE [kmcPhoton] ->Fill(ecluster,eprim-ecluster);
2907         fhMCDeltaPt[kmcPhoton] ->Fill(ptcluster,ptprim-ptcluster);
2908
2909         if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion) &&
2910            GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton)     &&
2911            fhMCE[kmcConversion])
2912         {
2913           fhMCE  [kmcConversion] ->Fill(ecluster);
2914           fhMCPt [kmcConversion] ->Fill(ptcluster);
2915           fhMCPhi[kmcConversion] ->Fill(ecluster,phicluster);
2916           fhMCEta[kmcConversion] ->Fill(ecluster,etacluster);
2917
2918           fhMC2E     [kmcConversion] ->Fill(ecluster, eprim);
2919           fhMC2Pt    [kmcConversion] ->Fill(ptcluster, ptprim);
2920           fhMCDeltaE [kmcConversion] ->Fill(ecluster,eprim-ecluster);
2921           fhMCDeltaPt[kmcConversion] ->Fill(ptcluster,ptprim-ptcluster);
2922         }
2923
2924         if     (GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt) && fhMCE[kmcPrompt])
2925         {
2926           mcParticleTag = kmcPrompt;
2927         }
2928         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation)&& fhMCE[kmcFragmentation])
2929         {
2930           mcParticleTag = kmcFragmentation;
2931         }
2932         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCISR)&& fhMCE[kmcISR])
2933         {
2934           mcParticleTag = kmcISR;
2935         }
2936         else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay) &&
2937                 !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) && fhMCE[kmcPi0Decay])
2938         {
2939           mcParticleTag = kmcPi0Decay;
2940         }
2941         else if((( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay) &&
2942                   !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta)        ) ||
2943                    GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay) ) && fhMCE[kmcOtherDecay])
2944         {
2945           mcParticleTag = kmcOtherDecay;
2946         }
2947         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) && fhMCE[kmcPi0])
2948         {
2949           mcParticleTag = kmcPi0;
2950         }
2951         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta) && fhMCE[kmcEta])
2952         {
2953           mcParticleTag = kmcEta;
2954         }
2955       }
2956       else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiNeutron) && fhMCE[kmcAntiNeutron])
2957       {
2958         mcParticleTag = kmcAntiNeutron;
2959       }
2960       else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiProton) && fhMCE[kmcAntiProton])
2961       {
2962         mcParticleTag = kmcAntiProton;
2963       }
2964       else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron) && fhMCE[kmcElectron])
2965       {
2966         mcParticleTag = kmcElectron;
2967       }
2968       else if( fhMCE[kmcOther])
2969       {
2970         mcParticleTag = kmcOther;
2971
2972         //               printf(" AliAnaPhoton::MakeAnalysisFillHistograms() - Label %d, pT %2.3f Unknown, bits set: ",
2973         //                                      ph->GetLabel(),ph->Pt());
2974         //                for(Int_t i = 0; i < 20; i++) {
2975         //                        if(GetMCAnalysisUtils()->CheckTagBit(tag,i)) printf(" %d, ",i);
2976         //                }
2977         //                printf("\n");
2978
2979       }
2980
2981       fhMCE  [mcParticleTag] ->Fill(ecluster);
2982       fhMCPt [mcParticleTag] ->Fill(ptcluster);
2983       fhMCPhi[mcParticleTag] ->Fill(ecluster,phicluster);
2984       fhMCEta[mcParticleTag] ->Fill(ecluster,etacluster);
2985
2986       fhMC2E[mcParticleTag]     ->Fill(ecluster, eprim);
2987       fhMC2Pt[mcParticleTag]    ->Fill(ptcluster, ptprim);
2988       fhMCDeltaE[mcParticleTag] ->Fill(ecluster,eprim-ecluster);
2989       fhMCDeltaPt[mcParticleTag]->Fill(ptcluster,ptprim-ptcluster);
2990
2991     }//Histograms with MC
2992
2993   }// aod loop
2994
2995 }
2996
2997
2998 //__________________________________________________________________
2999 void AliAnaPhoton::Print(const Option_t * opt) const
3000 {
3001   //Print some relevant parameters set for the analysis
3002
3003   if(! opt)
3004     return;
3005
3006   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
3007   AliAnaCaloTrackCorrBaseClass::Print(" ");
3008
3009   printf("Calorimeter            =     %s\n", fCalorimeter.Data()) ;
3010   printf("Min Distance to Bad Channel   = %2.1f\n",fMinDist);
3011   printf("Min Distance to Bad Channel 2 = %2.1f\n",fMinDist2);
3012   printf("Min Distance to Bad Channel 3 = %2.1f\n",fMinDist3);
3013   printf("Reject clusters with a track matched = %d\n",fRejectTrackMatch);
3014   printf("Time Cut: %3.1f < TOF  < %3.1f\n", fTimeCutMin, fTimeCutMax);
3015   printf("Number of cells in cluster is        > %d \n", fNCellsCut);
3016   printf("    \n") ;
3017
3018 }