PWGGA/CaloTrackCorrelations/AliAnaPhoton.cxx

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