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