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