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