PWGGA/CaloTrackCorrelations/AliAnaPhoton.cxx

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