PWGGA/CaloTrackCorrelations/AliAnaPi0EbE.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 is 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 // Class for the analysis of high pT pi0 event by event
  18 // Pi0/Eta identified by one of the following:
  19 //  -Invariant mass of 2 cluster in calorimeter
  20 //  -Shower shape analysis in calorimeter
  21 //  -Invariant mass of one cluster in calorimeter and one photon reconstructed in CTS
  22 //
  23 // -- Author: Gustavo Conesa (LNF-INFN) &  Raphaelle Ichou (SUBATECH)
  24 //////////////////////////////////////////////////////////////////////////////
  25
  26
  27 // --- ROOT system ---
  28 #include <TList.h>
  29 #include <TClonesArray.h>
  30 #include <TObjString.h>
  31
  32 // --- Analysis system ---
  33 #include "AliAnaPi0EbE.h"
  34 #include "AliCaloTrackReader.h"
  35 #include "AliIsolationCut.h"
  36 #include "AliNeutralMesonSelection.h"
  37 #include "AliCaloPID.h"
  38 #include "AliMCAnalysisUtils.h"
  39 #include "AliStack.h"
  40 #include "AliFiducialCut.h"
  41 #include "TParticle.h"
  42 #include "AliVCluster.h"
  43 #include "AliESDEvent.h"
  44 #include "AliAODEvent.h"
  45 #include "AliAODMCParticle.h"
  46
  47 ClassImp(AliAnaPi0EbE)
  48
  49 //____________________________
  50 AliAnaPi0EbE::AliAnaPi0EbE() :
  51 AliAnaCaloTrackCorrBaseClass(),fAnaType(kIMCalo),            fCalorimeter(""),
  52 fMinDist(0.),fMinDist2(0.),    fMinDist3(0.),
  53 fNLMCutMin(-1),                fNLMCutMax(10),
  54 fTimeCutMin(-10000),           fTimeCutMax(10000),
  55 fRejectTrackMatch(kTRUE),
  56 fFillPileUpHistograms(0),
  57 fFillWeightHistograms(kFALSE), fFillTMHisto(0),
  58 fFillSelectClHisto(0),         fFillOnlySimpleSSHisto(1),    fFillEMCALBCHistograms(0),
  59 fInputAODGammaConvName(""),
  60 // Histograms
  61 fhPt(0),                       fhE(0),
  62 fhEEta(0),                     fhEPhi(0),
  63 fhPtEta(0),                    fhPtPhi(0),                   fhEtaPhi(0),
  64 fhEtaPhiEMCALBC0(0),           fhEtaPhiEMCALBC1(0),          fhEtaPhiEMCALBCN(0),
  65 fhTimeTriggerEMCALBC0UMReMatchOpenTime(0),
  66 fhTimeTriggerEMCALBC0UMReMatchCheckNeigh(0),
  67 fhTimeTriggerEMCALBC0UMReMatchBoth(0),
  68 fhPtCentrality(),              fhPtEventPlane(0),
  69 fhPtReject(0),                 fhEReject(0),
  70 fhEEtaReject(0),               fhEPhiReject(0),              fhEtaPhiReject(0),
  71 fhMass(0),                     fhMassPt(0),                  fhMassSplitPt(0),
  72 fhSelectedMass(0),             fhSelectedMassPt(0),          fhSelectedMassSplitPt(0),
  73 fhAsymmetry(0),                fhSelectedAsymmetry(0),
  74 fhSplitE(0),                   fhSplitPt(0),
  75 fhSplitPtEta(0),               fhSplitPtPhi(0),
  76 fhNLocMaxSplitPt(0),
  77 fhPtDecay(0),                  fhEDecay(0),
  78 // Shower shape histos
  79 fhEDispersion(0),              fhELambda0(0),                fhELambda1(0),
  80 fhELambda0NoTRD(0),            fhELambda0FracMaxCellCut(0),
  81 fhEFracMaxCell(0),             fhEFracMaxCellNoTRD(0),
  82 fhENCells(0),                  fhETime(0),                   fhEPairDiffTime(0),
  83 fhDispEtaE(0),                 fhDispPhiE(0),
  84 fhSumEtaE(0),                  fhSumPhiE(0),                 fhSumEtaPhiE(0),
  85 fhDispEtaPhiDiffE(0),          fhSphericityE(0),
  86
  87 // MC histos
  88 fhMCE(),                       fhMCPt(),
  89 fhMCPhi(),                     fhMCEta(),
  90 fhMCEReject(),                 fhMCPtReject(),
  91 fhMCPtCentrality(),
  92 fhMCPi0PtGenRecoFraction(0),   fhMCEtaPtGenRecoFraction(0),
  93 fhMCPi0DecayPt(0),             fhMCPi0DecayPtFraction(0),
  94 fhMCEtaDecayPt(0),             fhMCEtaDecayPtFraction(0),
  95 fhMCOtherDecayPt(0),
  96 fhMassPairMCPi0(0),            fhMassPairMCEta(0),
  97 fhAnglePairMCPi0(0),           fhAnglePairMCEta(0),
  98 // Weight studies
  99 fhECellClusterRatio(0),        fhECellClusterLogRatio(0),
 100 fhEMaxCellClusterRatio(0),     fhEMaxCellClusterLogRatio(0),
 101 fhTrackMatchedDEta(0),         fhTrackMatchedDPhi(0),        fhTrackMatchedDEtaDPhi(0),
 102 fhTrackMatchedDEtaPos(0),      fhTrackMatchedDPhiPos(0),     fhTrackMatchedDEtaDPhiPos(0),
 103 fhTrackMatchedDEtaNeg(0),      fhTrackMatchedDPhiNeg(0),     fhTrackMatchedDEtaDPhiNeg(0),
 104 fhTrackMatchedMCParticleE(0),
 105 fhTrackMatchedMCParticleDEta(0), fhTrackMatchedMCParticleDPhi(0),
 106 fhdEdx(0),                     fhEOverP(0),                 fhEOverPNoTRD(0),
 107 // Number of local maxima in cluster
 108 fhNLocMaxE(0),                 fhNLocMaxPt(0),
 109 // PileUp
 110 fhTimePtNoCut(0),                    fhTimePtSPD(0),           fhTimePtSPDMulti(0),
 111 fhTimeNPileUpVertSPD(0),            fhTimeNPileUpVertTrack(0),
 112 fhTimeNPileUpVertContributors(0),
 113 fhTimePileUpMainVertexZDistance(0), fhTimePileUpMainVertexZDiamond(0),
 114 fhPtNPileUpSPDVtx(0),               fhPtNPileUpTrkVtx(0),
 115 fhPtNPileUpSPDVtxTimeCut(0),        fhPtNPileUpTrkVtxTimeCut(0),
 116 fhPtNPileUpSPDVtxTimeCut2(0),       fhPtNPileUpTrkVtxTimeCut2(0)
 117 {
 118   //default ctor
 119
 120   for(Int_t i = 0; i < 6; i++)
 121   {
 122     fhMCE              [i] = 0;
 123     fhMCPt             [i] = 0;
 124     fhMCNLocMaxPt      [i] = 0;
 125     fhMCPhi            [i] = 0;
 126     fhMCEta            [i] = 0;
 127     fhMCPtCentrality   [i] = 0;
 128
 129     fhMCSplitE         [i] = 0;
 130     fhMCSplitPt        [i] = 0;
 131     fhMCSplitPtPhi     [i] = 0;
 132     fhMCSplitPtEta     [i] = 0;
 133     fhMCNLocMaxSplitPt [i] = 0;
 134
 135     fhEMCLambda0       [i] = 0;
 136     fhEMCLambda0NoTRD  [i] = 0;
 137     fhEMCLambda0FracMaxCellCut[i]= 0;
 138     fhEMCFracMaxCell   [i] = 0;
 139     fhEMCLambda1       [i] = 0;
 140     fhEMCDispersion    [i] = 0;
 141
 142     fhMCEDispEta       [i] = 0;
 143     fhMCEDispPhi       [i] = 0;
 144     fhMCESumEtaPhi     [i] = 0;
 145     fhMCEDispEtaPhiDiff[i] = 0;
 146     fhMCESphericity    [i] = 0;
 147     fhMCEAsymmetry     [i] = 0;
 148
 149     fhMCMassPt             [i]=0;
 150     fhMCMassSplitPt        [i]=0;
 151     fhMCSelectedMassPt     [i]=0;
 152     fhMCSelectedMassSplitPt[i]=0;
 153
 154     for(Int_t j = 0; j < 7; j++)
 155     {
 156       fhMCLambda0DispEta    [j][i] = 0;
 157       fhMCLambda0DispPhi    [j][i] = 0;
 158       fhMCDispEtaDispPhi    [j][i] = 0;
 159       fhMCAsymmetryLambda0  [j][i] = 0;
 160       fhMCAsymmetryDispEta  [j][i] = 0;
 161       fhMCAsymmetryDispPhi  [j][i] = 0;
 162     }
 163   }
 164
 165   for(Int_t j = 0; j < 7; j++)
 166   {
 167     fhLambda0DispEta    [j] = 0;
 168     fhLambda0DispPhi    [j] = 0;
 169     fhDispEtaDispPhi    [j] = 0;
 170     fhAsymmetryLambda0  [j] = 0;
 171     fhAsymmetryDispEta  [j] = 0;
 172     fhAsymmetryDispPhi  [j] = 0;
 173
 174     fhPtPileUp       [j] = 0;
 175   }
 176
 177   for(Int_t i = 0; i < 3; i++)
 178   {
 179     fhELambda0LocMax       [i] = 0;
 180     fhELambda1LocMax       [i] = 0;
 181     fhEDispersionLocMax    [i] = 0;
 182     fhEDispEtaLocMax       [i] = 0;
 183     fhEDispPhiLocMax       [i] = 0;
 184     fhESumEtaPhiLocMax     [i] = 0;
 185     fhEDispEtaPhiDiffLocMax[i] = 0;
 186     fhESphericityLocMax    [i] = 0;
 187     fhEAsymmetryLocMax     [i] = 0;
 188   }
 189
 190   //Weight studies
 191   for(Int_t i =0; i < 14; i++){
 192     fhLambda0ForW0[i] = 0;
 193     //fhLambda1ForW0[i] = 0;
 194     if(i<8)fhMassPairLocMax[i] = 0;
 195   }
 196
 197   for(Int_t i = 0; i < 11; i++)
 198   {
 199     fhEtaPhiTriggerEMCALBC       [i] = 0 ;
 200     fhTimeTriggerEMCALBC         [i] = 0 ;
 201     fhTimeTriggerEMCALBCPileUpSPD[i] = 0 ;
 202
 203     fhEtaPhiTriggerEMCALBCUM     [i] = 0 ;
 204     fhTimeTriggerEMCALBCUM       [i] = 0 ;
 205
 206   }
 207
 208   //Initialize parameters
 209   InitParameters();
 210
 211 }
 212
 213 //_______________________________________________________________________________________________
 214 void AliAnaPi0EbE::FillPileUpHistograms(const Float_t pt, const Float_t time, AliVCluster * calo)
 215 {
 216   // Fill some histograms to understand pile-up
 217   if(!fFillPileUpHistograms) return;
 218
 219   //printf("E %f, time %f\n",energy,time);
 220   AliVEvent * event = GetReader()->GetInputEvent();
 221
 222   fhTimePtNoCut->Fill(pt,time);
 223   if(GetReader()->IsPileUpFromSPD())
 224
 225   if(GetReader()->IsPileUpFromSPD())             { fhPtPileUp[0]->Fill(pt); fhTimePtSPD     ->Fill(pt,time); }
 226   if(GetReader()->IsPileUpFromEMCal())             fhPtPileUp[1]->Fill(pt);
 227   if(GetReader()->IsPileUpFromSPDOrEMCal())        fhPtPileUp[2]->Fill(pt);
 228   if(GetReader()->IsPileUpFromSPDAndEMCal())       fhPtPileUp[3]->Fill(pt);
 229   if(GetReader()->IsPileUpFromSPDAndNotEMCal())    fhPtPileUp[4]->Fill(pt);
 230   if(GetReader()->IsPileUpFromEMCalAndNotSPD())    fhPtPileUp[5]->Fill(pt);
 231   if(GetReader()->IsPileUpFromNotSPDAndNotEMCal()) fhPtPileUp[6]->Fill(pt);
 232
 233   if(event->IsPileupFromSPDInMultBins()) fhTimePtSPDMulti->Fill(pt,time);
 234
 235   // cells in cluster
 236
 237   AliVCaloCells* cells = 0;
 238   if(fCalorimeter == "EMCAL") cells = GetEMCALCells();
 239   else                        cells = GetPHOSCells();
 240
 241   Float_t maxCellFraction = 0.;
 242   Int_t absIdMax = GetCaloUtils()->GetMaxEnergyCell(cells,calo,maxCellFraction);
 243
 244   Double_t tmax  = cells->GetCellTime(absIdMax);
 245   GetCaloUtils()->RecalibrateCellTime(tmax, fCalorimeter, absIdMax,GetReader()->GetInputEvent()->GetBunchCrossNumber());
 246   tmax*=1.e9;
 247
 248   //Loop on cells inside cluster, max cell must be over 100 MeV and time in BC=0
 249   if(cells->GetCellAmplitude(absIdMax) > 0.1 && TMath::Abs(tmax) < 30)
 250   {
 251     for (Int_t ipos = 0; ipos < calo->GetNCells(); ipos++)
 252     {
 253       Int_t absId  = calo->GetCellsAbsId()[ipos];
 254
 255       if( absId == absIdMax ) continue ;
 256
 257       Double_t timecell  = cells->GetCellTime(absId);
 258       Float_t  amp       = cells->GetCellAmplitude(absId);
 259       Int_t    bc        = GetReader()->GetInputEvent()->GetBunchCrossNumber();
 260       GetCaloUtils()->GetEMCALRecoUtils()->AcceptCalibrateCell(absId,bc,amp,timecell,cells);
 261       timecell*=1e9;
 262
 263       Float_t diff = (tmax-timecell);
 264
 265       if( cells->GetCellAmplitude(absIdMax) < 0.05 ) continue ;
 266
 267       if(GetReader()->IsPileUpFromSPD())
 268       {
 269         fhPtCellTimePileUp[0]->Fill(pt, timecell);
 270         fhPtTimeDiffPileUp[0]->Fill(pt, diff);
 271        }
 272
 273       if(GetReader()->IsPileUpFromEMCal())
 274       {
 275         fhPtCellTimePileUp[1]->Fill(pt, timecell);
 276         fhPtTimeDiffPileUp[1]->Fill(pt, diff);
 277       }
 278
 279       if(GetReader()->IsPileUpFromSPDOrEMCal())
 280       {
 281         fhPtCellTimePileUp[2]->Fill(pt, timecell);
 282         fhPtTimeDiffPileUp[2]->Fill(pt, diff);
 283       }
 284
 285       if(GetReader()->IsPileUpFromSPDAndEMCal())
 286       {
 287         fhPtCellTimePileUp[3]->Fill(pt, timecell);
 288         fhPtTimeDiffPileUp[3]->Fill(pt, diff);
 289       }
 290
 291       if(GetReader()->IsPileUpFromSPDAndNotEMCal())
 292       {
 293         fhPtCellTimePileUp[4]->Fill(pt, timecell);
 294         fhPtTimeDiffPileUp[4]->Fill(pt, diff);
 295       }
 296
 297       if(GetReader()->IsPileUpFromEMCalAndNotSPD())
 298       {
 299         fhPtCellTimePileUp[5]->Fill(pt, timecell);
 300         fhPtTimeDiffPileUp[5]->Fill(pt, diff);
 301       }
 302
 303       if(GetReader()->IsPileUpFromNotSPDAndNotEMCal())
 304       {
 305         fhPtCellTimePileUp[6]->Fill(pt, timecell);
 306         fhPtTimeDiffPileUp[6]->Fill(pt, diff);
 307       }
 308     }//loop
 309   }
 310
 311   if(pt < 8) return; // Fill time figures for high energy clusters not too close to trigger threshold
 312
 313   AliESDEvent* esdEv = dynamic_cast<AliESDEvent*> (event);
 314   AliAODEvent* aodEv = dynamic_cast<AliAODEvent*> (event);
 315
 316   // N pile up vertices
 317   Int_t nVtxSPD = -1;
 318   Int_t nVtxTrk = -1;
 319
 320   if      (esdEv)
 321   {
 322     nVtxSPD = esdEv->GetNumberOfPileupVerticesSPD();
 323     nVtxTrk = esdEv->GetNumberOfPileupVerticesTracks();
 324
 325   }//ESD
 326   else if (aodEv)
 327   {
 328     nVtxSPD = aodEv->GetNumberOfPileupVerticesSPD();
 329     nVtxTrk = aodEv->GetNumberOfPileupVerticesTracks();
 330   }//AOD
 331
 332   fhTimeNPileUpVertSPD  ->Fill(time,nVtxSPD);
 333   fhTimeNPileUpVertTrack->Fill(time,nVtxTrk);
 334
 335         fhPtNPileUpSPDVtx->Fill(pt,nVtxSPD);
 336         fhPtNPileUpTrkVtx->Fill(pt,nVtxTrk);
 337
 338         if(TMath::Abs(time) < 25)
 339         {
 340                 fhPtNPileUpSPDVtxTimeCut ->Fill(pt,nVtxSPD);
 341                 fhPtNPileUpTrkVtxTimeCut ->Fill(pt,nVtxTrk);
 342   }
 343
 344   if(time < 75 && time > -25)
 345   {
 346     fhPtNPileUpSPDVtxTimeCut2->Fill(pt,nVtxSPD);
 347     fhPtNPileUpTrkVtxTimeCut2->Fill(pt,nVtxTrk);
 348   }
 349
 350   //printf("Is SPD %d, Is SPD Multi %d, n spd %d, n track %d\n",
 351   //       GetReader()->IsPileUpFromSPD(),event->IsPileupFromSPDInMultBins(),nVtxSPD,nVtxTracks);
 352
 353   Int_t ncont = -1;
 354   Float_t z1 = -1, z2 = -1;
 355   Float_t diamZ = -1;
 356   for(Int_t iVert=0; iVert<nVtxSPD;iVert++)
 357   {
 358     if      (esdEv)
 359     {
 360       const AliESDVertex* pv=esdEv->GetPileupVertexSPD(iVert);
 361       ncont=pv->GetNContributors();
 362       z1 = esdEv->GetPrimaryVertexSPD()->GetZ();
 363       z2 = pv->GetZ();
 364       diamZ = esdEv->GetDiamondZ();
 365     }//ESD
 366     else if (aodEv)
 367     {
 368       AliAODVertex *pv=aodEv->GetVertex(iVert);
 369       if(pv->GetType()!=AliAODVertex::kPileupSPD) continue;
 370       ncont=pv->GetNContributors();
 371       z1=aodEv->GetPrimaryVertexSPD()->GetZ();
 372       z2=pv->GetZ();
 373       diamZ = aodEv->GetDiamondZ();
 374     }// AOD
 375
 376     Double_t distZ  = TMath::Abs(z2-z1);
 377     diamZ  = TMath::Abs(z2-diamZ);
 378
 379     fhTimeNPileUpVertContributors  ->Fill(time,ncont);
 380     fhTimePileUpMainVertexZDistance->Fill(time,distZ);
 381     fhTimePileUpMainVertexZDiamond ->Fill(time,diamZ);
 382
 383   }// vertex loop
 384 }
 385
 386
 387 //___________________________________________________________________________________________
 388 void AliAnaPi0EbE::FillRejectedClusterHistograms(const TLorentzVector mom, const Int_t mctag)
 389 {
 390   // Fill histograms that do not pass the identification (SS case only)
 391
 392   Float_t ener  = mom.E();
 393   Float_t pt    = mom.Pt();
 394   Float_t phi   = mom.Phi();
 395   if(phi < 0) phi+=TMath::TwoPi();
 396   Float_t eta = mom.Eta();
 397
 398   fhPtReject     ->Fill(pt);
 399   fhEReject      ->Fill(ener);
 400
 401   fhEEtaReject   ->Fill(ener,eta);
 402   fhEPhiReject   ->Fill(ener,phi);
 403   fhEtaPhiReject ->Fill(eta,phi);
 404
 405   if(IsDataMC())
 406   {
 407     Int_t mcIndex = GetMCIndex(mctag);
 408     fhMCEReject  [mcIndex] ->Fill(ener);
 409     fhMCPtReject [mcIndex] ->Fill(pt);
 410   }
 411 }
 412
 413 //_____________________________________________________________________________________
 414 void AliAnaPi0EbE::FillSelectedClusterHistograms(AliVCluster* cluster,
 415                                                  const Int_t nMaxima,
 416                                                  const Int_t tag,
 417                                                  const Float_t asy)
 418 {
 419   // Fill shower shape, timing and other histograms for selected clusters from decay
 420
 421   Float_t e    = cluster->E();
 422   Float_t disp = cluster->GetDispersion()*cluster->GetDispersion();
 423   Float_t l0   = cluster->GetM02();
 424   Float_t l1   = cluster->GetM20();
 425   Int_t   nSM  = GetModuleNumber(cluster);
 426
 427   Int_t ebin = -1;
 428   if      (e < 2 ) ebin = 0;
 429   else if (e < 4 ) ebin = 1;
 430   else if (e < 6 ) ebin = 2;
 431   else if (e < 10) ebin = 3;
 432   else if (e < 15) ebin = 4;
 433   else if (e < 20) ebin = 5;
 434   else             ebin = 6;
 435
 436   Int_t indexMax = -1;
 437   if     (nMaxima==1) indexMax = 0 ;
 438   else if(nMaxima==2) indexMax = 1 ;
 439   else                indexMax = 2 ;
 440
 441
 442   AliVCaloCells * cell = 0x0;
 443   if(fCalorimeter == "PHOS")
 444     cell = GetPHOSCells();
 445   else
 446     cell = GetEMCALCells();
 447
 448   Float_t maxCellFraction = 0;
 449   GetCaloUtils()->GetMaxEnergyCell(cell, cluster, maxCellFraction);
 450   fhEFracMaxCell->Fill(e,maxCellFraction);
 451
 452   FillWeightHistograms(cluster);
 453
 454   fhEDispersion->Fill(e, disp);
 455   fhELambda0   ->Fill(e, l0  );
 456   fhELambda1   ->Fill(e, l1  );
 457
 458   Float_t ll0  = 0., ll1  = 0.;
 459   Float_t dispp= 0., dEta = 0., dPhi    = 0.;
 460   Float_t sEta = 0., sPhi = 0., sEtaPhi = 0.;
 461   if(fCalorimeter == "EMCAL" && !fFillOnlySimpleSSHisto)
 462   {
 463     GetCaloUtils()->GetEMCALRecoUtils()->RecalculateClusterShowerShapeParameters(GetEMCALGeometry(), GetReader()->GetInputEvent()->GetEMCALCells(), cluster,
 464                                                                                  ll0, ll1, dispp, dEta, dPhi, sEta, sPhi, sEtaPhi);
 465
 466     fhDispEtaE        -> Fill(e,dEta);
 467     fhDispPhiE        -> Fill(e,dPhi);
 468     fhSumEtaE         -> Fill(e,sEta);
 469     fhSumPhiE         -> Fill(e,sPhi);
 470     fhSumEtaPhiE      -> Fill(e,sEtaPhi);
 471     fhDispEtaPhiDiffE -> Fill(e,dPhi-dEta);
 472     if(dEta+dPhi>0)fhSphericityE -> Fill(e,(dPhi-dEta)/(dEta+dPhi));
 473
 474     fhDispEtaDispPhi[ebin]->Fill(dEta,dPhi);
 475     fhLambda0DispEta[ebin]->Fill(l0  ,dEta);
 476     fhLambda0DispPhi[ebin]->Fill(l0  ,dPhi);
 477
 478     if (fAnaType==kSSCalo)
 479     {
 480       // Asymmetry histograms
 481       fhAsymmetryLambda0[ebin]->Fill(l0 ,asy);
 482       fhAsymmetryDispEta[ebin]->Fill(dEta,asy);
 483       fhAsymmetryDispPhi[ebin]->Fill(dPhi,asy);
 484     }
 485   }
 486
 487   fhNLocMaxE ->Fill(e ,nMaxima);
 488
 489   fhELambda0LocMax   [indexMax]->Fill(e,l0);
 490   fhELambda1LocMax   [indexMax]->Fill(e,l1);
 491   fhEDispersionLocMax[indexMax]->Fill(e,disp);
 492
 493   if(fCalorimeter=="EMCAL" && !fFillOnlySimpleSSHisto)
 494   {
 495     fhEDispEtaLocMax       [indexMax]-> Fill(e,dEta);
 496     fhEDispPhiLocMax       [indexMax]-> Fill(e,dPhi);
 497     fhESumEtaPhiLocMax     [indexMax]-> Fill(e,sEtaPhi);
 498     fhEDispEtaPhiDiffLocMax[indexMax]-> Fill(e,dPhi-dEta);
 499     if(dEta+dPhi>0)       fhESphericityLocMax[indexMax]->Fill(e,(dPhi-dEta)/(dEta+dPhi));
 500     if(fAnaType==kSSCalo) fhEAsymmetryLocMax [indexMax]->Fill(e  ,asy);
 501
 502   }
 503
 504   if(fCalorimeter=="EMCAL" && nSM < 6)
 505   {
 506     fhELambda0NoTRD->Fill(e, l0  );
 507     fhEFracMaxCellNoTRD->Fill(e,maxCellFraction);
 508   }
 509
 510   if(maxCellFraction < 0.5)
 511     fhELambda0FracMaxCellCut->Fill(e, l0  );
 512
 513   fhETime  ->Fill(e, cluster->GetTOF()*1.e9);
 514   fhENCells->Fill(e, cluster->GetNCells());
 515
 516   // Fill Track matching control histograms
 517   if(fFillTMHisto)
 518   {
 519     Float_t dZ  = cluster->GetTrackDz();
 520     Float_t dR  = cluster->GetTrackDx();
 521
 522     if(cluster->IsEMCAL() && GetCaloUtils()->IsRecalculationOfClusterTrackMatchingOn())
 523     {
 524       dR = 2000., dZ = 2000.;
 525       GetCaloUtils()->GetEMCALRecoUtils()->GetMatchedResiduals(cluster->GetID(),dZ,dR);
 526     }
 527     //printf("Pi0EbE: dPhi %f, dEta %f\n",dR,dZ);
 528
 529     AliVTrack *track = GetCaloUtils()->GetMatchedTrack(cluster, GetReader()->GetInputEvent());
 530
 531     Bool_t positive = kFALSE;
 532     if(track) positive = (track->Charge()>0);
 533
 534     if(fhTrackMatchedDEta && TMath::Abs(dR) < 999)
 535     {
 536       fhTrackMatchedDEta->Fill(e,dZ);
 537       fhTrackMatchedDPhi->Fill(e,dR);
 538       if(e > 0.5) fhTrackMatchedDEtaDPhi->Fill(dZ,dR);
 539
 540       if(track)
 541       {
 542         if(positive)
 543         {
 544           fhTrackMatchedDEtaPos->Fill(cluster->E(),dZ);
 545           fhTrackMatchedDPhiPos->Fill(cluster->E(),dR);
 546           if(cluster->E() > 0.5) fhTrackMatchedDEtaDPhiPos->Fill(dZ,dR);
 547         }
 548         else
 549         {
 550           fhTrackMatchedDEtaNeg->Fill(cluster->E(),dZ);
 551           fhTrackMatchedDPhiNeg->Fill(cluster->E(),dR);
 552           if(cluster->E() > 0.5) fhTrackMatchedDEtaDPhiNeg->Fill(dZ,dR);
 553         }
 554     }
 555     }
 556     // Check dEdx and E/p of matched clusters
 557
 558     if(TMath::Abs(dZ) < 0.05 && TMath::Abs(dR) < 0.05)
 559     {
 560       if(track)
 561       {
 562         Float_t dEdx = track->GetTPCsignal();
 563         fhdEdx->Fill(e, dEdx);
 564
 565         Float_t eOverp = e/track->P();
 566         fhEOverP->Fill(e,  eOverp);
 567
 568         if(fCalorimeter=="EMCAL" && nSM < 6) fhEOverPNoTRD->Fill(e,  eOverp);
 569
 570       }
 571       //else
 572       //  printf("AliAnaPi0EbE::FillSelectedClusterHistograms() - Residual OK but (dR, dZ)= (%2.4f,%2.4f) no track associated WHAT? \n", dR,dZ);
 573
 574       if(IsDataMC())
 575       {
 576         Float_t mctag = -1;
 577         if  ( !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion)  )
 578         {
 579           if       ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0)      ||
 580                     GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta)       ) mctag =  2.5 ;
 581           else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton)    ) mctag =  0.5 ;
 582           else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron)  ) mctag =  1.5 ;
 583           else                                                                                 mctag =  3.5 ;
 584
 585         }
 586         else
 587         {
 588           if       ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0)      ||
 589                     GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta)       ) mctag =  6.5 ;
 590           else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton)    ) mctag =  4.5 ;
 591           else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron)  ) mctag =  5.5 ;
 592           else                                                                                 mctag =  7.5 ;
 593         }
 594
 595         fhTrackMatchedMCParticleE   ->Fill(e , mctag);
 596         fhTrackMatchedMCParticleDEta->Fill(dZ, mctag);
 597         fhTrackMatchedMCParticleDPhi->Fill(dR, mctag);
 598
 599       }  // MC
 600     }
 601   }// Track matching histograms
 602
 603   if(IsDataMC())
 604   {
 605     Int_t mcIndex = GetMCIndex(tag);
 606
 607     fhEMCLambda0[mcIndex]    ->Fill(e, l0);
 608     fhEMCLambda1[mcIndex]    ->Fill(e, l1);
 609     fhEMCDispersion[mcIndex] ->Fill(e, disp);
 610     fhEMCFracMaxCell[mcIndex]->Fill(e,maxCellFraction);
 611
 612     if(fCalorimeter=="EMCAL" && nSM < 6)
 613       fhEMCLambda0NoTRD[mcIndex]->Fill(e, l0  );
 614
 615     if(maxCellFraction < 0.5)
 616       fhEMCLambda0FracMaxCellCut[mcIndex]->Fill(e, l0  );
 617
 618     if(fCalorimeter == "EMCAL" && !fFillOnlySimpleSSHisto)
 619     {
 620       fhMCEDispEta        [mcIndex]-> Fill(e,dEta);
 621       fhMCEDispPhi        [mcIndex]-> Fill(e,dPhi);
 622       fhMCESumEtaPhi      [mcIndex]-> Fill(e,sEtaPhi);
 623       fhMCEDispEtaPhiDiff [mcIndex]-> Fill(e,dPhi-dEta);
 624       if(dEta+dPhi>0)fhMCESphericity[mcIndex]-> Fill(e,(dPhi-dEta)/(dEta+dPhi));
 625
 626       if (fAnaType==kSSCalo)
 627       {
 628         fhMCAsymmetryLambda0[ebin][mcIndex]->Fill(l0 ,asy);
 629         fhMCAsymmetryDispEta[ebin][mcIndex]->Fill(dEta,asy);
 630         fhMCAsymmetryDispPhi[ebin][mcIndex]->Fill(dPhi,asy);
 631       }
 632
 633       fhMCDispEtaDispPhi[ebin][mcIndex]->Fill(dEta,dPhi);
 634       fhMCLambda0DispEta[ebin][mcIndex]->Fill(l0  ,dEta);
 635       fhMCLambda0DispPhi[ebin][mcIndex]->Fill(l0  ,dPhi);
 636
 637     }
 638
 639   }//MC
 640
 641 }
 642
 643 //________________________________________________________
 644 void AliAnaPi0EbE::FillWeightHistograms(AliVCluster *clus)
 645 {
 646   // Calculate weights and fill histograms
 647
 648   if(!fFillWeightHistograms || GetMixedEvent()) return;
 649
 650   AliVCaloCells* cells = 0;
 651   if(fCalorimeter == "EMCAL") cells = GetEMCALCells();
 652   else                        cells = GetPHOSCells();
 653
 654   // First recalculate energy in case non linearity was applied
 655   Float_t  energy = 0;
 656   Float_t  ampMax = 0;
 657   for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++)
 658   {
 659
 660     Int_t id       = clus->GetCellsAbsId()[ipos];
 661
 662     //Recalibrate cell energy if needed
 663     Float_t amp = cells->GetCellAmplitude(id);
 664     GetCaloUtils()->RecalibrateCellAmplitude(amp,fCalorimeter, id);
 665
 666     energy    += amp;
 667
 668     if(amp> ampMax)
 669       ampMax = amp;
 670
 671   } // energy loop
 672
 673   if(energy <=0 )
 674   {
 675     printf("AliAnaPi0EbE::WeightHistograms()- Wrong calculated energy %f\n",energy);
 676     return;
 677   }
 678
 679   fhEMaxCellClusterRatio   ->Fill(energy,ampMax/energy);
 680   fhEMaxCellClusterLogRatio->Fill(energy,TMath::Log(ampMax/energy));
 681
 682   //Get the ratio and log ratio to all cells in cluster
 683   for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++)
 684   {
 685     Int_t id       = clus->GetCellsAbsId()[ipos];
 686
 687     //Recalibrate cell energy if needed
 688     Float_t amp = cells->GetCellAmplitude(id);
 689     GetCaloUtils()->RecalibrateCellAmplitude(amp,fCalorimeter, id);
 690
 691     fhECellClusterRatio   ->Fill(energy,amp/energy);
 692     fhECellClusterLogRatio->Fill(energy,TMath::Log(amp/energy));
 693   }
 694
 695   //Recalculate shower shape for different W0
 696   if(fCalorimeter=="EMCAL"){
 697
 698     Float_t l0org = clus->GetM02();
 699     Float_t l1org = clus->GetM20();
 700     Float_t dorg  = clus->GetDispersion();
 701
 702     for(Int_t iw = 0; iw < 14; iw++)
 703     {
 704       GetCaloUtils()->GetEMCALRecoUtils()->SetW0(1+iw*0.5);
 705       GetCaloUtils()->GetEMCALRecoUtils()->RecalculateClusterShowerShapeParameters(GetEMCALGeometry(), cells, clus);
 706
 707       fhLambda0ForW0[iw]->Fill(energy,clus->GetM02());
 708       //fhLambda1ForW0[iw]->Fill(energy,clus->GetM20());
 709
 710     } // w0 loop
 711
 712     // Set the original values back
 713     clus->SetM02(l0org);
 714     clus->SetM20(l1org);
 715     clus->SetDispersion(dorg);
 716
 717   }// EMCAL
 718 }
 719
 720 //__________________________________________
 721 TObjString * AliAnaPi0EbE::GetAnalysisCuts()
 722 {
 723         //Save parameters used for analysis
 724   TString parList ; //this will be list of parameters used for this analysis.
 725   const Int_t buffersize = 255;
 726   char onePar[buffersize] ;
 727
 728   snprintf(onePar,buffersize,"--- AliAnaPi0EbE ---\n") ;
 729   parList+=onePar ;
 730   snprintf(onePar,buffersize,"fAnaType=%d (Pi0 selection type) \n",fAnaType) ;
 731   parList+=onePar ;
 732
 733   if(fAnaType == kSSCalo)
 734   {
 735     snprintf(onePar,buffersize,"Calorimeter: %s\n",fCalorimeter.Data()) ;
 736     parList+=onePar ;
 737     snprintf(onePar,buffersize,"fMinDist =%2.2f (Minimal distance to bad channel to accept cluster) \n",fMinDist) ;
 738     parList+=onePar ;
 739     snprintf(onePar,buffersize,"fMinDist2=%2.2f (Cuts on Minimal distance to study acceptance evaluation) \n",fMinDist2) ;
 740     parList+=onePar ;
 741     snprintf(onePar,buffersize,"fMinDist3=%2.2f (One more cut on distance used for acceptance-efficiency study) \n",fMinDist3) ;
 742     parList+=onePar ;
 743   }
 744
 745   //Get parameters set in base class.
 746   parList += GetBaseParametersList() ;
 747
 748   //Get parameters set in PID class.
 749   if(fAnaType == kSSCalo) parList += GetCaloPID()->GetPIDParametersList() ;
 750
 751   return new TObjString(parList) ;
 752 }
 753
 754 //_____________________________________________
 755 TList *  AliAnaPi0EbE::GetCreateOutputObjects()
 756 {
 757   // Create histograms to be saved in output file and
 758   // store them in outputContainer
 759   TList * outputContainer = new TList() ;
 760   outputContainer->SetName("Pi0EbEHistos") ;
 761
 762   Int_t nptbins  = GetHistogramRanges()->GetHistoPtBins();           Float_t ptmax  = GetHistogramRanges()->GetHistoPtMax();           Float_t ptmin  = GetHistogramRanges()->GetHistoPtMin();
 763   Int_t nphibins = GetHistogramRanges()->GetHistoPhiBins();          Float_t phimax = GetHistogramRanges()->GetHistoPhiMax();          Float_t phimin = GetHistogramRanges()->GetHistoPhiMin();
 764   Int_t netabins = GetHistogramRanges()->GetHistoEtaBins();          Float_t etamax = GetHistogramRanges()->GetHistoEtaMax();          Float_t etamin = GetHistogramRanges()->GetHistoEtaMin();
 765   Int_t ssbins   = GetHistogramRanges()->GetHistoShowerShapeBins();  Float_t ssmax  = GetHistogramRanges()->GetHistoShowerShapeMax();  Float_t ssmin  = GetHistogramRanges()->GetHistoShowerShapeMin();
 766   Int_t tdbins   = GetHistogramRanges()->GetHistoDiffTimeBins() ;    Float_t tdmax  = GetHistogramRanges()->GetHistoDiffTimeMax();     Float_t tdmin  = GetHistogramRanges()->GetHistoDiffTimeMin();
 767   Int_t tbins    = GetHistogramRanges()->GetHistoTimeBins() ;        Float_t tmax   = GetHistogramRanges()->GetHistoTimeMax();         Float_t tmin   = GetHistogramRanges()->GetHistoTimeMin();
 768   Int_t nbins    = GetHistogramRanges()->GetHistoNClusterCellBins(); Int_t   nmax   = GetHistogramRanges()->GetHistoNClusterCellMax(); Int_t   nmin   = GetHistogramRanges()->GetHistoNClusterCellMin();
 769
 770   Int_t   nmassbins   = GetHistogramRanges()->GetHistoMassBins();
 771   Float_t massmin     = GetHistogramRanges()->GetHistoMassMin();
 772   Float_t massmax     = GetHistogramRanges()->GetHistoMassMax();
 773
 774   Int_t   nresetabins = GetHistogramRanges()->GetHistoTrackResidualEtaBins();
 775   Float_t resetamax   = GetHistogramRanges()->GetHistoTrackResidualEtaMax();
 776   Float_t resetamin   = GetHistogramRanges()->GetHistoTrackResidualEtaMin();
 777   Int_t   nresphibins = GetHistogramRanges()->GetHistoTrackResidualPhiBins();
 778   Float_t resphimax   = GetHistogramRanges()->GetHistoTrackResidualPhiMax();
 779   Float_t resphimin   = GetHistogramRanges()->GetHistoTrackResidualPhiMin();
 780
 781   Int_t   ndedxbins   = GetHistogramRanges()->GetHistodEdxBins();
 782   Float_t dedxmax     = GetHistogramRanges()->GetHistodEdxMax();
 783   Float_t dedxmin     = GetHistogramRanges()->GetHistodEdxMin();
 784   Int_t   nPoverEbins = GetHistogramRanges()->GetHistoPOverEBins();
 785   Float_t pOverEmax   = GetHistogramRanges()->GetHistoPOverEMax();
 786   Float_t pOverEmin   = GetHistogramRanges()->GetHistoPOverEMin();
 787
 788   Int_t   ntimebins= GetHistogramRanges()->GetHistoTimeBins();
 789   Float_t timemax  = GetHistogramRanges()->GetHistoTimeMax();
 790   Float_t timemin  = GetHistogramRanges()->GetHistoTimeMin();
 791
 792   TString nlm[]   ={"1 Local Maxima","2 Local Maxima", "NLM > 2"};
 793   TString ptype[] ={"#gamma","#gamma->e^{#pm}","#pi^{0}","#eta","e^{#pm}", "hadron"};
 794   TString pname[] ={"Photon","Conversion",     "Pi0",    "Eta", "Electron","Hadron"};
 795   Int_t   bin[]   = {0,2,4,6,10,15,20,100}; // energy bins
 796
 797   fhPt  = new TH1F("hPt","Number of identified  #pi^{0} (#eta) decay",nptbins,ptmin,ptmax);
 798   fhPt->SetYTitle("N");
 799   fhPt->SetXTitle("p_{T} (GeV/c)");
 800   outputContainer->Add(fhPt) ;
 801
 802   fhE  = new TH1F("hE","Number of identified  #pi^{0} (#eta) decay pairs",nptbins,ptmin,ptmax);
 803   fhE->SetYTitle("N");
 804   fhE->SetXTitle("E (GeV)");
 805   outputContainer->Add(fhE) ;
 806
 807   fhEPhi  = new TH2F
 808   ("hEPhi","Selected #pi^{0} (#eta) pairs: E vs #phi",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
 809   fhEPhi->SetYTitle("#phi (rad)");
 810   fhEPhi->SetXTitle("E (GeV)");
 811   outputContainer->Add(fhEPhi) ;
 812
 813   fhEEta  = new TH2F
 814   ("hEEta","Selected #pi^{0} (#eta) pairs: E vs #eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 815   fhEEta->SetYTitle("#eta");
 816   fhEEta->SetXTitle("E (GeV)");
 817   outputContainer->Add(fhEEta) ;
 818
 819   fhPtPhi  = new TH2F
 820   ("hPtPhi","Selected #pi^{0} (#eta) pairs: p_{T} vs #phi",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
 821   fhPtPhi->SetYTitle("#phi (rad)");
 822   fhPtPhi->SetXTitle("p_{T} (GeV/c)");
 823   outputContainer->Add(fhPtPhi) ;
 824
 825   fhPtEta  = new TH2F
 826   ("hPtEta","Selected #pi^{0} (#eta) pairs: p_{T} vs #eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 827   fhPtEta->SetYTitle("#eta");
 828   fhPtEta->SetXTitle("p_{T} (GeV/c)");
 829   outputContainer->Add(fhPtEta) ;
 830
 831   fhEtaPhi  = new TH2F
 832   ("hEtaPhi","Selected #pi^{0} (#eta) pairs: #eta vs #phi",netabins,etamin,etamax, nphibins,phimin,phimax);
 833   fhEtaPhi->SetYTitle("#phi (rad)");
 834   fhEtaPhi->SetXTitle("#eta");
 835   outputContainer->Add(fhEtaPhi) ;
 836
 837   if(fCalorimeter=="EMCAL" && fFillEMCALBCHistograms)
 838   {
 839     fhEtaPhiEMCALBC0  = new TH2F
 840     ("hEtaPhiEMCALBC0","cluster,E > 2 GeV, #eta vs #phi, for clusters with |time| < 25 ns, EMCAL-BC=0",netabins,etamin,etamax,nphibins,phimin,phimax);
 841     fhEtaPhiEMCALBC0->SetYTitle("#phi (rad)");
 842     fhEtaPhiEMCALBC0->SetXTitle("#eta");
 843     outputContainer->Add(fhEtaPhiEMCALBC0) ;
 844
 845     fhEtaPhiEMCALBC1  = new TH2F
 846     ("hEtaPhiEMCALBC1","cluster,E > 2 GeV, #eta vs #phi, for clusters with 25 < |time| < 75 ns, EMCAL-BC=1",netabins,etamin,etamax,nphibins,phimin,phimax);
 847     fhEtaPhiEMCALBC1->SetYTitle("#phi (rad)");
 848     fhEtaPhiEMCALBC1->SetXTitle("#eta");
 849     outputContainer->Add(fhEtaPhiEMCALBC1) ;
 850
 851     fhEtaPhiEMCALBCN  = new TH2F
 852     ("hEtaPhiEMCALBCN","cluster,E > 2 GeV, #eta vs #phi, for clusters with |time| > 75 ns, EMCAL-BC>1",netabins,etamin,etamax,nphibins,phimin,phimax);
 853     fhEtaPhiEMCALBCN->SetYTitle("#phi (rad)");
 854     fhEtaPhiEMCALBCN->SetXTitle("#eta");
 855     outputContainer->Add(fhEtaPhiEMCALBCN) ;
 856
 857     for(Int_t i = 0; i < 11; i++)
 858     {
 859       fhEtaPhiTriggerEMCALBC[i] = new TH2F
 860       (Form("hEtaPhiTriggerEMCALBC%d",i-5),
 861        Form("meson E > 2 GeV, #eta vs #phi, Trigger EMCAL-BC=%d",i-5),
 862        netabins,etamin,etamax,nphibins,phimin,phimax);
 863       fhEtaPhiTriggerEMCALBC[i]->SetYTitle("#phi (rad)");
 864       fhEtaPhiTriggerEMCALBC[i]->SetXTitle("#eta");
 865       outputContainer->Add(fhEtaPhiTriggerEMCALBC[i]) ;
 866
 867       fhTimeTriggerEMCALBC[i] = new TH2F
 868       (Form("hTimeTriggerEMCALBC%d",i-5),
 869        Form("meson time vs E, Trigger EMCAL-BC=%d",i-5),
 870        nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
 871       fhTimeTriggerEMCALBC[i]->SetXTitle("E (GeV)");
 872       fhTimeTriggerEMCALBC[i]->SetYTitle("time (ns)");
 873       outputContainer->Add(fhTimeTriggerEMCALBC[i]);
 874
 875       fhTimeTriggerEMCALBCPileUpSPD[i] = new TH2F
 876       (Form("hTimeTriggerEMCALBC%dPileUpSPD",i-5),
 877        Form("meson time vs E, Trigger EMCAL-BC=%d",i-5),
 878        nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
 879       fhTimeTriggerEMCALBCPileUpSPD[i]->SetXTitle("E (GeV)");
 880       fhTimeTriggerEMCALBCPileUpSPD[i]->SetYTitle("time (ns)");
 881       outputContainer->Add(fhTimeTriggerEMCALBCPileUpSPD[i]);
 882
 883       fhEtaPhiTriggerEMCALBCUM[i] = new TH2F
 884       (Form("hEtaPhiTriggerEMCALBC%d_UnMatch",i-5),
 885        Form("meson E > 2 GeV, #eta vs #phi, unmatched trigger EMCAL-BC=%d",i-5),
 886        netabins,etamin,etamax,nphibins,phimin,phimax);
 887       fhEtaPhiTriggerEMCALBCUM[i]->SetYTitle("#phi (rad)");
 888       fhEtaPhiTriggerEMCALBCUM[i]->SetXTitle("#eta");
 889       outputContainer->Add(fhEtaPhiTriggerEMCALBCUM[i]) ;
 890
 891       fhTimeTriggerEMCALBCUM[i] = new TH2F
 892       (Form("hTimeTriggerEMCALBC%d_UnMatch",i-5),
 893        Form("meson time vs E, unmatched trigger EMCAL-BC=%d",i-5),
 894        nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
 895       fhTimeTriggerEMCALBCUM[i]->SetXTitle("E (GeV)");
 896       fhTimeTriggerEMCALBCUM[i]->SetYTitle("time (ns)");
 897       outputContainer->Add(fhTimeTriggerEMCALBCUM[i]);
 898
 899     }
 900
 901     fhTimeTriggerEMCALBC0UMReMatchOpenTime = new TH2F("hTimeTriggerBC0_UnMatch_ReMatch_OpenTime",
 902                                                       "cluster time vs E of clusters, no match, rematch open time",
 903                                                       nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
 904     fhTimeTriggerEMCALBC0UMReMatchOpenTime->SetXTitle("E (GeV)");
 905     fhTimeTriggerEMCALBC0UMReMatchOpenTime->SetYTitle("time (ns)");
 906     outputContainer->Add(fhTimeTriggerEMCALBC0UMReMatchOpenTime);
 907
 908
 909     fhTimeTriggerEMCALBC0UMReMatchCheckNeigh = new TH2F("hTimeTriggerBC0_UnMatch_ReMatch_CheckNeighbours",
 910                                                         "cluster time vs E of clusters, no match, rematch with neigbour parches",
 911                                                         nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
 912     fhTimeTriggerEMCALBC0UMReMatchCheckNeigh->SetXTitle("E (GeV)");
 913     fhTimeTriggerEMCALBC0UMReMatchCheckNeigh->SetYTitle("time (ns)");
 914     outputContainer->Add(fhTimeTriggerEMCALBC0UMReMatchCheckNeigh);
 915
 916     fhTimeTriggerEMCALBC0UMReMatchBoth = new TH2F("hTimeTriggerBC0_UnMatch_ReMatch_Both",
 917                                                   "cluster time vs E of clusters, no match, rematch open time and neigbour",
 918                                                   nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
 919     fhTimeTriggerEMCALBC0UMReMatchBoth->SetXTitle("E (GeV)");
 920     fhTimeTriggerEMCALBC0UMReMatchBoth->SetYTitle("time (ns)");
 921     outputContainer->Add(fhTimeTriggerEMCALBC0UMReMatchBoth);
 922
 923   }
 924
 925   fhPtCentrality  = new TH2F("hPtCentrality","centrality vs p_{T}",nptbins,ptmin,ptmax, 100,0,100);
 926   fhPtCentrality->SetYTitle("centrality");
 927   fhPtCentrality->SetXTitle("p_{T}(GeV/c)");
 928   outputContainer->Add(fhPtCentrality) ;
 929
 930   fhPtEventPlane  = new TH2F("hPtEventPlane","event plane angle vs p_{T}",nptbins,ptmin,ptmax, 100,0,TMath::Pi());
 931   fhPtEventPlane->SetYTitle("Event plane angle (rad)");
 932   fhPtEventPlane->SetXTitle("p_{T} (GeV/c)");
 933   outputContainer->Add(fhPtEventPlane) ;
 934
 935   if(fAnaType == kSSCalo)
 936   {
 937     fhPtReject  = new TH1F("hPtReject","Number of rejected as #pi^{0} (#eta) decay",nptbins,ptmin,ptmax);
 938     fhPtReject->SetYTitle("N");
 939     fhPtReject->SetXTitle("p_{T} (GeV/c)");
 940     outputContainer->Add(fhPtReject) ;
 941
 942     fhEReject  = new TH1F("hEReject","Number of rejected as  #pi^{0} (#eta) decay pairs",nptbins,ptmin,ptmax);
 943     fhEReject->SetYTitle("N");
 944     fhEReject->SetXTitle("E (GeV)");
 945     outputContainer->Add(fhEReject) ;
 946
 947     fhEPhiReject  = new TH2F
 948     ("hEPhiReject","Rejected #pi^{0} (#eta) cluster: E vs #phi",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
 949     fhEPhiReject->SetYTitle("#phi (rad)");
 950     fhEPhiReject->SetXTitle("E (GeV)");
 951     outputContainer->Add(fhEPhiReject) ;
 952
 953     fhEEtaReject  = new TH2F
 954     ("hEEtaReject","Rejected #pi^{0} (#eta) cluster: E vs #eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 955     fhEEtaReject->SetYTitle("#eta");
 956     fhEEtaReject->SetXTitle("E (GeV)");
 957     outputContainer->Add(fhEEtaReject) ;
 958
 959     fhEtaPhiReject  = new TH2F
 960     ("hEtaPhiReject","Rejected #pi^{0} (#eta) cluster: #eta vs #phi",netabins,etamin,etamax, nphibins,phimin,phimax);
 961     fhEtaPhiReject->SetYTitle("#phi (rad)");
 962     fhEtaPhiReject->SetXTitle("#eta");
 963     outputContainer->Add(fhEtaPhiReject) ;
 964   }
 965
 966   fhMass  = new TH2F
 967   ("hMass","all pairs mass: E vs mass",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
 968   fhMass->SetYTitle("mass (GeV/c^{2})");
 969   fhMass->SetXTitle("E (GeV)");
 970   outputContainer->Add(fhMass) ;
 971
 972   fhSelectedMass  = new TH2F
 973   ("hSelectedMass","Selected #pi^{0} (#eta) pairs mass: E vs mass",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
 974   fhSelectedMass->SetYTitle("mass (GeV/c^{2})");
 975   fhSelectedMass->SetXTitle("E (GeV)");
 976   outputContainer->Add(fhSelectedMass) ;
 977
 978   fhMassPt  = new TH2F
 979   ("hMassPt","all pairs mass: p_{T} vs mass",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
 980   fhMassPt->SetYTitle("mass (GeV/c^{2})");
 981   fhMassPt->SetXTitle("p_{T} (GeV/c)");
 982   outputContainer->Add(fhMassPt) ;
 983
 984   fhSelectedMassPt  = new TH2F
 985   ("hSelectedMassPt","Selected #pi^{0} (#eta) pairs mass: p_{T} vs mass",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
 986   fhSelectedMassPt->SetYTitle("mass (GeV/c^{2})");
 987   fhSelectedMassPt->SetXTitle("p_{T} (GeV/c)");
 988   outputContainer->Add(fhSelectedMassPt) ;
 989
 990   if(fAnaType != kSSCalo)
 991   {
 992     fhPtDecay  = new TH1F("hPtDecay","Number of identified  #pi^{0} (#eta) decay photons",nptbins,ptmin,ptmax);
 993     fhPtDecay->SetYTitle("N");
 994     fhPtDecay->SetXTitle("p_{T} (GeV/c)");
 995     outputContainer->Add(fhPtDecay) ;
 996
 997     fhEDecay  = new TH1F("hEDecay","Number of identified  #pi^{0} (#eta) decay photons",nptbins,ptmin,ptmax);
 998     fhEDecay->SetYTitle("N");
 999     fhEDecay->SetXTitle("E (GeV)");
1000     outputContainer->Add(fhEDecay) ;
1001   }
1002
1003   ////////
1004
1005   if( fFillSelectClHisto )
1006   {
1007
1008     fhEDispersion  = new TH2F
1009     ("hEDispersion","Selected #pi^{0} (#eta) pairs: E vs dispersion",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1010     fhEDispersion->SetYTitle("D^{2}");
1011     fhEDispersion->SetXTitle("E (GeV)");
1012     outputContainer->Add(fhEDispersion) ;
1013
1014     fhELambda0  = new TH2F
1015     ("hELambda0","Selected #pi^{0} (#eta) pairs: E vs #lambda_{0}",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1016     fhELambda0->SetYTitle("#lambda_{0}^{2}");
1017     fhELambda0->SetXTitle("E (GeV)");
1018     outputContainer->Add(fhELambda0) ;
1019
1020     fhELambda1  = new TH2F
1021     ("hELambda1","Selected #pi^{0} (#eta) pairs: E vs #lambda_{1}",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1022     fhELambda1->SetYTitle("#lambda_{1}^{2}");
1023     fhELambda1->SetXTitle("E (GeV)");
1024     outputContainer->Add(fhELambda1) ;
1025
1026     fhELambda0FracMaxCellCut  = new TH2F
1027     ("hELambda0FracMaxCellCut","Selected #pi^{0} (#eta) pairs: E vs #lambda_{0}, Max cell fraction of energy < 0.5",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1028     fhELambda0FracMaxCellCut->SetYTitle("#lambda_{0}^{2}");
1029     fhELambda0FracMaxCellCut->SetXTitle("E (GeV)");
1030     outputContainer->Add(fhELambda0FracMaxCellCut) ;
1031
1032     fhEFracMaxCell  = new TH2F
1033     ("hEFracMaxCell","Selected #pi^{0} (#eta) pairs: E vs #lambda_{0}, Max cell fraction of energy",nptbins,ptmin,ptmax,100,0,1);
1034     fhEFracMaxCell->SetYTitle("Fraction");
1035     fhEFracMaxCell->SetXTitle("E (GeV)");
1036     outputContainer->Add(fhEFracMaxCell) ;
1037
1038     if(fCalorimeter=="EMCAL")
1039     {
1040       fhELambda0NoTRD  = new TH2F
1041       ("hELambda0NoTRD","Selected #pi^{0} (#eta) pairs: E vs #lambda_{0}, not behind TRD",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1042       fhELambda0NoTRD->SetYTitle("#lambda_{0}^{2}");
1043       fhELambda0NoTRD->SetXTitle("E (GeV)");
1044       outputContainer->Add(fhELambda0NoTRD) ;
1045
1046       fhEFracMaxCellNoTRD  = new TH2F
1047       ("hEFracMaxCellNoTRD","Selected #pi^{0} (#eta) pairs: E vs #lambda_{0}, Max cell fraction of energy, not behind TRD",nptbins,ptmin,ptmax,100,0,1);
1048       fhEFracMaxCellNoTRD->SetYTitle("Fraction");
1049       fhEFracMaxCellNoTRD->SetXTitle("E (GeV)");
1050       outputContainer->Add(fhEFracMaxCellNoTRD) ;
1051
1052       if(!fFillOnlySimpleSSHisto)
1053       {
1054         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);
1055         fhDispEtaE->SetXTitle("E (GeV)");
1056         fhDispEtaE->SetYTitle("#sigma^{2}_{#eta #eta}");
1057         outputContainer->Add(fhDispEtaE);
1058
1059         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);
1060         fhDispPhiE->SetXTitle("E (GeV)");
1061         fhDispPhiE->SetYTitle("#sigma^{2}_{#phi #phi}");
1062         outputContainer->Add(fhDispPhiE);
1063
1064         fhSumEtaE  = new TH2F ("hSumEtaE","#sigma^{2}_{#eta #eta} = #Sigma w_{i}(#eta_{i})^{2} / #Sigma w_{i} - <#eta>^{2} vs E",  nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
1065         fhSumEtaE->SetXTitle("E (GeV)");
1066         fhSumEtaE->SetYTitle("#delta^{2}_{#eta #eta}");
1067         outputContainer->Add(fhSumEtaE);
1068
1069         fhSumPhiE  = new TH2F ("hSumPhiE","#sigma^{2}_{#phi #phi} = #Sigma w_{i}(#phi_{i})^{2}/ #Sigma w_{i} - <#phi>^{2} vs E",
1070                                nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
1071         fhSumPhiE->SetXTitle("E (GeV)");
1072         fhSumPhiE->SetYTitle("#delta^{2}_{#phi #phi}");
1073         outputContainer->Add(fhSumPhiE);
1074
1075         fhSumEtaPhiE  = new TH2F ("hSumEtaPhiE","#delta^{2}_{#eta #phi} = #Sigma w_{i}(#phi_{i} #eta_{i} ) / #Sigma w_{i} - <#phi><#eta> vs E",
1076                                   nptbins,ptmin,ptmax, 2*ssbins,-ssmax,ssmax);
1077         fhSumEtaPhiE->SetXTitle("E (GeV)");
1078         fhSumEtaPhiE->SetYTitle("#delta^{2}_{#eta #phi}");
1079         outputContainer->Add(fhSumEtaPhiE);
1080
1081         fhDispEtaPhiDiffE  = new TH2F ("hDispEtaPhiDiffE","#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta} vs E",
1082                                        nptbins,ptmin,ptmax,200, -10,10);
1083         fhDispEtaPhiDiffE->SetXTitle("E (GeV)");
1084         fhDispEtaPhiDiffE->SetYTitle("#sigma^{2}_{#phi #phi}-#sigma^{2}_{#eta #eta}");
1085         outputContainer->Add(fhDispEtaPhiDiffE);
1086
1087         fhSphericityE  = new TH2F ("hSphericityE","(#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi}) vs E",
1088                                    nptbins,ptmin,ptmax, 200, -1,1);
1089         fhSphericityE->SetXTitle("E (GeV)");
1090         fhSphericityE->SetYTitle("s = (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi})");
1091         outputContainer->Add(fhSphericityE);
1092
1093         for(Int_t i = 0; i < 7; i++)
1094         {
1095           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]),
1096                                           ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
1097           fhDispEtaDispPhi[i]->SetXTitle("#sigma^{2}_{#eta #eta}");
1098           fhDispEtaDispPhi[i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1099           outputContainer->Add(fhDispEtaDispPhi[i]);
1100
1101           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]),
1102                                           ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
1103           fhLambda0DispEta[i]->SetXTitle("#lambda^{2}_{0}");
1104           fhLambda0DispEta[i]->SetYTitle("#sigma^{2}_{#eta #eta}");
1105           outputContainer->Add(fhLambda0DispEta[i]);
1106
1107           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]),
1108                                           ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
1109           fhLambda0DispPhi[i]->SetXTitle("#lambda^{2}_{0}");
1110           fhLambda0DispPhi[i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1111           outputContainer->Add(fhLambda0DispPhi[i]);
1112
1113         }
1114       }
1115     }
1116
1117     fhNLocMaxE = new TH2F("hNLocMaxE","Number of local maxima in cluster",
1118                           nptbins,ptmin,ptmax,10,0,10);
1119     fhNLocMaxE ->SetYTitle("N maxima");
1120     fhNLocMaxE ->SetXTitle("E (GeV)");
1121     outputContainer->Add(fhNLocMaxE) ;
1122
1123     if(fAnaType == kSSCalo)
1124     {
1125       fhNLocMaxPt = new TH2F("hNLocMaxPt","Number of local maxima in cluster",
1126                              nptbins,ptmin,ptmax,10,0,10);
1127       fhNLocMaxPt ->SetYTitle("N maxima");
1128       fhNLocMaxPt ->SetXTitle("p_{T} (GeV/c)");
1129       outputContainer->Add(fhNLocMaxPt) ;
1130     }
1131
1132     for (Int_t i = 0; i < 3; i++)
1133     {
1134       fhELambda0LocMax[i]  = new TH2F(Form("hELambda0LocMax%d",i+1),
1135                                       Form("Selected #pi^{0} (#eta) pairs: E vs #lambda_{0}, %s",nlm[i].Data()),
1136                                       nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1137       fhELambda0LocMax[i]->SetYTitle("#lambda_{0}^{2}");
1138       fhELambda0LocMax[i]->SetXTitle("E (GeV)");
1139       outputContainer->Add(fhELambda0LocMax[i]) ;
1140
1141       fhELambda1LocMax[i]  = new TH2F(Form("hELambda1LocMax%d",i+1),
1142                                       Form("Selected #pi^{0} (#eta) pairs: E vs #lambda_{1}, %s",nlm[i].Data()),
1143                                       nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1144       fhELambda1LocMax[i]->SetYTitle("#lambda_{1}^{2}");
1145       fhELambda1LocMax[i]->SetXTitle("E (GeV)");
1146       outputContainer->Add(fhELambda1LocMax[i]) ;
1147
1148       fhEDispersionLocMax[i]  = new TH2F(Form("hEDispersionLocMax%d",i+1),
1149                                          Form("Selected #pi^{0} (#eta) pairs: E vs dispersion^{2}, %s",nlm[i].Data()),
1150                                          nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1151       fhEDispersionLocMax[i]->SetYTitle("dispersion^{2}");
1152       fhEDispersionLocMax[i]->SetXTitle("E (GeV)");
1153       outputContainer->Add(fhEDispersionLocMax[i]) ;
1154
1155       if(fCalorimeter == "EMCAL" && !fFillOnlySimpleSSHisto)
1156       {
1157         fhEDispEtaLocMax[i]  = new TH2F(Form("hEDispEtaLocMax%d",i+1),
1158                                         Form("Selected #pi^{0} (#eta) pairs: E vs #sigma_{#eta #eta}, %s",nlm[i].Data()),
1159                                         nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1160         fhEDispEtaLocMax[i]->SetYTitle("#sigma_{#eta #eta}");
1161         fhEDispEtaLocMax[i]->SetXTitle("E (GeV)");
1162         outputContainer->Add(fhEDispEtaLocMax[i]) ;
1163
1164         fhEDispPhiLocMax[i]  = new TH2F(Form("hEDispPhiLocMax%d",i+1),
1165                                         Form("Selected #pi^{0} (#eta) pairs: E vs #sigma_{#phi #phi}, %s",nlm[i].Data()),
1166                                         nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1167         fhEDispPhiLocMax[i]->SetYTitle("#sigma_{#phi #phi}");
1168         fhEDispPhiLocMax[i]->SetXTitle("E (GeV)");
1169         outputContainer->Add(fhEDispPhiLocMax[i]) ;
1170
1171         fhESumEtaPhiLocMax[i]  = new TH2F(Form("hESumEtaPhiLocMax%d",i+1),
1172                                           Form("Selected #pi^{0} (#eta) pairs: E vs #sigma_{#eta #phi}, %s",nlm[i].Data()),
1173                                           nptbins,ptmin,ptmax,2*ssbins,-ssmax,ssmax);
1174         fhESumEtaPhiLocMax[i]->SetYTitle("#sigma_{#eta #phi}");
1175         fhESumEtaPhiLocMax[i]->SetXTitle("E (GeV)");
1176         outputContainer->Add(fhESumEtaPhiLocMax[i]) ;
1177
1178         fhEDispEtaPhiDiffLocMax[i]  = new TH2F(Form("hEDispEtaPhiDiffLocMax%d",i+1),
1179                                                Form("Selected #pi^{0} (#eta) pairs: E vs #sigma_{#phi #phi} - #sigma_{#eta #eta}, %s",nlm[i].Data()),
1180                                                nptbins,ptmin,ptmax,200, -10,10);
1181         fhEDispEtaPhiDiffLocMax[i]->SetYTitle("#sigma_{#phi #phi} - #sigma_{#eta #eta}");
1182         fhEDispEtaPhiDiffLocMax[i]->SetXTitle("E (GeV)");
1183         outputContainer->Add(fhEDispEtaPhiDiffLocMax[i]) ;
1184
1185         fhESphericityLocMax[i]  = new TH2F(Form("hESphericityLocMax%d",i+1),
1186                                            Form("Selected #pi^{0} (#eta) pairs: E vs #sigma_{#phi #phi} - #sigma_{#eta #eta} / (#sigma_{#phi #phi} + #sigma_{#eta #eta}), %s",nlm[i].Data()),
1187                                            nptbins,ptmin,ptmax,200, -1,1);
1188         fhESphericityLocMax[i]->SetYTitle("#sigma_{#phi #phi} - #sigma_{#eta #eta} / (#sigma_{#phi #phi} + #sigma_{#eta #eta})");
1189         fhESphericityLocMax[i]->SetXTitle("E (GeV)");
1190         outputContainer->Add(fhESphericityLocMax[i]) ;
1191       }
1192
1193     }
1194
1195     fhENCells  = new TH2F ("hENCells","N cells in cluster vs E ", nptbins,ptmin,ptmax, nbins,nmin,nmax);
1196     fhENCells->SetXTitle("E (GeV)");
1197     fhENCells->SetYTitle("# of cells in cluster");
1198     outputContainer->Add(fhENCells);
1199
1200     fhETime = new TH2F("hETime","cluster time vs pair E",nptbins,ptmin,ptmax, tbins,tmin,tmax);
1201     fhETime->SetXTitle("E (GeV)");
1202     fhETime->SetYTitle("t (ns)");
1203     outputContainer->Add(fhETime);
1204
1205   }
1206
1207
1208   fhEPairDiffTime = new TH2F("hEPairDiffTime","cluster pair time difference vs E",nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
1209   fhEPairDiffTime->SetXTitle("E_{pair} (GeV)");
1210   fhEPairDiffTime->SetYTitle("#Delta t (ns)");
1211   outputContainer->Add(fhEPairDiffTime);
1212
1213   if(fAnaType == kIMCalo)
1214   {
1215     TString combiName [] = {"1LocMax","2LocMax","NLocMax","1LocMax2LocMax","1LocMaxNLocMax","2LocMaxNLocMax","1LocMaxSSBad","NLocMaxSSGood"};
1216     TString combiTitle[] = {"1 Local Maxima in both clusters","2 Local Maxima in both clusters","more than 2 Local Maxima in both clusters",
1217       "1 Local Maxima paired with 2 Local Maxima","1 Local Maxima paired with more than 2 Local Maxima",
1218       "2 Local Maxima paired with more than 2 Local Maxima",
1219       "1 Local Maxima paired with #lambda_{0}^{2}>0.3","N Local Maxima paired with 0.1<#lambda_{0}^{2}<0.3"};
1220
1221     for (Int_t i = 0; i < 8 ; i++)
1222     {
1223
1224       if (fAnaType == kIMCaloTracks && i > 2 ) continue ;
1225
1226       fhMassPairLocMax[i]  = new TH2F
1227       (Form("MassPairLocMax%s",combiName[i].Data()),
1228        Form("Mass for decay #gamma pair vs E_{pair}, origin #pi^{0}, %s", combiTitle[i].Data()),
1229        nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
1230       fhMassPairLocMax[i]->SetYTitle("Mass (MeV/c^{2})");
1231       fhMassPairLocMax[i]->SetXTitle("E_{pair} (GeV)");
1232       outputContainer->Add(fhMassPairLocMax[i]) ;
1233     }
1234   }
1235
1236   if(fFillTMHisto)
1237   {
1238     fhTrackMatchedDEta  = new TH2F
1239     ("hTrackMatchedDEta",
1240      "d#eta of cluster-track vs cluster energy",
1241      nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
1242     fhTrackMatchedDEta->SetYTitle("d#eta");
1243     fhTrackMatchedDEta->SetXTitle("E_{cluster} (GeV)");
1244
1245     fhTrackMatchedDPhi  = new TH2F
1246     ("hTrackMatchedDPhi",
1247      "d#phi of cluster-track vs cluster energy",
1248      nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
1249     fhTrackMatchedDPhi->SetYTitle("d#phi (rad)");
1250     fhTrackMatchedDPhi->SetXTitle("E_{cluster} (GeV)");
1251
1252     fhTrackMatchedDEtaDPhi  = new TH2F
1253     ("hTrackMatchedDEtaDPhi",
1254      "d#eta vs d#phi of cluster-track vs cluster energy",
1255      nresetabins,resetamin,resetamax,nresphibins,resphimin,resphimax);
1256     fhTrackMatchedDEtaDPhi->SetYTitle("d#phi (rad)");
1257     fhTrackMatchedDEtaDPhi->SetXTitle("d#eta");
1258
1259     outputContainer->Add(fhTrackMatchedDEta) ;
1260     outputContainer->Add(fhTrackMatchedDPhi) ;
1261     outputContainer->Add(fhTrackMatchedDEtaDPhi) ;
1262
1263     fhTrackMatchedDEtaPos  = new TH2F
1264     ("hTrackMatchedDEtaPos",
1265      "d#eta of cluster-track vs cluster energy",
1266      nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
1267     fhTrackMatchedDEtaPos->SetYTitle("d#eta");
1268     fhTrackMatchedDEtaPos->SetXTitle("E_{cluster} (GeV)");
1269
1270     fhTrackMatchedDPhiPos  = new TH2F
1271     ("hTrackMatchedDPhiPos",
1272      "d#phi of cluster-track vs cluster energy",
1273      nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
1274     fhTrackMatchedDPhiPos->SetYTitle("d#phi (rad)");
1275     fhTrackMatchedDPhiPos->SetXTitle("E_{cluster} (GeV)");
1276
1277     fhTrackMatchedDEtaDPhiPos  = new TH2F
1278     ("hTrackMatchedDEtaDPhiPos",
1279      "d#eta vs d#phi of cluster-track vs cluster energy",
1280      nresetabins,resetamin,resetamax,nresphibins,resphimin,resphimax);
1281     fhTrackMatchedDEtaDPhiPos->SetYTitle("d#phi (rad)");
1282     fhTrackMatchedDEtaDPhiPos->SetXTitle("d#eta");
1283
1284     outputContainer->Add(fhTrackMatchedDEtaPos) ;
1285     outputContainer->Add(fhTrackMatchedDPhiPos) ;
1286     outputContainer->Add(fhTrackMatchedDEtaDPhiPos) ;
1287
1288     fhTrackMatchedDEtaNeg  = new TH2F
1289     ("hTrackMatchedDEtaNeg",
1290      "d#eta of cluster-track vs cluster energy",
1291      nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
1292     fhTrackMatchedDEtaNeg->SetYTitle("d#eta");
1293     fhTrackMatchedDEtaNeg->SetXTitle("E_{cluster} (GeV)");
1294
1295     fhTrackMatchedDPhiNeg  = new TH2F
1296     ("hTrackMatchedDPhiNeg",
1297      "d#phi of cluster-track vs cluster energy",
1298      nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
1299     fhTrackMatchedDPhiNeg->SetYTitle("d#phi (rad)");
1300     fhTrackMatchedDPhiNeg->SetXTitle("E_{cluster} (GeV)");
1301
1302     fhTrackMatchedDEtaDPhiNeg  = new TH2F
1303     ("hTrackMatchedDEtaDPhiNeg",
1304      "d#eta vs d#phi of cluster-track vs cluster energy",
1305      nresetabins,resetamin,resetamax,nresphibins,resphimin,resphimax);
1306     fhTrackMatchedDEtaDPhiNeg->SetYTitle("d#phi (rad)");
1307     fhTrackMatchedDEtaDPhiNeg->SetXTitle("d#eta");
1308
1309     outputContainer->Add(fhTrackMatchedDEtaNeg) ;
1310     outputContainer->Add(fhTrackMatchedDPhiNeg) ;
1311     outputContainer->Add(fhTrackMatchedDEtaDPhiNeg) ;
1312
1313     fhdEdx  = new TH2F ("hdEdx","matched track <dE/dx> vs cluster E ", nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
1314     fhdEdx->SetXTitle("E (GeV)");
1315     fhdEdx->SetYTitle("<dE/dx>");
1316     outputContainer->Add(fhdEdx);
1317
1318     fhEOverP  = new TH2F ("hEOverP","matched track E/p vs cluster E ", nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
1319     fhEOverP->SetXTitle("E (GeV)");
1320     fhEOverP->SetYTitle("E/p");
1321     outputContainer->Add(fhEOverP);
1322
1323     if(fCalorimeter=="EMCAL")
1324     {
1325       fhEOverPNoTRD  = new TH2F ("hEOverPNoTRD","matched track E/p vs cluster E, SM not behind TRD ", nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
1326       fhEOverPNoTRD->SetXTitle("E (GeV)");
1327       fhEOverPNoTRD->SetYTitle("E/p");
1328       outputContainer->Add(fhEOverPNoTRD);
1329     }
1330
1331     if(IsDataMC() && fFillTMHisto)
1332     {
1333       fhTrackMatchedMCParticleE  = new TH2F
1334       ("hTrackMatchedMCParticleE",
1335        "Origin of particle vs energy",
1336        nptbins,ptmin,ptmax,8,0,8);
1337       fhTrackMatchedMCParticleE->SetXTitle("E (GeV)");
1338       //fhTrackMatchedMCParticleE->SetYTitle("Particle type");
1339
1340       fhTrackMatchedMCParticleE->GetYaxis()->SetBinLabel(1 ,"Photon");
1341       fhTrackMatchedMCParticleE->GetYaxis()->SetBinLabel(2 ,"Electron");
1342       fhTrackMatchedMCParticleE->GetYaxis()->SetBinLabel(3 ,"Meson Merged");
1343       fhTrackMatchedMCParticleE->GetYaxis()->SetBinLabel(4 ,"Rest");
1344       fhTrackMatchedMCParticleE->GetYaxis()->SetBinLabel(5 ,"Conv. Photon");
1345       fhTrackMatchedMCParticleE->GetYaxis()->SetBinLabel(6 ,"Conv. Electron");
1346       fhTrackMatchedMCParticleE->GetYaxis()->SetBinLabel(7 ,"Conv. Merged");
1347       fhTrackMatchedMCParticleE->GetYaxis()->SetBinLabel(8 ,"Conv. Rest");
1348
1349       outputContainer->Add(fhTrackMatchedMCParticleE);
1350
1351       fhTrackMatchedMCParticleDEta  = new TH2F
1352       ("hTrackMatchedMCParticleDEta",
1353        "Origin of particle vs #eta residual",
1354        nresetabins,resetamin,resetamax,8,0,8);
1355       fhTrackMatchedMCParticleDEta->SetXTitle("#Delta #eta");
1356       //fhTrackMatchedMCParticleDEta->SetYTitle("Particle type");
1357
1358       fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(1 ,"Photon");
1359       fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(2 ,"Electron");
1360       fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(3 ,"Meson Merged");
1361       fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(4 ,"Rest");
1362       fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(5 ,"Conv. Photon");
1363       fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(6 ,"Conv. Electron");
1364       fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(7 ,"Conv. Merged");
1365       fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(8 ,"Conv. Rest");
1366
1367       outputContainer->Add(fhTrackMatchedMCParticleDEta);
1368
1369       fhTrackMatchedMCParticleDPhi  = new TH2F
1370       ("hTrackMatchedMCParticleDPhi",
1371        "Origin of particle vs #phi residual",
1372        nresphibins,resphimin,resphimax,8,0,8);
1373       fhTrackMatchedMCParticleDPhi->SetXTitle("#Delta #phi");
1374       //fhTrackMatchedMCParticleDPhi->SetYTitle("Particle type");
1375
1376       fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(1 ,"Photon");
1377       fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(2 ,"Electron");
1378       fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(3 ,"Meson Merged");
1379       fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(4 ,"Rest");
1380       fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(5 ,"Conv. Photon");
1381       fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(6 ,"Conv. Electron");
1382       fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(7 ,"Conv. Merged");
1383       fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(8 ,"Conv. Rest");
1384
1385       outputContainer->Add(fhTrackMatchedMCParticleDPhi);
1386
1387
1388     }
1389   }
1390
1391   if(fFillWeightHistograms)
1392   {
1393     fhECellClusterRatio  = new TH2F ("hECellClusterRatio"," cell energy / cluster energy vs cluster energy, for selected decay photons from neutral meson",
1394                                      nptbins,ptmin,ptmax, 100,0,1.);
1395     fhECellClusterRatio->SetXTitle("E_{cluster} (GeV) ");
1396     fhECellClusterRatio->SetYTitle("E_{cell i}/E_{cluster}");
1397     outputContainer->Add(fhECellClusterRatio);
1398
1399     fhECellClusterLogRatio  = new TH2F ("hECellClusterLogRatio"," Log(cell energy / cluster energy) vs cluster energy, for selected decay photons from neutral meson",
1400                                         nptbins,ptmin,ptmax, 100,-10,0);
1401     fhECellClusterLogRatio->SetXTitle("E_{cluster} (GeV) ");
1402     fhECellClusterLogRatio->SetYTitle("Log (E_{max cell}/E_{cluster})");
1403     outputContainer->Add(fhECellClusterLogRatio);
1404
1405     fhEMaxCellClusterRatio  = new TH2F ("hEMaxCellClusterRatio"," max cell energy / cluster energy vs cluster energy, for selected decay photons from neutral meson",
1406                                         nptbins,ptmin,ptmax, 100,0,1.);
1407     fhEMaxCellClusterRatio->SetXTitle("E_{cluster} (GeV) ");
1408     fhEMaxCellClusterRatio->SetYTitle("E_{max cell}/E_{cluster}");
1409     outputContainer->Add(fhEMaxCellClusterRatio);
1410
1411     fhEMaxCellClusterLogRatio  = new TH2F ("hEMaxCellClusterLogRatio"," Log(max cell energy / cluster energy) vs cluster energy, for selected decay photons from neutral meson",
1412                                            nptbins,ptmin,ptmax, 100,-10,0);
1413     fhEMaxCellClusterLogRatio->SetXTitle("E_{cluster} (GeV) ");
1414     fhEMaxCellClusterLogRatio->SetYTitle("Log (E_{max cell}/E_{cluster})");
1415     outputContainer->Add(fhEMaxCellClusterLogRatio);
1416
1417     for(Int_t iw = 0; iw < 14; iw++)
1418     {
1419       fhLambda0ForW0[iw]  = new TH2F (Form("hLambda0ForW0%d",iw),Form("shower shape, #lambda^{2}_{0} vs E, w0 = %1.1f, for selected decay photons from neutral meson",1+0.5*iw),
1420                                       nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1421       fhLambda0ForW0[iw]->SetXTitle("E_{cluster}");
1422       fhLambda0ForW0[iw]->SetYTitle("#lambda^{2}_{0}");
1423       outputContainer->Add(fhLambda0ForW0[iw]);
1424
1425       //      fhLambda1ForW0[iw]  = new TH2F (Form("hLambda1ForW0%d",iw),Form("shower shape, #lambda^{2}_{1} vs E, w0 = %1.1f, for selected decay photons from neutral meson",0.5+0.5*iw),
1426       //                                      nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1427       //      fhLambda1ForW0[iw]->SetXTitle("E_{cluster}");
1428       //      fhLambda1ForW0[iw]->SetYTitle("#lambda^{2}_{1}");
1429       //      outputContainer->Add(fhLambda1ForW0[iw]);
1430
1431     }
1432   }
1433
1434   if(IsDataMC())
1435   {
1436     if(GetReader()->GetDataType() != AliCaloTrackReader::kMC && fAnaType==kSSCalo)
1437     {
1438       fhMCPi0PtGenRecoFraction = new TH2F("hMCPi0PtGenRecoFraction","Number of clusters from #pi^{0} (2 #gamma) identified as #pi^{0} (#eta), pT versus E primary #pi^{0} / E reco",
1439                                           nptbins,ptmin,ptmax,200,0,2);
1440       fhMCPi0PtGenRecoFraction->SetXTitle("p^{rec}_{T} (GeV/c)");
1441       fhMCPi0PtGenRecoFraction->SetYTitle("E^{ #pi^{0} mother} / E^{rec}");
1442       outputContainer->Add(fhMCPi0PtGenRecoFraction) ;
1443
1444       fhMCEtaPtGenRecoFraction = new TH2F("hMCEtaPtGenRecoFraction","Number of clusters from #eta (2 #gamma) identified as #pi^{0} (#eta),pT versus E primary #eta / E reco",
1445                                           nptbins,ptmin,ptmax,200,0,2);
1446       fhMCEtaPtGenRecoFraction->SetXTitle("p^{rec}_{T} (GeV/c)");
1447       fhMCEtaPtGenRecoFraction->SetYTitle("E^{ #eta mother} / E^{rec}");
1448       outputContainer->Add(fhMCEtaPtGenRecoFraction) ;
1449
1450       fhMCPi0DecayPt = new TH1F("hMCPi0DecayPt","Number of #gamma from #pi^{0} decay  identified as #pi^{0} (#eta)",nptbins,ptmin,ptmax);
1451       fhMCPi0DecayPt->SetYTitle("N");
1452       fhMCPi0DecayPt->SetXTitle("p^{rec}_{T} (GeV/c)");
1453       outputContainer->Add(fhMCPi0DecayPt) ;
1454
1455       fhMCPi0DecayPtFraction = new TH2F("hMCPi0DecayPtFraction","Number of #gamma from #pi^{0} decay  identified as #pi^{0} (#eta), pT versus E primary  #gamma / E primary #pi^{0}",
1456                                         nptbins,ptmin,ptmax,100,0,1);
1457       fhMCPi0DecayPtFraction->SetXTitle("p^{rec}_{T} (GeV/c)");
1458       fhMCPi0DecayPtFraction->SetYTitle("E^{gen} / E^{gen-mother}");
1459       outputContainer->Add(fhMCPi0DecayPtFraction) ;
1460
1461       fhMCEtaDecayPt = new TH1F("hMCEtaDecayPt","Number of #gamma from #eta decay  identified as #pi^{0} (#eta)",nptbins,ptmin,ptmax);
1462       fhMCEtaDecayPt->SetYTitle("N");
1463       fhMCEtaDecayPt->SetXTitle("p^{rec}_{T} (GeV/c)");
1464       outputContainer->Add(fhMCEtaDecayPt) ;
1465
1466       fhMCEtaDecayPtFraction = new TH2F("hMCEtaDecayPtFraction","Number of #gamma from #eta decay  identified as #pi^{0} (#eta), pT versus E primary  #gamma / E primary #eta",
1467                                         nptbins,ptmin,ptmax,100,0,1);
1468       fhMCEtaDecayPtFraction->SetXTitle("p^{rec}_{T} (GeV/c)");
1469       fhMCEtaDecayPtFraction->SetYTitle("E^{gen} / E^{gen-mother}");
1470       outputContainer->Add(fhMCEtaDecayPtFraction) ;
1471
1472       fhMCOtherDecayPt = new TH1F("hMCOtherDecayPt","Number of #gamma decay (not #eta or #pi^{0})  identified as #pi^{0} (#eta)",nptbins,ptmin,ptmax);
1473       fhMCOtherDecayPt->SetYTitle("N");
1474       fhMCOtherDecayPt->SetXTitle("p^{rec}_{T} (GeV/c)");
1475       outputContainer->Add(fhMCOtherDecayPt) ;
1476
1477     }
1478
1479     if((GetReader()->GetDataType() == AliCaloTrackReader::kMC && fAnaType!=kSSCalo) ||
1480        GetReader()->GetDataType() != AliCaloTrackReader::kMC)
1481     {
1482
1483       fhAnglePairMCPi0  = new TH2F
1484       ("AnglePairMCPi0",
1485        "Angle between decay #gamma pair vs E_{pair}, origin #pi^{0}",nptbins,ptmin,ptmax,250,0,0.5);
1486       fhAnglePairMCPi0->SetYTitle("#alpha (rad)");
1487       fhAnglePairMCPi0->SetXTitle("E_{pair} (GeV)");
1488       outputContainer->Add(fhAnglePairMCPi0) ;
1489
1490       if (fAnaType!= kSSCalo)
1491       {
1492         fhAnglePairMCEta  = new TH2F
1493         ("AnglePairMCEta",
1494          "Angle between decay #gamma pair vs E_{pair}, origin #eta",nptbins,ptmin,ptmax,250,0,0.5);
1495         fhAnglePairMCEta->SetYTitle("#alpha (rad)");
1496         fhAnglePairMCEta->SetXTitle("E_{pair} (GeV)");
1497         outputContainer->Add(fhAnglePairMCEta) ;
1498
1499         fhMassPairMCPi0  = new TH2F
1500         ("MassPairMCPi0",
1501          "Mass for decay #gamma pair vs E_{pair}, origin #pi^{0}",nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
1502         fhMassPairMCPi0->SetYTitle("Mass (MeV/c^{2})");
1503         fhMassPairMCPi0->SetXTitle("E_{pair} (GeV)");
1504         outputContainer->Add(fhMassPairMCPi0) ;
1505
1506         fhMassPairMCEta  = new TH2F
1507         ("MassPairMCEta",
1508          "Mass for decay #gamma pair vs E_{pair}, origin #eta",nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
1509         fhMassPairMCEta->SetYTitle("Mass (MeV/c^{2})");
1510         fhMassPairMCEta->SetXTitle("E_{pair} (GeV)");
1511         outputContainer->Add(fhMassPairMCEta) ;
1512       }
1513
1514       for(Int_t i = 0; i < 6; i++)
1515       {
1516
1517         fhMCE[i]  = new TH1F
1518         (Form("hE_MC%s",pname[i].Data()),
1519          Form("Identified as #pi^{0} (#eta), cluster from %s",
1520               ptype[i].Data()),
1521          nptbins,ptmin,ptmax);
1522         fhMCE[i]->SetYTitle("N");
1523         fhMCE[i]->SetXTitle("E (GeV)");
1524         outputContainer->Add(fhMCE[i]) ;
1525
1526         fhMCPt[i]  = new TH1F
1527         (Form("hPt_MC%s",pname[i].Data()),
1528          Form("Identified as #pi^{0} (#eta), cluster from %s",
1529               ptype[i].Data()),
1530          nptbins,ptmin,ptmax);
1531         fhMCPt[i]->SetYTitle("N");
1532         fhMCPt[i]->SetXTitle("p_{T} (GeV/c)");
1533         outputContainer->Add(fhMCPt[i]) ;
1534
1535         fhMCPtCentrality[i]  = new TH2F
1536         (Form("hPtCentrality_MC%s",pname[i].Data()),
1537          Form("Identified as #pi^{0} (#eta), cluster from %s",
1538               ptype[i].Data()),
1539          nptbins,ptmin,ptmax, 100,0,100);
1540         fhMCPtCentrality[i]->SetYTitle("centrality");
1541         fhMCPtCentrality[i]->SetXTitle("p_{T} (GeV/c)");
1542         outputContainer->Add(fhMCPtCentrality[i]) ;
1543
1544         if(fAnaType == kSSCalo)
1545         {
1546
1547           fhMCNLocMaxPt[i] = new TH2F
1548           (Form("hNLocMaxPt_MC%s",pname[i].Data()),
1549            Form("cluster from %s, pT of cluster, for NLM",ptype[i].Data()),
1550            nptbins,ptmin,ptmax,10,0,10);
1551           fhMCNLocMaxPt[i] ->SetYTitle("N maxima");
1552           fhMCNLocMaxPt[i] ->SetXTitle("p_{T} (GeV/c)");
1553           outputContainer->Add(fhMCNLocMaxPt[i]) ;
1554
1555           fhMCEReject[i]  = new TH1F
1556           (Form("hEReject_MC%s",pname[i].Data()),
1557            Form("Rejected as #pi^{0} (#eta), cluster from %s",
1558                 ptype[i].Data()),
1559            nptbins,ptmin,ptmax);
1560           fhMCEReject[i]->SetYTitle("N");
1561           fhMCEReject[i]->SetXTitle("E (GeV)");
1562           outputContainer->Add(fhMCEReject[i]) ;
1563
1564           fhMCPtReject[i]  = new TH1F
1565           (Form("hPtReject_MC%s",pname[i].Data()),
1566            Form("Rejected as #pi^{0} (#eta), cluster from %s",
1567                 ptype[i].Data()),
1568            nptbins,ptmin,ptmax);
1569           fhMCPtReject[i]->SetYTitle("N");
1570           fhMCPtReject[i]->SetXTitle("p_{T} (GeV/c)");
1571           outputContainer->Add(fhMCPtReject[i]) ;
1572         }
1573
1574         fhMCPhi[i]  = new TH2F
1575         (Form("hPhi_MC%s",pname[i].Data()),
1576          Form("Identified as #pi^{0} (#eta), cluster from %s",ptype[i].Data()),
1577          nptbins,ptmin,ptmax,nphibins,phimin,phimax);
1578         fhMCPhi[i]->SetYTitle("#phi");
1579         fhMCPhi[i]->SetXTitle("p_{T} (GeV/c)");
1580         outputContainer->Add(fhMCPhi[i]) ;
1581
1582         fhMCEta[i]  = new TH2F
1583         (Form("hEta_MC%s",pname[i].Data()),
1584          Form("Identified as #pi^{0} (#eta), cluster from %s",
1585               ptype[i].Data()),nptbins,ptmin,ptmax,netabins,etamin,etamax);
1586         fhMCEta[i]->SetYTitle("#eta");
1587         fhMCEta[i]->SetXTitle("p_{T} (GeV/c)");
1588         outputContainer->Add(fhMCEta[i]) ;
1589
1590         fhMCMassPt[i]  = new TH2F
1591         (Form("hMassPt_MC%s",pname[i].Data()),
1592          Form("all pairs mass: p_{T} vs massfrom %s",ptype[i].Data()),
1593          nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1594         fhMCMassPt[i]->SetYTitle("mass (GeV/c^{2})");
1595         fhMCMassPt[i]->SetXTitle("p_{T} (GeV/c)");
1596         outputContainer->Add(fhMCMassPt[i]) ;
1597
1598         fhMCSelectedMassPt[i]  = new TH2F
1599         (Form("hSelectedMassPt_MC%s",pname[i].Data()),
1600          Form("Selected #pi^{0} (#eta) pairs mass: p_{T} vs massfrom %s",ptype[i].Data()),
1601          nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1602         fhMCSelectedMassPt[i]->SetYTitle("mass (GeV/c^{2})");
1603         fhMCSelectedMassPt[i]->SetXTitle("p_{T} (GeV/c)");
1604         outputContainer->Add(fhMCSelectedMassPt[i]) ;
1605
1606
1607         if( fFillSelectClHisto )
1608         {
1609           fhEMCLambda0[i]  = new TH2F(Form("hELambda0_MC%s",pname[i].Data()),
1610                                       Form("Selected pair, cluster from %s : E vs #lambda_{0}^{2}",ptype[i].Data()),
1611                                       nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1612           fhEMCLambda0[i]->SetYTitle("#lambda_{0}^{2}");
1613           fhEMCLambda0[i]->SetXTitle("E (GeV)");
1614           outputContainer->Add(fhEMCLambda0[i]) ;
1615
1616           fhEMCLambda1[i]  = new TH2F(Form("hELambda1_MC%s",pname[i].Data()),
1617                                       Form("Selected pair, cluster from %s : E vs #lambda_{1}^{2}",ptype[i].Data()),
1618                                       nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1619           fhEMCLambda1[i]->SetYTitle("#lambda_{1}^{2}");
1620           fhEMCLambda1[i]->SetXTitle("E (GeV)");
1621           outputContainer->Add(fhEMCLambda1[i]) ;
1622
1623           fhEMCDispersion[i]  = new TH2F(Form("hEDispersion_MC%s",pname[i].Data()),
1624                                          Form("Selected pair, cluster from %s : E vs dispersion^{2}",ptype[i].Data()),
1625                                          nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1626           fhEMCDispersion[i]->SetYTitle("D^{2}");
1627           fhEMCDispersion[i]->SetXTitle("E (GeV)");
1628           outputContainer->Add(fhEMCDispersion[i]) ;
1629
1630           if(fCalorimeter=="EMCAL")
1631           {
1632             fhEMCLambda0NoTRD[i]  = new TH2F(Form("hELambda0NoTRD_MC%s",pname[i].Data()),
1633                                              Form("Selected pair, cluster from %s : E vs #lambda_{0}^{2}, NoTRD",ptype[i].Data()),
1634                                              nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1635             fhEMCLambda0NoTRD[i]->SetYTitle("#lambda_{0}^{2}");
1636             fhEMCLambda0NoTRD[i]->SetXTitle("E (GeV)");
1637             outputContainer->Add(fhEMCLambda0NoTRD[i]) ;
1638
1639             if(!fFillOnlySimpleSSHisto)
1640             {
1641               fhMCEDispEta[i]  = new TH2F (Form("hEDispEtaE_MC%s",pname[i].Data()),
1642                                            Form("cluster from %s : #sigma^{2}_{#eta #eta} = #Sigma w_{i}(#eta_{i} - <#eta>)^{2}/ #Sigma w_{i} vs E",ptype[i].Data()),
1643                                            nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
1644               fhMCEDispEta[i]->SetXTitle("E (GeV)");
1645               fhMCEDispEta[i]->SetYTitle("#sigma^{2}_{#eta #eta}");
1646               outputContainer->Add(fhMCEDispEta[i]);
1647
1648               fhMCEDispPhi[i]  = new TH2F (Form("hEDispPhiE_MC%s",pname[i].Data()),
1649                                            Form("cluster from %s : #sigma^{2}_{#phi #phi} = #Sigma w_{i}(#phi_{i} - <#phi>)^{2} / #Sigma w_{i} vs E",ptype[i].Data()),
1650                                            nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
1651               fhMCEDispPhi[i]->SetXTitle("E (GeV)");
1652               fhMCEDispPhi[i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1653               outputContainer->Add(fhMCEDispPhi[i]);
1654
1655               fhMCESumEtaPhi[i]  = new TH2F (Form("hESumEtaPhiE_MC%s",pname[i].Data()),
1656                                              Form("cluster from %s : #delta^{2}_{#eta #phi} = #Sigma w_{i}(#phi_{i} #eta_{i} ) / #Sigma w_{i} - <#phi><#eta> vs E",ptype[i].Data()),
1657                                              nptbins,ptmin,ptmax, 2*ssbins,-ssmax,ssmax);
1658               fhMCESumEtaPhi[i]->SetXTitle("E (GeV)");
1659               fhMCESumEtaPhi[i]->SetYTitle("#delta^{2}_{#eta #phi}");
1660               outputContainer->Add(fhMCESumEtaPhi[i]);
1661
1662               fhMCEDispEtaPhiDiff[i]  = new TH2F (Form("hEDispEtaPhiDiffE_MC%s",pname[i].Data()),
1663                                                   Form("cluster from %s : #sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta} vs E",ptype[i].Data()),
1664                                                   nptbins,ptmin,ptmax,200,-10,10);
1665               fhMCEDispEtaPhiDiff[i]->SetXTitle("E (GeV)");
1666               fhMCEDispEtaPhiDiff[i]->SetYTitle("#sigma^{2}_{#phi #phi}-#sigma^{2}_{#eta #eta}");
1667               outputContainer->Add(fhMCEDispEtaPhiDiff[i]);
1668
1669               fhMCESphericity[i]  = new TH2F (Form("hESphericity_MC%s",pname[i].Data()),
1670                                               Form("cluster from %s : (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi}) vs E",ptype[i].Data()),
1671                                               nptbins,ptmin,ptmax, 200,-1,1);
1672               fhMCESphericity[i]->SetXTitle("E (GeV)");
1673               fhMCESphericity[i]->SetYTitle("s = (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi})");
1674               outputContainer->Add(fhMCESphericity[i]);
1675
1676               for(Int_t ie = 0; ie < 7; ie++)
1677               {
1678                 fhMCDispEtaDispPhi[ie][i] = new TH2F (Form("hMCDispEtaDispPhi_EBin%d_MC%s",ie,pname[i].Data()),
1679                                                       Form("cluster from %s : #sigma^{2}_{#phi #phi} vs #sigma^{2}_{#eta #eta} for %d < E < %d GeV",pname[i].Data(),bin[ie],bin[ie+1]),
1680                                                       ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
1681                 fhMCDispEtaDispPhi[ie][i]->SetXTitle("#sigma^{2}_{#eta #eta}");
1682                 fhMCDispEtaDispPhi[ie][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1683                 outputContainer->Add(fhMCDispEtaDispPhi[ie][i]);
1684
1685                 fhMCLambda0DispEta[ie][i] = new TH2F (Form("hMCLambda0DispEta_EBin%d_MC%s",ie,pname[i].Data()),
1686                                                       Form("cluster from %s : #lambda^{2}_{0} vs #sigma^{2}_{#eta #eta} for %d < E < %d GeV",pname[i].Data(),bin[ie],bin[ie+1]),
1687                                                       ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
1688                 fhMCLambda0DispEta[ie][i]->SetXTitle("#lambda^{2}_{0}");
1689                 fhMCLambda0DispEta[ie][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1690                 outputContainer->Add(fhMCLambda0DispEta[ie][i]);
1691
1692                 fhMCLambda0DispPhi[ie][i] = new TH2F (Form("hMCLambda0DispPhi_EBin%d_MC%s",ie,pname[i].Data()),
1693                                                       Form("cluster from %s :#lambda^{2}_{0} vs #sigma^{2}_{#phi #phi} for %d < E < %d GeV",pname[i].Data(),bin[ie],bin[ie+1]),
1694                                                       ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
1695                 fhMCLambda0DispPhi[ie][i]->SetXTitle("#lambda^{2}_{0}");
1696                 fhMCLambda0DispPhi[ie][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1697                 outputContainer->Add(fhMCLambda0DispPhi[ie][i]);
1698
1699               }
1700             }
1701           }
1702
1703           fhEMCLambda0FracMaxCellCut[i]  = new TH2F(Form("hELambda0FracMaxCellCut_MC%s",pname[i].Data()),
1704                                                     Form("Selected pair, cluster from %s : E vs #lambda_{0}^{2}, Max cell fraction of energy < 0.5 ",ptype[i].Data()),
1705                                                     nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1706           fhEMCLambda0FracMaxCellCut[i]->SetYTitle("#lambda_{0}^{2}");
1707           fhEMCLambda0FracMaxCellCut[i]->SetXTitle("E (GeV)");
1708           outputContainer->Add(fhEMCLambda0FracMaxCellCut[i]) ;
1709
1710           fhEMCFracMaxCell[i]  = new TH2F(Form("hEFracMaxCell_MC%s",pname[i].Data()),
1711                                           Form("Selected pair, cluster from %s : E vs Max cell fraction of energy",ptype[i].Data()),
1712                                           nptbins,ptmin,ptmax,100,0,1);
1713           fhEMCFracMaxCell[i]->SetYTitle("Fraction");
1714           fhEMCFracMaxCell[i]->SetXTitle("E (GeV)");
1715           outputContainer->Add(fhEMCFracMaxCell[i]) ;
1716
1717         }//
1718       } // shower shape histo
1719
1720     } //Not MC reader
1721   }//Histos with MC
1722
1723   if(fAnaType==kSSCalo)
1724   {
1725     fhAsymmetry  = new TH2F ("hAsymmetry","A = ( E1 - E2 ) / ( E1 + E2 ) vs E",
1726                              nptbins,ptmin,ptmax, 200, -1,1);
1727     fhAsymmetry->SetXTitle("E (GeV)");
1728     fhAsymmetry->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1729     outputContainer->Add(fhAsymmetry);
1730
1731     fhSelectedAsymmetry  = new TH2F ("hSelectedAsymmetry","A = ( E1 - E2 ) / ( E1 + E2 ) vs E",
1732                                      nptbins,ptmin,ptmax, 200, -1,1);
1733     fhSelectedAsymmetry->SetXTitle("E (GeV)");
1734     fhSelectedAsymmetry->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1735     outputContainer->Add(fhSelectedAsymmetry);
1736
1737     fhSplitE  = new TH1F
1738     ("hSplitE","Selected #pi^{0} (#eta) pairs energy sum of split sub-clusters",nptbins,ptmin,ptmax);
1739     fhSplitE->SetYTitle("counts");
1740     fhSplitE->SetXTitle("E (GeV)");
1741     outputContainer->Add(fhSplitE) ;
1742
1743     fhSplitPt  = new TH1F
1744     ("hSplitPt","Selected #pi^{0} (#eta) pairs pT sum of split sub-clusters",nptbins,ptmin,ptmax);
1745     fhSplitPt->SetYTitle("counts");
1746     fhSplitPt->SetXTitle("p_{T} (GeV/c)");
1747     outputContainer->Add(fhSplitPt) ;
1748
1749
1750     fhSplitPtPhi  = new TH2F
1751     ("hSplitPtPhi","Selected #pi^{0} (#eta) pairs: sum split sub-cluster p_{T} vs #phi",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
1752     fhSplitPtPhi->SetYTitle("#phi (rad)");
1753     fhSplitPtPhi->SetXTitle("p_{T} (GeV/c)");
1754     outputContainer->Add(fhSplitPtPhi) ;
1755
1756     fhSplitPtEta  = new TH2F
1757     ("hSplitPtEta","Selected #pi^{0} (#eta) pairs: sum split sub-cluster p_{T} vs #eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
1758     fhSplitPtEta->SetYTitle("#eta");
1759     fhSplitPtEta->SetXTitle("p_{T} (GeV/c)");
1760     outputContainer->Add(fhSplitPtEta) ;
1761
1762
1763     fhNLocMaxSplitPt = new TH2F("hNLocMaxSplitPt","Number of local maxima in cluster",
1764                                 nptbins,ptmin,ptmax,10,0,10);
1765     fhNLocMaxSplitPt ->SetYTitle("N maxima");
1766     fhNLocMaxSplitPt ->SetXTitle("p_{T} (GeV/c)");
1767     outputContainer->Add(fhNLocMaxSplitPt) ;
1768
1769
1770     fhMassSplitPt  = new TH2F
1771     ("hMassSplitPt","all pairs mass: sum split sub-cluster p_{T} vs mass",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1772     fhMassSplitPt->SetYTitle("mass (GeV/c^{2})");
1773     fhMassSplitPt->SetXTitle("p_{T} (GeV/c)");
1774     outputContainer->Add(fhMassSplitPt) ;
1775
1776     fhSelectedMassSplitPt  = new TH2F
1777     ("hSelectedMassSplitPt","Selected #pi^{0} (#eta) pairs mass: sum split sub-cluster p_{T} vs mass",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1778     fhSelectedMassSplitPt->SetYTitle("mass (GeV/c^{2})");
1779     fhSelectedMassSplitPt->SetXTitle("p_{T} (GeV/c)");
1780     outputContainer->Add(fhSelectedMassSplitPt) ;
1781
1782
1783
1784     if(IsDataMC())
1785     {
1786       for(Int_t i = 0; i< 6; i++)
1787       {
1788         fhMCEAsymmetry[i]  = new TH2F (Form("hEAsymmetry_MC%s",pname[i].Data()),
1789                                        Form("cluster from %s : A = ( E1 - E2 ) / ( E1 + E2 ) vs E",ptype[i].Data()),
1790                                        nptbins,ptmin,ptmax, 200,-1,1);
1791         fhMCEAsymmetry[i]->SetXTitle("E (GeV)");
1792         fhMCEAsymmetry[i]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1793         outputContainer->Add(fhMCEAsymmetry[i]);
1794
1795         fhMCSplitE[i]  = new TH1F
1796         (Form("hSplitE_MC%s",pname[i].Data()),
1797          Form("cluster from %s, energy sum of split sub-clusters",ptype[i].Data()),
1798          nptbins,ptmin,ptmax);
1799         fhMCSplitE[i]->SetYTitle("counts");
1800         fhMCSplitE[i]->SetXTitle("E (GeV)");
1801         outputContainer->Add(fhMCSplitE[i]) ;
1802
1803         fhMCSplitPt[i]  = new TH1F
1804         (Form("hSplitPt_MC%s",pname[i].Data()),
1805          Form("cluster from %s, pT sum of split sub-clusters",ptype[i].Data()),
1806          nptbins,ptmin,ptmax);
1807         fhMCSplitPt[i]->SetYTitle("counts");
1808         fhMCSplitPt[i]->SetXTitle("p_{T} (GeV/c)");
1809         outputContainer->Add(fhMCSplitPt[i]) ;
1810
1811
1812         fhMCSplitPtPhi[i]  = new TH2F
1813         (Form("hSplitPtPhi_MC%s",pname[i].Data()),
1814          Form("Identified as #pi^{0} (#eta), cluster from %s",ptype[i].Data()),
1815          nptbins,ptmin,ptmax,nphibins,phimin,phimax);
1816         fhMCSplitPtPhi[i]->SetYTitle("#phi");
1817         fhMCSplitPtPhi[i]->SetXTitle("p_{T} (GeV/c)");
1818         outputContainer->Add(fhMCSplitPtPhi[i]) ;
1819
1820         fhMCSplitPtEta[i]  = new TH2F
1821         (Form("hSplitPtEta_MC%s",pname[i].Data()),
1822          Form("Identified as #pi^{0} (#eta), cluster from %s",
1823               ptype[i].Data()),nptbins,ptmin,ptmax,netabins,etamin,etamax);
1824         fhMCSplitPtEta[i]->SetYTitle("#eta");
1825         fhMCSplitPtEta[i]->SetXTitle("p_{T} (GeV/c)");
1826         outputContainer->Add(fhMCSplitPtEta[i]) ;
1827
1828
1829         fhMCNLocMaxSplitPt[i] = new TH2F
1830         (Form("hNLocMaxSplitPt_MC%s",pname[i].Data()),
1831          Form("cluster from %s, pT sum of split sub-clusters, for NLM",ptype[i].Data()),
1832          nptbins,ptmin,ptmax,10,0,10);
1833         fhMCNLocMaxSplitPt[i] ->SetYTitle("N maxima");
1834         fhMCNLocMaxSplitPt[i] ->SetXTitle("p_{T} (GeV/c)");
1835         outputContainer->Add(fhMCNLocMaxSplitPt[i]) ;
1836
1837         fhMCMassSplitPt[i]  = new TH2F
1838         (Form("hMassSplitPt_MC%s",pname[i].Data()),
1839          Form("all pairs mass: split p_{T} vs mass from %s",ptype[i].Data()),
1840          nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1841         fhMCMassSplitPt[i]->SetYTitle("mass (GeV/c^{2})");
1842         fhMCMassSplitPt[i]->SetXTitle("p_{T} (GeV/c)");
1843         outputContainer->Add(fhMCMassSplitPt[i]) ;
1844
1845         fhMCSelectedMassSplitPt[i]  = new TH2F
1846         (Form("hSelectedMassSplitPt_MC%s",pname[i].Data()),
1847          Form("Selected #pi^{0} (#eta) pairs mass: split p_{T} vs mass from %s",ptype[i].Data()),
1848          nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1849         fhMCSelectedMassSplitPt[i]->SetYTitle("mass (GeV/c^{2})");
1850         fhMCSelectedMassSplitPt[i]->SetXTitle("p_{T} (GeV/c)");
1851         outputContainer->Add(fhMCSelectedMassSplitPt[i]) ;
1852
1853       }
1854     }
1855   }
1856
1857   if(fAnaType==kSSCalo && fFillSelectClHisto && !fFillOnlySimpleSSHisto )
1858   {
1859
1860
1861     for(Int_t i = 0; i< 3; i++)
1862     {
1863       fhEAsymmetryLocMax[i]  = new TH2F(Form("hEAsymmetryLocMax%d",i+1),
1864                                         Form("Selected #pi^{0} (#eta) pairs: E vs A = ( E1 - E2 ) / ( E1 + E2 ), %s",nlm[i].Data()),
1865                                         nptbins,ptmin,ptmax,200, -1,1);
1866       fhEAsymmetryLocMax[i]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1867       fhEAsymmetryLocMax[i]->SetXTitle("E (GeV)");
1868       outputContainer->Add(fhEAsymmetryLocMax[i]) ;
1869     }
1870
1871     for(Int_t ie = 0; ie< 7; ie++)
1872     {
1873
1874       fhAsymmetryLambda0[ie] = new TH2F (Form("hAsymmetryLambda0_EBin%d",ie),
1875                                          Form("#lambda_{0}^{2} vs A for %d < E < %d GeV",bin[ie],bin[ie+1]),
1876                                          ssbins,ssmin,ssmax , 200,-1,1);
1877       fhAsymmetryLambda0[ie]->SetXTitle("#lambda_{0}^{2}");
1878       fhAsymmetryLambda0[ie]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1879       outputContainer->Add(fhAsymmetryLambda0[ie]);
1880
1881       fhAsymmetryDispEta[ie] = new TH2F (Form("hAsymmetryDispEta_EBin%d",ie),
1882                                          Form("#sigma^{2}_{#eta #eta} vs A for %d < E < %d GeV",bin[ie],bin[ie+1]),
1883                                          ssbins,ssmin,ssmax , 200,-1,1);
1884       fhAsymmetryDispEta[ie]->SetXTitle("#sigma^{2}_{#eta #eta}");
1885       fhAsymmetryDispEta[ie]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1886       outputContainer->Add(fhAsymmetryDispEta[ie]);
1887
1888       fhAsymmetryDispPhi[ie] = new TH2F (Form("hAsymmetryDispPhi_EBin%d",ie),
1889                                          Form("#sigma^{2}_{#phi #phi} vs A for %d < E < %d GeV",bin[ie],bin[ie+1]),
1890                                          ssbins,ssmin,ssmax , 200,-1,1);
1891       fhAsymmetryDispPhi[ie]->SetXTitle("#sigma^{2}_{#phi #phi}");
1892       fhAsymmetryDispPhi[ie]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1893       outputContainer->Add(fhAsymmetryDispPhi[ie]);
1894     }
1895
1896
1897     if(IsDataMC())
1898     {
1899       for(Int_t i = 0; i< 6; i++)
1900       {
1901         for(Int_t ie = 0; ie < 7; ie++)
1902         {
1903           fhMCAsymmetryLambda0[ie][i] = new TH2F (Form("hMCAsymmetryLambda0_EBin%d_MC%s",ie,pname[i].Data()),
1904                                                   Form("cluster from %s : #lambda_{0}^{2} vs A for %d < E < %d GeV",pname[i].Data(),bin[ie],bin[ie+1]),
1905                                                   ssbins,ssmin,ssmax , 200,-1,1);
1906           fhMCAsymmetryLambda0[ie][i]->SetXTitle("#lambda_{0}^{2}");
1907           fhMCAsymmetryLambda0[ie][i]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1908           outputContainer->Add(fhMCAsymmetryLambda0[ie][i]);
1909
1910           fhMCAsymmetryDispEta[ie][i] = new TH2F (Form("hMCAsymmetryDispEta_EBin%d_MC%s",ie,pname[i].Data()),
1911                                                   Form("cluster from %s : #sigma^{2}_{#eta #eta} vs A for %d < E < %d GeV",pname[i].Data(),bin[ie],bin[ie+1]),
1912                                                   ssbins,ssmin,ssmax , 200,-1,1);
1913           fhMCAsymmetryDispEta[ie][i]->SetXTitle("#sigma^{2}_{#eta #eta}");
1914           fhMCAsymmetryDispEta[ie][i]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1915           outputContainer->Add(fhMCAsymmetryDispEta[ie][i]);
1916
1917           fhMCAsymmetryDispPhi[ie][i] = new TH2F (Form("hMCAsymmetryDispPhi_EBin%d_MC%s",ie,pname[i].Data()),
1918                                                   Form("cluster from %s : #sigma^{2}_{#phi #phi} vs A for %d < E < %d GeV",pname[i].Data(),bin[ie],bin[ie+1]),
1919                                                   ssbins,ssmin,ssmax , 200,-1,1);
1920           fhMCAsymmetryDispPhi[ie][i]->SetXTitle("#sigma^{2}_{#phi #phi}");
1921           fhMCAsymmetryDispPhi[ie][i]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1922           outputContainer->Add(fhMCAsymmetryDispPhi[ie][i]);
1923         }
1924       }
1925     }
1926   }
1927
1928   if(fFillPileUpHistograms)
1929   {
1930
1931     TString pileUpName[] = {"SPD","EMCAL","SPDOrEMCAL","SPDAndEMCAL","SPDAndNotEMCAL","EMCALAndNotSPD","NotSPDAndNotEMCAL"} ;
1932
1933     for(Int_t i = 0 ; i < 7 ; i++)
1934     {
1935       fhPtPileUp[i]  = new TH1F(Form("hPtPileUp%s",pileUpName[i].Data()),
1936                                    Form("Selected #pi^{0} (#eta) p_{T} distribution, %s Pile-Up event",pileUpName[i].Data()), nptbins,ptmin,ptmax);
1937       fhPtPileUp[i]->SetXTitle("p_{T} (GeV/c)");
1938       outputContainer->Add(fhPtPileUp[i]);
1939
1940       fhPtCellTimePileUp[i]  = new TH2F(Form("hPtCellTimePileUp%s",pileUpName[i].Data()),
1941                                              Form("Pt vs cell time in cluster, %s Pile-Up event",pileUpName[i].Data()),
1942                                              nptbins,ptmin,ptmax,ntimebins,timemin,timemax);
1943       fhPtCellTimePileUp[i]->SetXTitle("p_{T} (GeV/c)");
1944       fhPtCellTimePileUp[i]->SetYTitle("t_{cell} (ns)");
1945       outputContainer->Add(fhPtCellTimePileUp[i]);
1946
1947       fhPtTimeDiffPileUp[i]  = new TH2F(Form("hPtTimeDiffPileUp%s",pileUpName[i].Data()),
1948                                              Form("Pt vs t_{max}-t_{cell} in cluster, %s Pile-Up event",pileUpName[i].Data()),
1949                                              nptbins,ptmin,ptmax,200,-100,100);
1950       fhPtTimeDiffPileUp[i]->SetXTitle("p_{T} (GeV/c");
1951       fhPtTimeDiffPileUp[i]->SetYTitle("t_{max}-t_{cell} (ns)");
1952       outputContainer->Add(fhPtTimeDiffPileUp[i]);
1953
1954     }
1955
1956     fhTimePtNoCut  = new TH2F ("hTimePt_NoCut","time of cluster vs E of clusters, no cut", nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
1957     fhTimePtNoCut->SetXTitle("p_{T} (GeV/c)");
1958     fhTimePtNoCut->SetYTitle("time (ns)");
1959     outputContainer->Add(fhTimePtNoCut);
1960
1961     fhTimePtSPD  = new TH2F ("hTimePt_SPD","time of cluster vs E of clusters, SPD cut", nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
1962     fhTimePtSPD->SetXTitle("p_{T} (GeV/c)");
1963     fhTimePtSPD->SetYTitle("time (ns)");
1964     outputContainer->Add(fhTimePtSPD);
1965
1966     fhTimePtSPDMulti  = new TH2F ("hTimePt_SPDMulti","time of cluster vs E of clusters, SPD multi cut", nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
1967     fhTimePtSPDMulti->SetXTitle("p_{T} (GeV/c)");
1968     fhTimePtSPDMulti->SetYTitle("time (ns)");
1969     outputContainer->Add(fhTimePtSPDMulti);
1970
1971     fhTimeNPileUpVertSPD  = new TH2F ("hTime_NPileUpVertSPD","time of cluster vs N pile-up SPD vertex", ntimebins,timemin,timemax,50,0,50);
1972     fhTimeNPileUpVertSPD->SetYTitle("# vertex ");
1973     fhTimeNPileUpVertSPD->SetXTitle("time (ns)");
1974     outputContainer->Add(fhTimeNPileUpVertSPD);
1975
1976     fhTimeNPileUpVertTrack  = new TH2F ("hTime_NPileUpVertTracks","time of cluster vs N pile-up Tracks vertex", ntimebins,timemin,timemax, 50,0,50 );
1977     fhTimeNPileUpVertTrack->SetYTitle("# vertex ");
1978     fhTimeNPileUpVertTrack->SetXTitle("time (ns)");
1979     outputContainer->Add(fhTimeNPileUpVertTrack);
1980
1981     fhTimeNPileUpVertContributors  = new TH2F ("hTime_NPileUpVertContributors","time of cluster vs N constributors to pile-up SPD vertex", ntimebins,timemin,timemax,50,0,50);
1982     fhTimeNPileUpVertContributors->SetYTitle("# vertex ");
1983     fhTimeNPileUpVertContributors->SetXTitle("time (ns)");
1984     outputContainer->Add(fhTimeNPileUpVertContributors);
1985
1986     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);
1987     fhTimePileUpMainVertexZDistance->SetYTitle("distance Z (cm) ");
1988     fhTimePileUpMainVertexZDistance->SetXTitle("time (ns)");
1989     outputContainer->Add(fhTimePileUpMainVertexZDistance);
1990
1991     fhTimePileUpMainVertexZDiamond  = new TH2F ("hTime_PileUpMainVertexZDiamond","time of cluster vs distance in Z pile-up SPD vertex - z diamond",ntimebins,timemin,timemax,100,0,50);
1992     fhTimePileUpMainVertexZDiamond->SetYTitle("diamond distance Z (cm) ");
1993     fhTimePileUpMainVertexZDiamond->SetXTitle("time (ns)");
1994     outputContainer->Add(fhTimePileUpMainVertexZDiamond);
1995
1996                 fhPtNPileUpSPDVtx  = new TH2F ("hPt_NPileUpVertSPD","pT of cluster vs N pile-up SPD vertex",
1997                                                                                                                                          nptbins,ptmin,ptmax,20,0,20);
1998                 fhPtNPileUpSPDVtx->SetYTitle("# vertex ");
1999                 fhPtNPileUpSPDVtx->SetXTitle("p_{T} (GeV/c)");
2000                 outputContainer->Add(fhPtNPileUpSPDVtx);
2001
2002                 fhPtNPileUpTrkVtx  = new TH2F ("hPt_NPileUpVertTracks","pT of cluster vs N pile-up Tracks vertex",
2003                                                                                                                                          nptbins,ptmin,ptmax, 20,0,20 );
2004                 fhPtNPileUpTrkVtx->SetYTitle("# vertex ");
2005                 fhPtNPileUpTrkVtx->SetXTitle("p_{T} (GeV/c)");
2006                 outputContainer->Add(fhPtNPileUpTrkVtx);
2007
2008                 fhPtNPileUpSPDVtxTimeCut  = new TH2F ("hPt_NPileUpVertSPD_TimeCut","pT of cluster vs N pile-up SPD vertex, |tof| < 25 ns",
2009                                           nptbins,ptmin,ptmax,20,0,20);
2010                 fhPtNPileUpSPDVtxTimeCut->SetYTitle("# vertex ");
2011                 fhPtNPileUpSPDVtxTimeCut->SetXTitle("p_{T} (GeV/c)");
2012                 outputContainer->Add(fhPtNPileUpSPDVtxTimeCut);
2013
2014                 fhPtNPileUpTrkVtxTimeCut  = new TH2F ("hPt_NPileUpVertTracks_TimeCut","pT of cluster vs N pile-up Tracks vertex, |tof| < 25 ns",
2015                                                                                                                                                                         nptbins,ptmin,ptmax, 20,0,20 );
2016                 fhPtNPileUpTrkVtxTimeCut->SetYTitle("# vertex ");
2017                 fhPtNPileUpTrkVtxTimeCut->SetXTitle("p_{T} (GeV/c)");
2018                 outputContainer->Add(fhPtNPileUpTrkVtxTimeCut);
2019
2020     fhPtNPileUpSPDVtxTimeCut2  = new TH2F ("hPt_NPileUpVertSPD_TimeCut2","pT of cluster vs N pile-up SPD vertex, -25 < tof < 75 ns",
2021                                            nptbins,ptmin,ptmax,20,0,20);
2022                 fhPtNPileUpSPDVtxTimeCut2->SetYTitle("# vertex ");
2023                 fhPtNPileUpSPDVtxTimeCut2->SetXTitle("p_{T} (GeV/c)");
2024                 outputContainer->Add(fhPtNPileUpSPDVtxTimeCut2);
2025
2026                 fhPtNPileUpTrkVtxTimeCut2  = new TH2F ("hPt_NPileUpVertTracks_TimeCut2","pT of cluster vs N pile-up Tracks vertex, -25 < tof < 75 ns",
2027                                            nptbins,ptmin,ptmax, 20,0,20 );
2028                 fhPtNPileUpTrkVtxTimeCut2->SetYTitle("# vertex ");
2029                 fhPtNPileUpTrkVtxTimeCut2->SetXTitle("p_{T} (GeV/c)");
2030                 outputContainer->Add(fhPtNPileUpTrkVtxTimeCut2);
2031
2032   }
2033
2034   //Keep neutral meson selection histograms if requiered
2035   //Setting done in AliNeutralMesonSelection
2036
2037   if(fAnaType!=kSSCalo && GetNeutralMesonSelection())
2038   {
2039     TList * nmsHistos = GetNeutralMesonSelection()->GetCreateOutputObjects() ;
2040
2041     if(GetNeutralMesonSelection()->AreNeutralMesonSelectionHistosKept())
2042       for(Int_t i = 0; i < nmsHistos->GetEntries(); i++) outputContainer->Add(nmsHistos->At(i)) ;
2043
2044     delete nmsHistos;
2045   }
2046
2047   return outputContainer ;
2048
2049 }
2050
2051 //_____________________________________________
2052 Int_t AliAnaPi0EbE::GetMCIndex(const Int_t tag)
2053 {
2054
2055   // Assign mc index depending on MC bit set
2056
2057   if       ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0)  )
2058   {
2059     return kmcPi0 ;
2060   }//pi0
2061   else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta)  )
2062   {
2063     return kmcEta ;
2064   }//eta
2065   else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) &&
2066             GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion) )
2067   {
2068     return kmcConversion ;
2069   }//conversion photon
2070   else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) )
2071   {
2072     return kmcPhoton ;
2073   }//photon   no conversion
2074   else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron))
2075   {
2076     return kmcElectron ;
2077   }//electron
2078   else
2079   {
2080     return kmcHadron ;
2081   }//other particles
2082
2083 }
2084
2085 //__________________________________________________________________
2086 void AliAnaPi0EbE::HasPairSameMCMother(AliAODPWG4Particle * photon1,
2087                                        AliAODPWG4Particle * photon2,
2088                                        Int_t & label, Int_t & tag)
2089 {
2090   // Check the labels of pare in case mother was same pi0 or eta
2091   // Set the new AOD accordingly
2092
2093   Int_t  label1 = photon1->GetLabel();
2094   Int_t  label2 = photon2->GetLabel();
2095
2096   if(label1 < 0 || label2 < 0 ) return ;
2097
2098   //Int_t tag1 = GetMCAnalysisUtils()->CheckOrigin(label1, GetReader());
2099   //Int_t tag2 = GetMCAnalysisUtils()->CheckOrigin(label2, GetReader());
2100   Int_t tag1 = photon1->GetTag();
2101   Int_t tag2 = photon2->GetTag();
2102
2103   if(GetDebug() > 0) printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - Origin of: photon1 %d; photon2 %d \n",tag1, tag2);
2104   if( (GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCPi0Decay) &&
2105        GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCPi0Decay)    ) ||
2106      (GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCEtaDecay) &&
2107       GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCEtaDecay)    )
2108      )
2109   {
2110
2111     //Check if pi0/eta mother is the same
2112     if(GetReader()->ReadStack())
2113     {
2114       if(label1>=0)
2115       {
2116         TParticle * mother1 = GetMCStack()->Particle(label1);//photon in kine tree
2117         label1 = mother1->GetFirstMother();
2118         //mother1 = GetMCStack()->Particle(label1);//pi0
2119       }
2120       if(label2>=0)
2121       {
2122         TParticle * mother2 = GetMCStack()->Particle(label2);//photon in kine tree
2123         label2 = mother2->GetFirstMother();
2124         //mother2 = GetMCStack()->Particle(label2);//pi0
2125       }
2126     } // STACK
2127     else if(GetReader()->ReadAODMCParticles())
2128     {//&& (input > -1)){
2129       if(label1>=0)
2130       {
2131         AliAODMCParticle * mother1 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles())->At(label1);//photon in kine tree
2132         label1 = mother1->GetMother();
2133         //mother1 = GetMCStack()->Particle(label1);//pi0
2134       }
2135       if(label2>=0)
2136       {
2137         AliAODMCParticle * mother2 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles())->At(label2);//photon in kine tree
2138         label2 = mother2->GetMother();
2139         //mother2 = GetMCStack()->Particle(label2);//pi0
2140       }
2141     }// AOD
2142
2143     //printf("mother1 %d, mother2 %d\n",label1,label2);
2144     if( label1 == label2 && label1>=0 )
2145     {
2146
2147       label = label1;
2148
2149       TLorentzVector mom1 = *(photon1->Momentum());
2150       TLorentzVector mom2 = *(photon2->Momentum());
2151
2152       Double_t angle = mom2.Angle(mom1.Vect());
2153       Double_t mass  = (mom1+mom2).M();
2154       Double_t epair = (mom1+mom2).E();
2155
2156       if(GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCPi0Decay))
2157       {
2158         fhMassPairMCPi0 ->Fill(epair,mass);
2159         fhAnglePairMCPi0->Fill(epair,angle);
2160         GetMCAnalysisUtils()->SetTagBit(tag,AliMCAnalysisUtils::kMCPi0);
2161       }
2162       else
2163       {
2164         fhMassPairMCEta ->Fill(epair,mass);
2165         fhAnglePairMCEta->Fill(epair,angle);
2166         GetMCAnalysisUtils()->SetTagBit(tag,AliMCAnalysisUtils::kMCEta);
2167       }
2168
2169     } // same label
2170   } // both from eta or pi0 decay
2171
2172 }
2173
2174 //____________________________________________________________________________
2175 void AliAnaPi0EbE::Init()
2176 {
2177   //Init
2178   //Do some checks
2179   if(fCalorimeter == "PHOS" && !GetReader()->IsPHOSSwitchedOn() && NewOutputAOD()){
2180     printf("AliAnaPi0EbE::Init() - !!STOP: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!\n");
2181     abort();
2182   }
2183   else  if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn() && NewOutputAOD()){
2184     printf("AliAnaPi0EbE::Init() - !!STOP: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!\n");
2185     abort();
2186   }
2187
2188 }
2189
2190 //____________________________________________________________________________
2191 void AliAnaPi0EbE::InitParameters()
2192 {
2193   //Initialize the parameters of the analysis.
2194   AddToHistogramsName("AnaPi0EbE_");
2195
2196   fInputAODGammaConvName = "PhotonsCTS" ;
2197   fAnaType = kIMCalo ;
2198   fCalorimeter = "EMCAL" ;
2199   fMinDist  = 2.;
2200   fMinDist2 = 4.;
2201   fMinDist3 = 5.;
2202
2203   fNLMECutMin[0] = 10.;
2204   fNLMECutMin[1] = 6. ;
2205   fNLMECutMin[2] = 6. ;
2206
2207 }
2208
2209 //__________________________________________________________________
2210 void  AliAnaPi0EbE::MakeAnalysisFillAOD()
2211 {
2212   //Do analysis and fill aods
2213
2214   switch(fAnaType)
2215   {
2216     case kIMCalo:
2217       MakeInvMassInCalorimeter();
2218       break;
2219
2220     case kSSCalo:
2221       MakeShowerShapeIdentification();
2222       break;
2223
2224     case kIMCaloTracks:
2225       MakeInvMassInCalorimeterAndCTS();
2226       break;
2227
2228   }
2229 }
2230
2231 //____________________________________________
2232 void  AliAnaPi0EbE::MakeInvMassInCalorimeter()
2233 {
2234   //Do analysis and fill aods
2235   //Search for the photon decay in calorimeters
2236   //Read photon list from AOD, produced in class AliAnaPhoton
2237   //Check if 2 photons have the mass of the pi0.
2238
2239   TLorentzVector mom1;
2240   TLorentzVector mom2;
2241   TLorentzVector mom ;
2242
2243   Int_t tag   = 0;
2244   Int_t label = 0;
2245
2246   if(!GetInputAODBranch()){
2247     printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - No input calo photons in AOD with name branch < %s >, STOP \n",GetInputAODName().Data());
2248     abort();
2249   }
2250
2251   //Get shower shape information of clusters
2252   TObjArray *clusters = 0;
2253   if     (fCalorimeter=="EMCAL") clusters = GetEMCALClusters();
2254   else if(fCalorimeter=="PHOS")  clusters = GetPHOSClusters() ;
2255
2256   for(Int_t iphoton = 0; iphoton < GetInputAODBranch()->GetEntriesFast()-1; iphoton++){
2257     AliAODPWG4Particle * photon1 =  (AliAODPWG4Particle*) (GetInputAODBranch()->At(iphoton));
2258
2259     //Vertex cut in case of mixed events
2260     Int_t evtIndex1 = 0 ;
2261     if(GetMixedEvent())
2262       evtIndex1 = GetMixedEvent()->EventIndexForCaloCluster(photon1->GetCaloLabel(0)) ;
2263     if(TMath::Abs(GetVertex(evtIndex1)[2]) > GetZvertexCut()) continue ;  //vertex cut
2264     mom1 = *(photon1->Momentum());
2265
2266     //Get original cluster, to recover some information
2267     Int_t iclus = -1;
2268     AliVCluster *cluster1 = FindCluster(clusters,photon1->GetCaloLabel(0),iclus);
2269
2270     if(!cluster1){
2271       printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - First cluster not found\n");
2272       return;
2273     }
2274
2275     for(Int_t jphoton = iphoton+1; jphoton < GetInputAODBranch()->GetEntriesFast(); jphoton++)
2276     {
2277       AliAODPWG4Particle * photon2 =  (AliAODPWG4Particle*) (GetInputAODBranch()->At(jphoton));
2278
2279       Int_t evtIndex2 = 0 ;
2280       if(GetMixedEvent())
2281         evtIndex2 = GetMixedEvent()->EventIndexForCaloCluster(photon2->GetCaloLabel(0)) ;
2282
2283       if(GetMixedEvent() && (evtIndex1 == evtIndex2))
2284         continue ;
2285
2286       if(TMath::Abs(GetVertex(evtIndex2)[2]) > GetZvertexCut()) continue ;  //vertex cut
2287
2288       mom2 = *(photon2->Momentum());
2289
2290       //Get original cluster, to recover some information
2291       Int_t iclus2;
2292       AliVCluster *cluster2 = FindCluster(clusters,photon2->GetCaloLabel(0),iclus2,iclus+1);
2293
2294       if(!cluster2)
2295       {
2296         printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - Second cluster not found\n");
2297         return;
2298       }
2299
2300       Float_t e1    = photon1->E();
2301       Float_t e2    = photon2->E();
2302
2303       //Select clusters with good time window difference
2304       Float_t tof1  = cluster1->GetTOF()*1e9;;
2305       Float_t tof2  = cluster2->GetTOF()*1e9;;
2306       Double_t t12diff = tof1-tof2;
2307       fhEPairDiffTime->Fill(e1+e2,    t12diff);
2308       if(TMath::Abs(t12diff) > GetPairTimeCut()) continue;
2309
2310       //Play with the MC stack if available
2311       if(IsDataMC()) HasPairSameMCMother(photon1, photon2, label, tag) ;
2312
2313       // Check the invariant mass for different selection on the local maxima
2314       // Name of AOD method TO BE FIXED
2315       Int_t nMaxima1 = photon1->GetFiducialArea();
2316       Int_t nMaxima2 = photon2->GetFiducialArea();
2317
2318       Double_t mass  = (mom1+mom2).M();
2319       Double_t epair = (mom1+mom2).E();
2320
2321       if(nMaxima1==nMaxima2)
2322       {
2323         if     (nMaxima1==1) fhMassPairLocMax[0]->Fill(epair,mass);
2324         else if(nMaxima1==2) fhMassPairLocMax[1]->Fill(epair,mass);
2325         else                 fhMassPairLocMax[2]->Fill(epair,mass);
2326       }
2327       else if(nMaxima1==1 || nMaxima2==1)
2328       {
2329         if  (nMaxima1==2 || nMaxima2==2) fhMassPairLocMax[3]->Fill(epair,mass);
2330         else                             fhMassPairLocMax[4]->Fill(epair,mass);
2331       }
2332       else
2333         fhMassPairLocMax[5]->Fill(epair,mass);
2334
2335       // combinations with SS axis cut and NLM cut
2336       if(nMaxima1 == 1 && cluster2->GetM02() > 0.3) fhMassPairLocMax[6]->Fill(epair,mass);
2337       if(nMaxima2 == 1 && cluster1->GetM02() > 0.3) fhMassPairLocMax[6]->Fill(epair,mass);
2338       if(nMaxima1 >  1 && cluster2->GetM02() < 0.3 && cluster2->GetM02()> 0.1 ) fhMassPairLocMax[7]->Fill(epair,mass);
2339       if(nMaxima2 >  1 && cluster1->GetM02() < 0.3 && cluster1->GetM02()> 0.1 ) fhMassPairLocMax[7]->Fill(epair,mass);
2340
2341       //Skip events with too few or too many  NLM
2342       if((nMaxima1 < fNLMCutMin || nMaxima1 > fNLMCutMax) || (nMaxima2 < fNLMCutMin || nMaxima2 > fNLMCutMax)) continue ;
2343
2344       if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - NLM of out of range: cluster1 %d, cluster2 %d \n",nMaxima1, nMaxima2);
2345
2346       //Mass of all pairs
2347       fhMass->Fill(epair,(mom1+mom2).M());
2348
2349       //Select good pair (good phi, pt cuts, aperture and invariant mass)
2350       if(GetNeutralMesonSelection()->SelectPair(mom1, mom2,fCalorimeter))
2351       {
2352         if(GetDebug()>1)
2353           printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - Selected gamma pair: pt %f, phi %f, eta%f \n",(mom1+mom2).Pt(), (mom1+mom2).Phi()*180./3.1416, (mom1+mom2).Eta());
2354
2355         //Fill some histograms about shower shape
2356         if(fFillSelectClHisto && clusters && GetReader()->GetDataType()!=AliCaloTrackReader::kMC)
2357         {
2358           FillSelectedClusterHistograms(cluster1, nMaxima1, photon1->GetTag());
2359           FillSelectedClusterHistograms(cluster2, nMaxima2, photon2->GetTag());
2360         }
2361
2362         // Tag both photons as decay
2363         photon1->SetTagged(kTRUE);
2364         photon2->SetTagged(kTRUE);
2365
2366         fhPtDecay->Fill(photon1->Pt());
2367         fhEDecay ->Fill(photon1->E() );
2368
2369         fhPtDecay->Fill(photon2->Pt());
2370         fhEDecay ->Fill(photon2->E() );
2371
2372         //Create AOD for analysis
2373         mom = mom1+mom2;
2374
2375         //Mass of selected pairs
2376         fhSelectedMass->Fill(epair,mom.M());
2377
2378         // Fill histograms to undertand pile-up before other cuts applied
2379         // Remember to relax time cuts in the reader
2380         FillPileUpHistograms(mom.Pt(),((cluster1->GetTOF()+cluster2->GetTOF())*1e9)/2,cluster1);
2381
2382         AliAODPWG4Particle pi0 = AliAODPWG4Particle(mom);
2383
2384         pi0.SetIdentifiedParticleType(AliCaloPID::kPi0);
2385         pi0.SetDetector(photon1->GetDetector());
2386
2387         // MC
2388         pi0.SetLabel(label);
2389         pi0.SetTag(tag);
2390
2391         //Set the indeces of the original caloclusters
2392         pi0.SetCaloLabel(photon1->GetCaloLabel(0), photon2->GetCaloLabel(0));
2393         //pi0.SetInputFileIndex(input);
2394
2395         AddAODParticle(pi0);
2396
2397       }//pi0
2398
2399     }//2n photon loop
2400
2401   }//1st photon loop
2402
2403   if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - End fill AODs \n");
2404
2405 }
2406
2407 //__________________________________________________
2408 void  AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS()
2409 {
2410   //Do analysis and fill aods
2411   //Search for the photon decay in calorimeters
2412   //Read photon list from AOD, produced in class AliAnaPhoton and AliGammaConversion
2413   //Check if 2 photons have the mass of the pi0.
2414
2415   TLorentzVector mom1;
2416   TLorentzVector mom2;
2417   TLorentzVector mom ;
2418   Int_t tag   = 0;
2419   Int_t label = 0;
2420   Int_t evtIndex = 0;
2421
2422   // Check calorimeter input
2423   if(!GetInputAODBranch()){
2424     printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - No input calo photons in AOD branch with name < %s > , STOP\n",GetInputAODName().Data());
2425     abort();
2426   }
2427
2428   // Get the array with conversion photons
2429   TClonesArray * inputAODGammaConv = (TClonesArray *) GetReader()->GetOutputEvent()->FindListObject(fInputAODGammaConvName);
2430   if(!inputAODGammaConv) {
2431
2432     inputAODGammaConv = (TClonesArray *) GetReader()->GetInputEvent()->FindListObject(fInputAODGammaConvName);
2433
2434     if(!inputAODGammaConv) {
2435       printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - No input gamma conversions in AOD branch with name < %s >\n",fInputAODGammaConvName.Data());
2436
2437       return;
2438     }
2439   }
2440
2441   //Get shower shape information of clusters
2442   TObjArray *clusters = 0;
2443   if     (fCalorimeter=="EMCAL") clusters = GetEMCALClusters();
2444   else if(fCalorimeter=="PHOS")  clusters = GetPHOSClusters() ;
2445
2446   Int_t nCTS  = inputAODGammaConv->GetEntriesFast();
2447   Int_t nCalo = GetInputAODBranch()->GetEntriesFast();
2448   if(nCTS<=0 || nCalo <=0)
2449   {
2450     if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - nCalo %d, nCTS %d, cannot loop\n",nCalo,nCTS);
2451     return;
2452   }
2453
2454   if(GetDebug() > 1)
2455     printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - Number of conversion photons %d\n",nCTS);
2456
2457   // Do the loop, first calo, second CTS
2458   for(Int_t iphoton = 0; iphoton < GetInputAODBranch()->GetEntriesFast(); iphoton++){
2459     AliAODPWG4Particle * photon1 =  (AliAODPWG4Particle*) (GetInputAODBranch()->At(iphoton));
2460     mom1 = *(photon1->Momentum());
2461
2462     //Get original cluster, to recover some information
2463     Int_t iclus = -1;
2464     AliVCluster *cluster = FindCluster(clusters,photon1->GetCaloLabel(0),iclus);
2465
2466     for(Int_t jphoton = 0; jphoton < nCTS; jphoton++){
2467       AliAODPWG4Particle * photon2 =  (AliAODPWG4Particle*) (inputAODGammaConv->At(jphoton));
2468       if(GetMixedEvent())
2469         evtIndex = GetMixedEvent()->EventIndexForCaloCluster(photon2->GetCaloLabel(0)) ;
2470       if(TMath::Abs(GetVertex(evtIndex)[2]) > GetZvertexCut()) continue ;  //vertex cut
2471
2472       mom2 = *(photon2->Momentum());
2473
2474       Double_t mass  = (mom1+mom2).M();
2475       Double_t epair = (mom1+mom2).E();
2476
2477       Int_t nMaxima = photon1->GetFiducialArea();
2478       if     (nMaxima==1) fhMassPairLocMax[0]->Fill(epair,mass);
2479       else if(nMaxima==2) fhMassPairLocMax[1]->Fill(epair,mass);
2480       else                fhMassPairLocMax[2]->Fill(epair,mass);
2481
2482       if(nMaxima < fNLMCutMin || nMaxima > fNLMCutMax) continue ;
2483       if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - NLM %d of out of range \n",nMaxima);
2484
2485       //Play with the MC stack if available
2486       if(IsDataMC())
2487       {
2488         Int_t   label2 = photon2->GetLabel();
2489         if(label2 >= 0 )photon2->SetTag(GetMCAnalysisUtils()->CheckOrigin(label2, GetReader()));
2490
2491         HasPairSameMCMother(photon1, photon2, label, tag) ;
2492       }
2493
2494       //Mass of selected pairs
2495       fhMass->Fill(epair,(mom1+mom2).M());
2496
2497       //Select good pair (good phi, pt cuts, aperture and invariant mass)
2498       if(GetNeutralMesonSelection()->SelectPair(mom1, mom2,fCalorimeter))
2499       {
2500         if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - Selected gamma pair: pt %f, phi %f, eta%f\n",(mom1+mom2).Pt(), (mom1+mom2).Phi()*180./3.1416, (mom1+mom2).Eta());
2501
2502         //Fill some histograms about shower shape
2503         if(fFillSelectClHisto && cluster && GetReader()->GetDataType()!=AliCaloTrackReader::kMC)
2504         {
2505           FillSelectedClusterHistograms(cluster, nMaxima, photon1->GetTag());
2506         }
2507
2508         // Tag both photons as decay
2509         photon1->SetTagged(kTRUE);
2510         photon2->SetTagged(kTRUE);
2511
2512         fhPtDecay->Fill(photon1->Pt());
2513         fhEDecay ->Fill(photon1->E() );
2514
2515         //Create AOD for analysis
2516
2517         mom = mom1+mom2;
2518
2519         //Mass of selected pairs
2520         fhSelectedMass->Fill(epair,mom.M());
2521
2522         // Fill histograms to undertand pile-up before other cuts applied
2523         // Remember to relax time cuts in the reader
2524         if(cluster) FillPileUpHistograms(mom.Pt(),cluster->GetTOF()*1e9,cluster);
2525
2526         AliAODPWG4Particle pi0 = AliAODPWG4Particle(mom);
2527
2528         pi0.SetIdentifiedParticleType(AliCaloPID::kPi0);
2529         pi0.SetDetector(photon1->GetDetector());
2530
2531         // MC
2532         pi0.SetLabel(label);
2533         pi0.SetTag(tag);
2534
2535         //Set the indeces of the original tracks or caloclusters
2536         pi0.SetCaloLabel(photon1->GetCaloLabel(0), -1);
2537         pi0.SetTrackLabel(photon2->GetTrackLabel(0), photon2->GetTrackLabel(1));
2538         //pi0.SetInputFileIndex(input);
2539
2540         AddAODParticle(pi0);
2541
2542       }//pi0
2543     }//2n photon loop
2544
2545   }//1st photon loop
2546
2547   if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - End fill AODs \n");
2548
2549 }
2550
2551
2552 //_________________________________________________
2553 void  AliAnaPi0EbE::MakeShowerShapeIdentification()
2554 {
2555   //Search for pi0 in fCalorimeter with shower shape analysis
2556
2557   TObjArray * pl        = 0x0;
2558   AliVCaloCells * cells = 0x0;
2559   //Select the Calorimeter of the photon
2560   if      (fCalorimeter == "PHOS" )
2561   {
2562     pl    = GetPHOSClusters();
2563     cells = GetPHOSCells();
2564   }
2565   else if (fCalorimeter == "EMCAL")
2566   {
2567     pl    = GetEMCALClusters();
2568     cells = GetEMCALCells();
2569   }
2570
2571   if(!pl)
2572   {
2573     Info("MakeShowerShapeIdentification","TObjArray with %s clusters is NULL!\n",fCalorimeter.Data());
2574     return;
2575   }
2576
2577   TLorentzVector mom ;
2578   for(Int_t icalo = 0; icalo < pl->GetEntriesFast(); icalo++)
2579   {
2580     AliVCluster * calo = (AliVCluster*) (pl->At(icalo));
2581
2582     Int_t evtIndex = 0 ;
2583     if (GetMixedEvent())
2584     {
2585       evtIndex=GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ;
2586     }
2587
2588     if(TMath::Abs(GetVertex(evtIndex)[2]) > GetZvertexCut()) continue ;  //vertex cut
2589
2590     //Get Momentum vector,
2591     Double_t vertex[]={0,0,0};
2592     if(GetReader()->GetDataType() != AliCaloTrackReader::kMC)
2593     {
2594       calo->GetMomentum(mom,GetVertex(evtIndex)) ;
2595     }//Assume that come from vertex in straight line
2596     else
2597     {
2598       calo->GetMomentum(mom,vertex) ;
2599     }
2600
2601     //If too small or big pt, skip it
2602     if(mom.E() < GetMinEnergy() || mom.E() > GetMaxEnergy() ) continue ;
2603
2604     //Check acceptance selection
2605     if(IsFiducialCutOn())
2606     {
2607       Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter) ;
2608       if(! in ) continue ;
2609     }
2610
2611     if(GetDebug() > 1)
2612       printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - FillAOD: Min pt cut and fiducial cut passed: pt %3.2f, phi %2.2f, eta %1.2f\n",mom.Pt(),mom.Phi(),mom.Eta());
2613
2614     //Check Distance to Bad channel, set bit.
2615     Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
2616     if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;
2617     if(distBad < fMinDist) //In bad channel (PHOS cristal size 2.2x2.2 cm)
2618       continue ;
2619
2620     if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - FillAOD: Bad channel cut passed %4.2f\n",distBad);
2621
2622     //.......................................
2623     // TOF cut, BE CAREFUL WITH THIS CUT
2624     Double_t tof = calo->GetTOF()*1e9;
2625     if(tof < fTimeCutMin || tof > fTimeCutMax) continue ;
2626
2627     //Play with the MC stack if available
2628     //Check origin of the candidates
2629     Int_t tag   = 0 ;
2630     if(IsDataMC())
2631     {
2632       tag = GetMCAnalysisUtils()->CheckOrigin(calo->GetLabels(),calo->GetNLabels(),GetReader());
2633       //GetMCAnalysisUtils()->CheckMultipleOrigin(calo->GetLabels(),calo->GetNLabels(), GetReader(), aodpi0.GetInputFileIndex(), tag);
2634       if(GetDebug() > 0) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - Origin of candidate %d\n",tag);
2635     }
2636
2637     //Skip matched clusters with tracks
2638     if(fRejectTrackMatch && IsTrackMatched(calo, GetReader()->GetInputEvent()))
2639     {
2640       FillRejectedClusterHistograms(mom,tag);
2641       continue ;
2642     }
2643
2644     //Check PID
2645     //PID selection or bit setting
2646     Int_t    nMaxima = 0 ;
2647     Double_t mass    = 0 , angle = 0;
2648     TLorentzVector    l1, l2;
2649     Int_t    absId1 = -1; Int_t absId2 = -1;
2650
2651     Int_t idPartType = GetCaloPID()->GetIdentifiedParticleTypeFromClusterSplitting(calo,cells,GetCaloUtils(),
2652                                                                                    GetVertex(evtIndex),nMaxima,
2653                                                                                    mass,angle,l1,l2,absId1,absId2) ;
2654
2655     if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - PDG of identified particle %d\n",idPartType);
2656
2657
2658     //Skip events with too few or too many  NLM
2659     if(nMaxima < fNLMCutMin || nMaxima > fNLMCutMax)
2660     {
2661       FillRejectedClusterHistograms(mom,tag);
2662       continue ;
2663     }
2664
2665     if(nMaxima == 1 && fNLMECutMin[0] > mom.E()) continue;
2666     if(nMaxima == 2 && fNLMECutMin[1] > mom.E()) continue;
2667     if(nMaxima >  2 && fNLMECutMin[2] > mom.E()) continue;
2668
2669     if(GetDebug() > 1)
2670       printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - NLM %d accepted \n",nMaxima);
2671
2672     Float_t e1 = l1.Energy();
2673     Float_t e2 = l2.Energy();
2674     TLorentzVector l12 = l1+l2;
2675     Float_t ptSplit = l12.Pt();
2676     Float_t  eSplit = e1+e2;
2677     Int_t mcIndex = GetMCIndex(tag);
2678
2679     //mass of all clusters
2680     fhMass       ->Fill(mom.E(),mass);
2681     fhMassPt     ->Fill(mom.Pt(),mass);
2682     fhMassSplitPt->Fill(ptSplit,mass);
2683
2684     if(IsDataMC())
2685     {
2686       fhMCMassPt[mcIndex]     ->Fill(mom.Pt(),mass);
2687       fhMCMassSplitPt[mcIndex]->Fill(ptSplit,mass);
2688     }
2689
2690     // Asymmetry of all clusters
2691     Float_t asy =-10;
2692
2693     if(e1+e2 > 0) asy = (e1-e2) / (e1+e2);
2694     fhAsymmetry->Fill(mom.E(),asy);
2695
2696
2697     if(IsDataMC())
2698     {
2699       fhMCEAsymmetry[mcIndex]->Fill(mom.E(),asy);
2700     }
2701
2702     // If cluster does not pass pid, not pi0/eta, skip it.
2703     if     (GetOutputAODName().Contains("Pi0") && idPartType != AliCaloPID::kPi0)
2704     {
2705       if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - Cluster is not Pi0\n");
2706       FillRejectedClusterHistograms(mom,tag);
2707       continue ;
2708     }
2709
2710     else if(GetOutputAODName().Contains("Eta") && idPartType != AliCaloPID::kEta)
2711     {
2712       if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - Cluster is not Eta\n");
2713       FillRejectedClusterHistograms(mom,tag);
2714       continue ;
2715     }
2716
2717     if(GetDebug() > 1)
2718       printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - Pi0/Eta selection cuts passed: pT %3.2f, pdg %d\n",
2719              mom.Pt(), idPartType);
2720
2721     //Mass and asymmetry of selected pairs
2722     fhSelectedAsymmetry  ->Fill(mom.E() ,asy );
2723     fhSelectedMass       ->Fill(mom.E() ,mass);
2724     fhSelectedMassPt     ->Fill(mom.Pt(),mass);
2725     fhSelectedMassSplitPt->Fill(ptSplit ,mass);
2726
2727     fhSplitE        ->Fill( eSplit);
2728     fhSplitPt       ->Fill(ptSplit);
2729     Float_t phi = mom.Phi();
2730     if(phi<0) phi+=TMath::TwoPi();
2731     fhSplitPtPhi    ->Fill(ptSplit,phi);
2732     fhSplitPtEta    ->Fill(ptSplit,mom.Eta());
2733     fhNLocMaxSplitPt->Fill(ptSplit ,nMaxima);
2734     fhNLocMaxPt     ->Fill(mom.Pt(),nMaxima);
2735
2736     //Check split-clusters with good time window difference
2737     Double_t tof1  = cells->GetCellTime(absId1);
2738     GetCaloUtils()->RecalibrateCellTime(tof1, fCalorimeter, absId1,GetReader()->GetInputEvent()->GetBunchCrossNumber());
2739     tof1*=1.e9;
2740
2741     Double_t tof2  = cells->GetCellTime(absId2);
2742     GetCaloUtils()->RecalibrateCellTime(tof2, fCalorimeter, absId2,GetReader()->GetInputEvent()->GetBunchCrossNumber());
2743     tof2*=1.e9;
2744
2745     Double_t t12diff = tof1-tof2;
2746     fhEPairDiffTime->Fill(e1+e2,    t12diff);
2747
2748     if(IsDataMC())
2749     {
2750       fhMCSplitE        [mcIndex]->Fill( eSplit);
2751       fhMCSplitPt       [mcIndex]->Fill(ptSplit);
2752       fhMCSplitPtPhi    [mcIndex]->Fill(ptSplit,phi);
2753       fhMCSplitPtEta    [mcIndex]->Fill(ptSplit,mom.Eta());
2754       fhMCNLocMaxSplitPt[mcIndex]->Fill(ptSplit ,nMaxima);
2755       fhMCNLocMaxPt     [mcIndex]->Fill(mom.Pt(),nMaxima);
2756
2757       fhMCSelectedMassPt     [mcIndex]->Fill(mom.Pt(),mass);
2758       fhMCSelectedMassSplitPt[mcIndex]->Fill(ptSplit,mass);
2759     }
2760
2761     //-----------------------
2762     //Create AOD for analysis
2763
2764     AliAODPWG4Particle aodpi0 = AliAODPWG4Particle(mom);
2765     aodpi0.SetLabel(calo->GetLabel());
2766
2767     //Set the indeces of the original caloclusters
2768     aodpi0.SetCaloLabel(calo->GetID(),-1);
2769     aodpi0.SetDetector(fCalorimeter);
2770
2771     if     (distBad > fMinDist3) aodpi0.SetDistToBad(2) ;
2772     else if(distBad > fMinDist2) aodpi0.SetDistToBad(1) ;
2773     else                         aodpi0.SetDistToBad(0) ;
2774
2775     // Check if cluster is pi0 via cluster splitting
2776     aodpi0.SetIdentifiedParticleType(idPartType);
2777
2778     // Add number of local maxima to AOD, method name in AOD to be FIXED
2779     aodpi0.SetFiducialArea(nMaxima);
2780
2781     aodpi0.SetTag(tag);
2782
2783     //Fill some histograms about shower shape
2784     if(fFillSelectClHisto && GetReader()->GetDataType()!=AliCaloTrackReader::kMC)
2785     {
2786       FillSelectedClusterHistograms(calo, nMaxima, tag, asy);
2787     }
2788
2789     // Fill histograms to undertand pile-up before other cuts applied
2790     // Remember to relax time cuts in the reader
2791     Double_t tofcluster   = calo->GetTOF()*1e9;
2792     Double_t tofclusterUS = TMath::Abs(tofcluster);
2793
2794     FillPileUpHistograms(aodpi0.Pt(),tofcluster,calo);
2795
2796     Int_t id = GetReader()->GetTriggerClusterId();
2797     if(fFillEMCALBCHistograms && fCalorimeter=="EMCAL" && id >=0 )
2798     {
2799       Float_t phicluster = aodpi0.Phi();
2800       if(phicluster < 0) phicluster+=TMath::TwoPi();
2801
2802       if(calo->E() > 2)
2803       {
2804         if      (tofclusterUS < 25) fhEtaPhiEMCALBC0->Fill(aodpi0.Eta(), phicluster);
2805         else if (tofclusterUS < 75) fhEtaPhiEMCALBC1->Fill(aodpi0.Eta(), phicluster);
2806         else                        fhEtaPhiEMCALBCN->Fill(aodpi0.Eta(), phicluster);
2807       }
2808
2809       Int_t bc = GetReader()->GetTriggerClusterBC();
2810       if(TMath::Abs(bc) < 6  && !GetReader()->IsBadCellTriggerEvent() && !GetReader()->IsExoticEvent() )
2811       {
2812         if(GetReader()->IsTriggerMatched())
2813         {
2814           if(calo->E() > 2) fhEtaPhiTriggerEMCALBC[bc+5]->Fill(aodpi0.Eta(), phicluster);
2815           fhTimeTriggerEMCALBC[bc+5]->Fill(calo->E(), tofcluster);
2816           if(GetReader()->IsPileUpFromSPD()) fhTimeTriggerEMCALBCPileUpSPD[bc+5]->Fill(calo->E(), tofcluster);
2817         }
2818         else
2819         {
2820           if(calo->E() > 2) fhEtaPhiTriggerEMCALBCUM[bc+5]->Fill(aodpi0.Eta(), phicluster);
2821           fhTimeTriggerEMCALBCUM[bc+5]->Fill(calo->E(), tofcluster);
2822
2823           if(bc==0)
2824           {
2825             if(GetReader()->IsTriggerMatchedOpenCuts(0)) fhTimeTriggerEMCALBC0UMReMatchOpenTime   ->Fill(calo->E(), tofcluster);
2826             if(GetReader()->IsTriggerMatchedOpenCuts(1)) fhTimeTriggerEMCALBC0UMReMatchCheckNeigh ->Fill(calo->E(), tofcluster);
2827             if(GetReader()->IsTriggerMatchedOpenCuts(2)) fhTimeTriggerEMCALBC0UMReMatchBoth       ->Fill(calo->E(), tofcluster);
2828           }
2829          }
2830       }
2831       else if(TMath::Abs(bc) >= 6)
2832         printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - Trigger BC not expected = %d\n",bc);
2833     }
2834
2835     //Add AOD with pi0 object to aod branch
2836     AddAODParticle(aodpi0);
2837
2838   }//loop
2839
2840   if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - End fill AODs \n");
2841
2842 }
2843 //______________________________________________
2844 void  AliAnaPi0EbE::MakeAnalysisFillHistograms()
2845 {
2846   //Do analysis and fill histograms
2847
2848   if(!GetOutputAODBranch())
2849   {
2850     printf("AliAnaPi0EbE::MakeAnalysisFillHistograms()  - No output pi0 in AOD branch with name < %s >,STOP \n",GetOutputAODName().Data());
2851     abort();
2852   }
2853   //Loop on stored AOD pi0
2854   Int_t naod = GetOutputAODBranch()->GetEntriesFast();
2855   if(GetDebug() > 0) printf("AliAnaPi0EbE::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod);
2856
2857   Float_t cen = GetEventCentrality();
2858   Float_t ep  = GetEventPlaneAngle();
2859
2860   for(Int_t iaod = 0; iaod < naod ; iaod++)
2861   {
2862
2863     AliAODPWG4Particle* pi0 =  (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
2864     Int_t pdg = pi0->GetIdentifiedParticleType();
2865
2866     if(IsCaloPIDOn() && pdg != AliCaloPID::kPi0) continue;
2867
2868     //Fill pi0 histograms
2869     Float_t ener  = pi0->E();
2870     Float_t pt    = pi0->Pt();
2871     Float_t phi   = pi0->Phi();
2872     if(phi < 0) phi+=TMath::TwoPi();
2873     Float_t eta = pi0->Eta();
2874
2875     fhPt     ->Fill(pt  );
2876     fhE      ->Fill(ener);
2877
2878     fhEEta   ->Fill(ener,eta);
2879     fhEPhi   ->Fill(ener,phi);
2880     fhPtEta  ->Fill(pt  ,eta);
2881     fhPtPhi  ->Fill(pt  ,phi);
2882     fhEtaPhi ->Fill(eta ,phi);
2883
2884     fhPtCentrality ->Fill(pt,cen) ;
2885     fhPtEventPlane ->Fill(pt,ep ) ;
2886
2887     if(IsDataMC())
2888     {
2889       Int_t tag     = pi0->GetTag();
2890       Int_t mcIndex = GetMCIndex(tag);
2891
2892       fhMCE  [mcIndex] ->Fill(ener);
2893       fhMCPt [mcIndex] ->Fill(pt);
2894       fhMCPhi[mcIndex] ->Fill(pt,phi);
2895       fhMCEta[mcIndex] ->Fill(pt,eta);
2896
2897       fhMCPtCentrality[mcIndex]->Fill(pt,cen);
2898
2899       if((mcIndex==kmcPhoton || mcIndex==kmcPi0 || mcIndex==kmcEta) && fAnaType==kSSCalo)
2900       {
2901         Float_t efracMC = 0;
2902         Int_t label = pi0->GetLabel();
2903
2904         Bool_t ok = kFALSE;
2905         TLorentzVector mom   = GetMCAnalysisUtils()->GetMother(label,GetReader(),ok);
2906         if(!ok) continue;
2907
2908         if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0))
2909         {
2910           TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(label,111,GetReader(),ok);
2911           if(grandmom.E() > 0 && ok)
2912           {
2913             efracMC =  grandmom.E()/ener;
2914             fhMCPi0PtGenRecoFraction ->Fill(pt,efracMC);
2915           }
2916         }
2917         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay))
2918         {
2919           fhMCPi0DecayPt->Fill(pt);
2920           TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(label,111,GetReader(),ok);
2921           if(grandmom.E() > 0 && ok)
2922           {
2923             efracMC =  mom.E()/grandmom.E();
2924             fhMCPi0DecayPtFraction ->Fill(pt,efracMC);
2925           }
2926         }
2927         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta))
2928         {
2929           TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(label,221,GetReader(),ok);
2930           if(grandmom.E() > 0 && ok)
2931           {
2932             efracMC =  grandmom.E()/ener;
2933             fhMCEtaPtGenRecoFraction ->Fill(pt,efracMC);
2934           }
2935         }
2936         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay))
2937         {
2938           fhMCEtaDecayPt->Fill(pt);
2939           TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(label,221,GetReader(),ok);
2940           if(grandmom.E() > 0 && ok)
2941           {
2942             efracMC =  mom.E()/grandmom.E();
2943             fhMCEtaDecayPtFraction ->Fill(pt,efracMC);
2944           }
2945         }
2946         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay))
2947         {
2948           fhMCOtherDecayPt->Fill(pt);
2949         }
2950
2951       }
2952
2953     }//Histograms with MC
2954
2955   }// aod loop
2956
2957 }
2958
2959 //__________________________________________________________________
2960 void AliAnaPi0EbE::Print(const Option_t * opt) const
2961 {
2962   //Print some relevant parameters set for the analysis
2963   if(! opt)
2964     return;
2965
2966   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
2967   AliAnaCaloTrackCorrBaseClass::Print("");
2968   printf("Analysis Type = %d \n",  fAnaType) ;
2969   if(fAnaType == kSSCalo)
2970   {
2971     printf("Calorimeter            =     %s\n", fCalorimeter.Data()) ;
2972     printf("Min Distance to Bad Channel   = %2.1f\n",fMinDist);
2973     printf("Min Distance to Bad Channel 2 = %2.1f\n",fMinDist2);
2974     printf("Min Distance to Bad Channel 3 = %2.1f\n",fMinDist3);
2975   }
2976   printf("    \n") ;
2977
2978 }
2979
2980