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