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