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