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