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