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),                  fCalorimeter(""),
  53 fMinDist(0.),fMinDist2(0.),         fMinDist3(0.),
  54 fNLMCutMin(-1),                     fNLMCutMax(10),
  55 fTimeCutMin(-10000),                fTimeCutMax(10000),
  56 fRejectTrackMatch(kTRUE),           fSelectIsolatedDecay(kFALSE),
  57 fFillPileUpHistograms(0),
  58 fFillWeightHistograms(kFALSE),      fFillTMHisto(0),
  59 fFillSelectClHisto(0),              fFillOnlySimpleSSHisto(1),
  60 fFillEMCALBCHistograms(0),          fFillHighMultHistograms(0),
  61 fFillAllNLMHistograms(0),
  62 fInputAODGammaConvName(""),
  63 fCheckSplitDistToBad(0),
  64 // Histograms
  65 fhPt(0),                            fhE(0),
  66 fhPtEta(0),                         fhPtPhi(0),                         fhEtaPhi(0),
  67 fhEtaPhiEMCALBC0(0),                fhEtaPhiEMCALBC1(0),                fhEtaPhiEMCALBCN(0),
  68 fhTimeTriggerEMCALBC0UMReMatchOpenTime(0),
  69 fhTimeTriggerEMCALBC0UMReMatchCheckNeigh(0),
  70 fhTimeTriggerEMCALBC0UMReMatchBoth(0),
  71 fhPtCentrality(),                   fhPtEventPlane(0),                  fhMCPtCentrality(),
  72 fhPtReject(0),                      fhEReject(0),
  73 fhPtEtaReject(0),                   fhPtPhiReject(0),                   fhEtaPhiReject(0),
  74 fhMass(0),                          fhMassPt(0),                        fhMassSplitPt(0),
  75 fhSelectedMass(),                   fhSelectedMassPt(0),                fhSelectedMassSplitPt(0),
  76 fhMassNoOverlap(0),                 fhMassPtNoOverlap(0),               fhMassSplitPtNoOverlap(0),
  77 fhSelectedMassNoOverlap(0),         fhSelectedMassPtNoOverlap(0),       fhSelectedMassSplitPtNoOverlap(0),
  78 fhMCPi0PtRecoPtPrim(0),                       fhMCEtaPtRecoPtPrim(0),
  79 fhMCPi0PtRecoPtPrimNoOverlap(0),              fhMCEtaPtRecoPtPrimNoOverlap(0),
  80 fhMCPi0SplitPtRecoPtPrim(0),                  fhMCEtaSplitPtRecoPtPrim(0),
  81 fhMCPi0SplitPtRecoPtPrimNoOverlap(0),         fhMCEtaSplitPtRecoPtPrimNoOverlap(0),
  82 fhMCPi0SelectedPtRecoPtPrim(0),               fhMCEtaSelectedPtRecoPtPrim(0),
  83 fhMCPi0SelectedPtRecoPtPrimNoOverlap(0),      fhMCEtaSelectedPtRecoPtPrimNoOverlap(0),
  84 fhMCPi0SelectedSplitPtRecoPtPrim(0),          fhMCEtaSelectedSplitPtRecoPtPrim(0),
  85 fhMCPi0SelectedSplitPtRecoPtPrimNoOverlap(0), fhMCEtaSelectedSplitPtRecoPtPrimNoOverlap(0),
  86 fhAsymmetry(0),                     fhSelectedAsymmetry(0),
  87 fhSplitE(0),                        fhSplitPt(0),
  88 fhSplitPtEta(0),                    fhSplitPtPhi(0),
  89 fhNLocMaxSplitPt(0),
  90 fhPtDecay(0),
  91 // Shower shape histos
  92 fhPtDispersion(0),                  fhPtLambda0(0),                     fhPtLambda0NoSplitCut(0),
  93 fhPtLambda1(0),                     fhPtLambda0NoTRD(0),                fhPtLambda0FracMaxCellCut(0),
  94 fhPtFracMaxCell(0),                 fhPtFracMaxCellNoTRD(0),
  95 fhPtNCells(0),                      fhPtTime(0),                        fhEPairDiffTime(0),
  96 fhPtDispEta(0),                     fhPtDispPhi(0),
  97 fhPtSumEta(0),                      fhPtSumPhi(0),                      fhPtSumEtaPhi(0),
  98 fhPtDispEtaPhiDiff(0),              fhPtSphericity(0),
  99
 100 // MC histos
 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(!fFillPileUpHistograms) 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(fCalorimeter == "EMCAL") 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, fCalorimeter, 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(fCalorimeter == "EMCAL" && !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(fCalorimeter=="EMCAL" &&  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(fCalorimeter=="EMCAL" &&  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(fCalorimeter=="EMCAL" && GetFirstSMCoveredByTRD() >= 0 &&
 719        GetModuleNumber(cluster) < GetFirstSMCoveredByTRD() )
 720       fhMCPtLambda0NoTRD[mcIndex]->Fill(pt, l0  );
 721
 722     if(fCalorimeter == "EMCAL" && !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(fCalorimeter == "EMCAL") 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,fCalorimeter, 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,fCalorimeter, 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(fCalorimeter=="EMCAL"){
 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;",fCalorimeter.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 [] = {"#gamma (#pi^{0})", "#gamma (#eta)", "#gamma (other)",  "#pi^{0}", "#eta", "#gamma (direct)","e^{#pm}"  , "hadron/other combinations"};
 930   TString pname [] = {"Pi0Decay"        , "EtaDecay"     , "OtherDecay"    ,  "Pi0"    , "Eta" , "Photon"         , "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(fCalorimeter=="EMCAL" && 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(fFillHighMultHistograms)
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       for(Int_t ipart = 0; ipart < fgkNmcTypes; ipart++)
1218       {
1219         fhMCPtDecay[ipart]  = new TH1F(Form("hPtDecay_MC%s",pname[ipart].Data()),
1220                                        Form("Selected  #pi^{0} (#eta) decay photons, from MC %s",ptype[ipart].Data()),
1221                                        nptbins,ptmin,ptmax);
1222         fhMCPtDecay[ipart]->SetYTitle("#it{N}");
1223         fhMCPtDecay[ipart]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1224         outputContainer->Add(fhMCPtDecay[ipart]) ;
1225       }
1226     }
1227   }
1228
1229   ////////
1230
1231   if( fFillSelectClHisto )
1232   {
1233     fhPtLambda0  = new TH2F
1234     ("hPtLambda0","Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #lambda_{0}",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1235     fhPtLambda0->SetYTitle("#lambda_{0}^{2}");
1236     fhPtLambda0->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1237     outputContainer->Add(fhPtLambda0) ;
1238
1239     fhPtLambda1  = new TH2F
1240     ("hPtLambda1","Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #lambda_{1}",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1241     fhPtLambda1->SetYTitle("#lambda_{1}^{2}");
1242     fhPtLambda1->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1243     outputContainer->Add(fhPtLambda1) ;
1244
1245     if(fCalorimeter=="EMCAL" &&  GetFirstSMCoveredByTRD() >=0 )
1246     {
1247       fhPtLambda0NoTRD  = new TH2F
1248       ("hPtLambda0NoTRD","Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #lambda_{0}, not behind TRD",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1249       fhPtLambda0NoTRD->SetYTitle("#lambda_{0}^{2}");
1250       fhPtLambda0NoTRD->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1251       outputContainer->Add(fhPtLambda0NoTRD) ;
1252
1253       if(!fFillOnlySimpleSSHisto)
1254       {
1255         fhPtFracMaxCellNoTRD  = new TH2F
1256         ("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);
1257         fhPtFracMaxCellNoTRD->SetYTitle("Fraction");
1258         fhPtFracMaxCellNoTRD->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1259         outputContainer->Add(fhPtFracMaxCellNoTRD) ;
1260       }
1261     }
1262
1263     if(!fFillOnlySimpleSSHisto)
1264     {
1265       fhPtDispersion  = new TH2F
1266       ("hPtDispersion","Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs dispersion",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1267       fhPtDispersion->SetYTitle("D^{2}");
1268       fhPtDispersion->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1269       outputContainer->Add(fhPtDispersion) ;
1270
1271       fhPtLambda0FracMaxCellCut  = new TH2F
1272       ("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);
1273       fhPtLambda0FracMaxCellCut->SetYTitle("#lambda_{0}^{2}");
1274       fhPtLambda0FracMaxCellCut->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1275       outputContainer->Add(fhPtLambda0FracMaxCellCut) ;
1276
1277       fhPtFracMaxCell  = new TH2F
1278       ("hPtFracMaxCell","Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #lambda_{0}, Max cell fraction of energy",nptbins,ptmin,ptmax,100,0,1);
1279       fhPtFracMaxCell->SetYTitle("Fraction");
1280       fhPtFracMaxCell->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1281       outputContainer->Add(fhPtFracMaxCell) ;
1282
1283       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);
1284       fhPtDispEta->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1285       fhPtDispEta->SetYTitle("#sigma^{2}_{#eta #eta}");
1286       outputContainer->Add(fhPtDispEta);
1287
1288       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);
1289       fhPtDispPhi->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1290       fhPtDispPhi->SetYTitle("#sigma^{2}_{#phi #phi}");
1291       outputContainer->Add(fhPtDispPhi);
1292
1293       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);
1294       fhPtSumEta->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1295       fhPtSumEta->SetYTitle("#delta^{2}_{#eta #eta}");
1296       outputContainer->Add(fhPtSumEta);
1297
1298       fhPtSumPhi  = new TH2F ("hPtSumPhi","#sigma^{2}_{#phi #phi} = #Sigma w_{i}(#phi_{i})^{2}/ #Sigma w_{i} - <#phi>^{2} vs #it{p}_{T}",
1299                               nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
1300       fhPtSumPhi->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1301       fhPtSumPhi->SetYTitle("#delta^{2}_{#phi #phi}");
1302       outputContainer->Add(fhPtSumPhi);
1303
1304       fhPtSumEtaPhi  = new TH2F ("hPtSumEtaPhi","#delta^{2}_{#eta #phi} = #Sigma w_{i}(#phi_{i} #eta_{i} ) / #Sigma w_{i} - <#phi><#eta> vs #it{p}_{T}",
1305                                  nptbins,ptmin,ptmax, 2*ssbins,-ssmax,ssmax);
1306       fhPtSumEtaPhi->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1307       fhPtSumEtaPhi->SetYTitle("#delta^{2}_{#eta #phi}");
1308       outputContainer->Add(fhPtSumEtaPhi);
1309
1310       fhPtDispEtaPhiDiff  = new TH2F ("hPtDispEtaPhiDiff","#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta} vs #it{p}_{T}",
1311                                       nptbins,ptmin,ptmax,200, -10,10);
1312       fhPtDispEtaPhiDiff->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1313       fhPtDispEtaPhiDiff->SetYTitle("#sigma^{2}_{#phi #phi}-#sigma^{2}_{#eta #eta}");
1314       outputContainer->Add(fhPtDispEtaPhiDiff);
1315
1316       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})",
1317                                   nptbins,ptmin,ptmax, 200, -1,1);
1318       fhPtSphericity->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1319       fhPtSphericity->SetYTitle("s = (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi})");
1320       outputContainer->Add(fhPtSphericity);
1321
1322       for(Int_t i = 0; i < 7; i++)
1323       {
1324         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]),
1325                                         ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
1326         fhDispEtaDispPhi[i]->SetXTitle("#sigma^{2}_{#eta #eta}");
1327         fhDispEtaDispPhi[i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1328         outputContainer->Add(fhDispEtaDispPhi[i]);
1329
1330         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]),
1331                                         ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
1332         fhLambda0DispEta[i]->SetXTitle("#lambda^{2}_{0}");
1333         fhLambda0DispEta[i]->SetYTitle("#sigma^{2}_{#eta #eta}");
1334         outputContainer->Add(fhLambda0DispEta[i]);
1335
1336         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]),
1337                                         ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
1338         fhLambda0DispPhi[i]->SetXTitle("#lambda^{2}_{0}");
1339         fhLambda0DispPhi[i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1340         outputContainer->Add(fhLambda0DispPhi[i]);
1341
1342       }
1343     }
1344
1345     fhNLocMaxPt = new TH2F("hNLocMaxPt","Number of local maxima in cluster, selected clusters",
1346                            nptbins,ptmin,ptmax,20,0,20);
1347     fhNLocMaxPt ->SetYTitle("N maxima");
1348     fhNLocMaxPt ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1349     outputContainer->Add(fhNLocMaxPt) ;
1350
1351     if(fFillAllNLMHistograms)
1352     {
1353       for(Int_t iSM = 0; iSM < GetCaloUtils()->GetNumberOfSuperModulesUsed(); iSM++)
1354       {
1355         fhNLocMaxPtSM[iSM] = new TH2F(Form("hNLocMaxPt_SM%d",iSM),Form("Number of local maxima in cluster, selected clusters in SM %d",iSM),
1356                                       nptbins,ptmin,ptmax,20,0,20);
1357         fhNLocMaxPtSM[iSM] ->SetYTitle("N maxima");
1358         fhNLocMaxPtSM[iSM] ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1359         outputContainer->Add(fhNLocMaxPtSM[iSM]) ;
1360       }
1361
1362       for (Int_t i = 0; i < 3; i++)
1363       {
1364         fhPtLambda0LocMax[i]  = new TH2F(Form("hPtLambda0LocMax%d",i+1),
1365                                          Form("Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #lambda_{0}, NLM=%s",nlm[i].Data()),
1366                                          nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1367         fhPtLambda0LocMax[i]->SetYTitle("#lambda_{0}^{2}");
1368         fhPtLambda0LocMax[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1369         outputContainer->Add(fhPtLambda0LocMax[i]) ;
1370
1371         if(IsDataMC())
1372         {
1373           for(Int_t ipart = 0; ipart < fgkNmcTypes; ipart++)
1374           {
1375             fhMCPtLambda0LocMax[ipart][i]  = new TH2F
1376             (Form("hPtLambda0LocMax%d_MC%s",i+1,pname[ipart].Data()),
1377              Form("Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #lambda_{0}, NLM=%s, MC %s",nlm[i].Data(),ptype[ipart].Data()),
1378              nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1379             fhMCPtLambda0LocMax[ipart][i]->SetYTitle("#lambda_{0}^{2}");
1380             fhMCPtLambda0LocMax[ipart][i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1381             outputContainer->Add(fhMCPtLambda0LocMax[ipart][i]) ;
1382           }
1383         }
1384
1385         fhPtLambda1LocMax[i]  = new TH2F(Form("hPtLambda1LocMax%d",i+1),
1386                                          Form("Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #lambda_{1}, %s",nlm[i].Data()),
1387                                          nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1388         fhPtLambda1LocMax[i]->SetYTitle("#lambda_{1}^{2}");
1389         fhPtLambda1LocMax[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1390         outputContainer->Add(fhPtLambda1LocMax[i]) ;
1391
1392         if(fCalorimeter == "EMCAL" && !fFillOnlySimpleSSHisto)
1393         {
1394           fhPtDispersionLocMax[i]  = new TH2F(Form("hPtDispersionLocMax%d",i+1),
1395                                               Form("Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs dispersion^{2}, %s",nlm[i].Data()),
1396                                               nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1397           fhPtDispersionLocMax[i]->SetYTitle("dispersion^{2}");
1398           fhPtDispersionLocMax[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1399           outputContainer->Add(fhPtDispersionLocMax[i]) ;
1400
1401           fhPtDispEtaLocMax[i]  = new TH2F(Form("hPtDispEtaLocMax%d",i+1),
1402                                            Form("Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #sigma_{#eta #eta}, %s",nlm[i].Data()),
1403                                            nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1404           fhPtDispEtaLocMax[i]->SetYTitle("#sigma_{#eta #eta}");
1405           fhPtDispEtaLocMax[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1406           outputContainer->Add(fhPtDispEtaLocMax[i]) ;
1407
1408           fhPtDispPhiLocMax[i]  = new TH2F(Form("hPtDispPhiLocMax%d",i+1),
1409                                            Form("Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #sigma_{#phi #phi}, %s",nlm[i].Data()),
1410                                            nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1411           fhPtDispPhiLocMax[i]->SetYTitle("#sigma_{#phi #phi}");
1412           fhPtDispPhiLocMax[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1413           outputContainer->Add(fhPtDispPhiLocMax[i]) ;
1414
1415           fhPtSumEtaPhiLocMax[i]  = new TH2F(Form("hPtSumEtaPhiLocMax%d",i+1),
1416                                              Form("Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #sigma_{#eta #phi}, %s",nlm[i].Data()),
1417                                              nptbins,ptmin,ptmax,2*ssbins,-ssmax,ssmax);
1418           fhPtSumEtaPhiLocMax[i]->SetYTitle("#sigma_{#eta #phi}");
1419           fhPtSumEtaPhiLocMax[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1420           outputContainer->Add(fhPtSumEtaPhiLocMax[i]) ;
1421
1422           fhPtDispEtaPhiDiffLocMax[i]  = new TH2F(Form("hPtDispEtaPhiDiffLocMax%d",i+1),
1423                                                   Form("Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #sigma_{#phi #phi} - #sigma_{#eta #eta}, %s",nlm[i].Data()),
1424                                                   nptbins,ptmin,ptmax,200, -10,10);
1425           fhPtDispEtaPhiDiffLocMax[i]->SetYTitle("#sigma_{#phi #phi} - #sigma_{#eta #eta}");
1426           fhPtDispEtaPhiDiffLocMax[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1427           outputContainer->Add(fhPtDispEtaPhiDiffLocMax[i]) ;
1428
1429           fhPtSphericityLocMax[i]  = new TH2F(Form("hPtSphericityLocMax%d",i+1),
1430                                               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()),
1431                                               nptbins,ptmin,ptmax,200, -1,1);
1432           fhPtSphericityLocMax[i]->SetYTitle("#sigma_{#phi #phi} - #sigma_{#eta #eta} / (#sigma_{#phi #phi} + #sigma_{#eta #eta})");
1433           fhPtSphericityLocMax[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1434           outputContainer->Add(fhPtSphericityLocMax[i]) ;
1435         }
1436       }
1437     } // all NLM histos
1438
1439     fhPtNCells  = new TH2F ("hPtNCells","N cells in cluster vs E ", nptbins,ptmin,ptmax, nbins,nmin,nmax);
1440     fhPtNCells->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1441     fhPtNCells->SetYTitle("# of cells in cluster");
1442     outputContainer->Add(fhPtNCells);
1443
1444     fhPtTime = new TH2F("hPtTime","cluster time vs pair E",nptbins,ptmin,ptmax, tbins,tmin,tmax);
1445     fhPtTime->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1446     fhPtTime->SetYTitle("t (ns)");
1447     outputContainer->Add(fhPtTime);
1448
1449   }
1450
1451   if(fAnaType != kIMCaloTracks)
1452   {
1453     fhEPairDiffTime = new TH2F("hEPairDiffTime","cluster pair time difference vs E",nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
1454     fhEPairDiffTime->SetXTitle("#it{E}_{pair} (GeV)");
1455     fhEPairDiffTime->SetYTitle("#Delta t (ns)");
1456     outputContainer->Add(fhEPairDiffTime);
1457   }
1458
1459   if(fAnaType == kIMCalo)
1460   {
1461     TString combiName [] = {"1LocMax","2LocMax","NLocMax","1LocMax2LocMax","1LocMaxNLocMax","2LocMaxNLocMax","1LocMaxSSBad","NLocMaxSSGood"};
1462     TString combiTitle[] = {"1 Local Maxima in both clusters","2 Local Maxima in both clusters","more than 2 Local Maxima in both clusters",
1463       "1 Local Maxima paired with 2 Local Maxima","1 Local Maxima paired with more than 2 Local Maxima",
1464       "2 Local Maxima paired with more than 2 Local Maxima",
1465       "1 Local Maxima paired with #lambda_{0}^{2}>0.3","N Local Maxima paired with 0.1<#lambda_{0}^{2}<0.3"};
1466
1467     if(fFillAllNLMHistograms)
1468     {
1469       for (Int_t i = 0; i < 8 ; i++)
1470       {
1471         if (fAnaType == kIMCaloTracks && i > 2 ) continue ;
1472
1473         fhMassPairLocMax[i]  = new TH2F
1474         (Form("MassPairLocMax%s",combiName[i].Data()),
1475          Form("#it{M} for decay #gamma pair vs #it{E}_{pair}, origin #pi^{0}, %s", combiTitle[i].Data()),
1476          nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
1477         fhMassPairLocMax[i]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
1478         fhMassPairLocMax[i]->SetXTitle("#it{E}_{pair} (GeV)");
1479         outputContainer->Add(fhMassPairLocMax[i]) ;
1480       }
1481     }
1482   }
1483
1484   if(fFillTMHisto)
1485   {
1486     fhTrackMatchedDEta  = new TH2F
1487     ("hTrackMatchedDEta",
1488      "d#eta of cluster-track vs cluster #it{p}_{T}",
1489      nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
1490     fhTrackMatchedDEta->SetYTitle("d#eta");
1491     fhTrackMatchedDEta->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1492
1493     fhTrackMatchedDPhi  = new TH2F
1494     ("hTrackMatchedDPhi",
1495      "d#phi of cluster-track vs cluster #it{p}_{T}",
1496      nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
1497     fhTrackMatchedDPhi->SetYTitle("d#phi (rad)");
1498     fhTrackMatchedDPhi->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1499
1500     fhTrackMatchedDEtaDPhi  = new TH2F
1501     ("hTrackMatchedDEtaDPhi",
1502      "d#eta vs d#phi of cluster-track",
1503      nresetabins,resetamin,resetamax,nresphibins,resphimin,resphimax);
1504     fhTrackMatchedDEtaDPhi->SetYTitle("d#phi (rad)");
1505     fhTrackMatchedDEtaDPhi->SetXTitle("d#eta");
1506
1507     outputContainer->Add(fhTrackMatchedDEta) ;
1508     outputContainer->Add(fhTrackMatchedDPhi) ;
1509     outputContainer->Add(fhTrackMatchedDEtaDPhi) ;
1510
1511     fhTrackMatchedDEtaPos  = new TH2F
1512     ("hTrackMatchedDEtaPos",
1513      "d#eta of cluster-track vs cluster #it{p}_{T}",
1514      nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
1515     fhTrackMatchedDEtaPos->SetYTitle("d#eta");
1516     fhTrackMatchedDEtaPos->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1517
1518     fhTrackMatchedDPhiPos  = new TH2F
1519     ("hTrackMatchedDPhiPos",
1520      "d#phi of cluster-track vs cluster #it{p}_{T}",
1521      nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
1522     fhTrackMatchedDPhiPos->SetYTitle("d#phi (rad)");
1523     fhTrackMatchedDPhiPos->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1524
1525     fhTrackMatchedDEtaDPhiPos  = new TH2F
1526     ("hTrackMatchedDEtaDPhiPos",
1527      "d#eta vs d#phi of cluster-track",
1528      nresetabins,resetamin,resetamax,nresphibins,resphimin,resphimax);
1529     fhTrackMatchedDEtaDPhiPos->SetYTitle("d#phi (rad)");
1530     fhTrackMatchedDEtaDPhiPos->SetXTitle("d#eta");
1531
1532     outputContainer->Add(fhTrackMatchedDEtaPos) ;
1533     outputContainer->Add(fhTrackMatchedDPhiPos) ;
1534     outputContainer->Add(fhTrackMatchedDEtaDPhiPos) ;
1535
1536     fhTrackMatchedDEtaNeg  = new TH2F
1537     ("hTrackMatchedDEtaNeg",
1538      "d#eta of cluster-track vs cluster #it{p}_{T}",
1539      nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
1540     fhTrackMatchedDEtaNeg->SetYTitle("d#eta");
1541     fhTrackMatchedDEtaNeg->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1542
1543     fhTrackMatchedDPhiNeg  = new TH2F
1544     ("hTrackMatchedDPhiNeg",
1545      "d#phi of cluster-track vs cluster #it{p}_{T}",
1546      nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
1547     fhTrackMatchedDPhiNeg->SetYTitle("d#phi (rad)");
1548     fhTrackMatchedDPhiNeg->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1549
1550     fhTrackMatchedDEtaDPhiNeg  = new TH2F
1551     ("hTrackMatchedDEtaDPhiNeg",
1552      "d#eta vs d#phi of cluster-track",
1553      nresetabins,resetamin,resetamax,nresphibins,resphimin,resphimax);
1554     fhTrackMatchedDEtaDPhiNeg->SetYTitle("d#phi (rad)");
1555     fhTrackMatchedDEtaDPhiNeg->SetXTitle("d#eta");
1556
1557     outputContainer->Add(fhTrackMatchedDEtaNeg) ;
1558     outputContainer->Add(fhTrackMatchedDPhiNeg) ;
1559     outputContainer->Add(fhTrackMatchedDEtaDPhiNeg) ;
1560
1561     fhdEdx  = new TH2F ("hdEdx","matched track <dE/dx> vs cluster #it{p}_{T}", nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
1562     fhdEdx->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1563     fhdEdx->SetYTitle("<#it{dE}/#it{dx}>");
1564     outputContainer->Add(fhdEdx);
1565
1566     fhEOverP  = new TH2F ("hEOverP","matched track E/p vs cluster #it{p}_{T}", nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
1567     fhEOverP->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1568     fhEOverP->SetYTitle("#it{E}/#it{p}");
1569     outputContainer->Add(fhEOverP);
1570
1571     if(fCalorimeter=="EMCAL" &&  GetFirstSMCoveredByTRD() >=0)
1572     {
1573       fhEOverPNoTRD  = new TH2F ("hEOverPNoTRD","matched track E/p vs cluster E, SM not behind TRD ", nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
1574       fhEOverPNoTRD->SetXTitle("#it{E} (GeV)");
1575       fhEOverPNoTRD->SetYTitle("#it{E}/#it{p}");
1576       outputContainer->Add(fhEOverPNoTRD);
1577     }
1578
1579     if(IsDataMC() && fFillTMHisto)
1580     {
1581       fhTrackMatchedMCParticlePt  = new TH2F
1582       ("hTrackMatchedMCParticlePt",
1583        "Origin of particle vs energy",
1584        nptbins,ptmin,ptmax,8,0,8);
1585       fhTrackMatchedMCParticlePt->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1586       //fhTrackMatchedMCParticlePt->SetYTitle("Particle type");
1587
1588       fhTrackMatchedMCParticlePt->GetYaxis()->SetBinLabel(1 ,"Photon");
1589       fhTrackMatchedMCParticlePt->GetYaxis()->SetBinLabel(2 ,"Electron");
1590       fhTrackMatchedMCParticlePt->GetYaxis()->SetBinLabel(3 ,"Meson Merged");
1591       fhTrackMatchedMCParticlePt->GetYaxis()->SetBinLabel(4 ,"Rest");
1592       fhTrackMatchedMCParticlePt->GetYaxis()->SetBinLabel(5 ,"Conv. Photon");
1593       fhTrackMatchedMCParticlePt->GetYaxis()->SetBinLabel(6 ,"Conv. Electron");
1594       fhTrackMatchedMCParticlePt->GetYaxis()->SetBinLabel(7 ,"Conv. Merged");
1595       fhTrackMatchedMCParticlePt->GetYaxis()->SetBinLabel(8 ,"Conv. Rest");
1596
1597       outputContainer->Add(fhTrackMatchedMCParticlePt);
1598
1599       fhTrackMatchedMCParticleDEta  = new TH2F
1600       ("hTrackMatchedMCParticleDEta",
1601        "Origin of particle vs #eta residual",
1602        nresetabins,resetamin,resetamax,8,0,8);
1603       fhTrackMatchedMCParticleDEta->SetXTitle("#Delta #eta");
1604       //fhTrackMatchedMCParticleDEta->SetYTitle("Particle type");
1605
1606       fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(1 ,"Photon");
1607       fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(2 ,"Electron");
1608       fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(3 ,"Meson Merged");
1609       fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(4 ,"Rest");
1610       fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(5 ,"Conv. Photon");
1611       fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(6 ,"Conv. Electron");
1612       fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(7 ,"Conv. Merged");
1613       fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(8 ,"Conv. Rest");
1614
1615       outputContainer->Add(fhTrackMatchedMCParticleDEta);
1616
1617       fhTrackMatchedMCParticleDPhi  = new TH2F
1618       ("hTrackMatchedMCParticleDPhi",
1619        "Origin of particle vs #phi residual",
1620        nresphibins,resphimin,resphimax,8,0,8);
1621       fhTrackMatchedMCParticleDPhi->SetXTitle("#Delta #phi");
1622       //fhTrackMatchedMCParticleDPhi->SetYTitle("Particle type");
1623
1624       fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(1 ,"Photon");
1625       fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(2 ,"Electron");
1626       fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(3 ,"Meson Merged");
1627       fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(4 ,"Rest");
1628       fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(5 ,"Conv. Photon");
1629       fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(6 ,"Conv. Electron");
1630       fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(7 ,"Conv. Merged");
1631       fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(8 ,"Conv. Rest");
1632
1633       outputContainer->Add(fhTrackMatchedMCParticleDPhi);
1634     }
1635   }
1636
1637   if(fFillWeightHistograms)
1638   {
1639     fhECellClusterRatio  = new TH2F ("hECellClusterRatio"," cell energy / cluster energy vs cluster energy, for selected decay photons from neutral meson",
1640                                      nptbins,ptmin,ptmax, 100,0,1.);
1641     fhECellClusterRatio->SetXTitle("#it{E}_{cluster} (GeV) ");
1642     fhECellClusterRatio->SetYTitle("#it{E}_{cell i}/#it{E}_{cluster}");
1643     outputContainer->Add(fhECellClusterRatio);
1644
1645     fhECellClusterLogRatio  = new TH2F ("hECellClusterLogRatio"," Log(cell energy / cluster energy) vs cluster energy, for selected decay photons from neutral meson",
1646                                         nptbins,ptmin,ptmax, 100,-10,0);
1647     fhECellClusterLogRatio->SetXTitle("#it{E}_{cluster} (GeV) ");
1648     fhECellClusterLogRatio->SetYTitle("Log (#it{E}_{max cell}/#it{E}_{cluster})");
1649     outputContainer->Add(fhECellClusterLogRatio);
1650
1651     fhEMaxCellClusterRatio  = new TH2F ("hEMaxCellClusterRatio"," max cell energy / cluster energy vs cluster energy, for selected decay photons from neutral meson",
1652                                         nptbins,ptmin,ptmax, 100,0,1.);
1653     fhEMaxCellClusterRatio->SetXTitle("#it{E}_{cluster} (GeV) ");
1654     fhEMaxCellClusterRatio->SetYTitle("#it{E}_{max cell}/#it{E}_{cluster}");
1655     outputContainer->Add(fhEMaxCellClusterRatio);
1656
1657     fhEMaxCellClusterLogRatio  = new TH2F ("hEMaxCellClusterLogRatio"," Log(max cell energy / cluster energy) vs cluster energy, for selected decay photons from neutral meson",
1658                                            nptbins,ptmin,ptmax, 100,-10,0);
1659     fhEMaxCellClusterLogRatio->SetXTitle("#it{E}_{cluster} (GeV) ");
1660     fhEMaxCellClusterLogRatio->SetYTitle("Log (#it{E}_{max cell}/#it{E}_{cluster})");
1661     outputContainer->Add(fhEMaxCellClusterLogRatio);
1662
1663     for(Int_t iw = 0; iw < 14; iw++)
1664     {
1665       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),
1666                                       nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1667       fhLambda0ForW0[iw]->SetXTitle("#it{E}_{cluster}");
1668       fhLambda0ForW0[iw]->SetYTitle("#lambda^{2}_{0}");
1669       outputContainer->Add(fhLambda0ForW0[iw]);
1670
1671       //      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),
1672       //                                      nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1673       //      fhLambda1ForW0[iw]->SetXTitle("#it{E}_{cluster}");
1674       //      fhLambda1ForW0[iw]->SetYTitle("#lambda^{2}_{1}");
1675       //      outputContainer->Add(fhLambda1ForW0[iw]);
1676
1677     }
1678   }
1679
1680   if(IsDataMC())
1681   {
1682     // Origin
1683
1684     fhMCPi0PtOrigin     = new TH2F("hMCPi0PtOrigin","Reconstructed pair from generated #pi^{0} #it{p}_{T} vs origin",nptbins,ptmin,ptmax,11,0,11) ;
1685     fhMCPi0PtOrigin->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1686     fhMCPi0PtOrigin->SetYTitle("Origin");
1687     fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(1 ,"Status 21");
1688     fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(2 ,"Quark");
1689     fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(3 ,"qq Resonances");
1690     fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(4 ,"Resonances");
1691     fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(5 ,"#rho");
1692     fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(6 ,"#omega");
1693     fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(7 ,"K");
1694     fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(8 ,"Other");
1695     fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(9 ,"#eta");
1696     fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(10 ,"#eta prime");
1697     outputContainer->Add(fhMCPi0PtOrigin) ;
1698
1699     fhMCEtaPtOrigin     = new TH2F("hMCEtaPtOrigin","Reconstructed pair from generated #pi^{0} #it{p}_{T} vs origin",nptbins,ptmin,ptmax,7,0,7) ;
1700     fhMCEtaPtOrigin->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1701     fhMCEtaPtOrigin->SetYTitle("Origin");
1702     fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(1 ,"Status 21");
1703     fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(2 ,"Quark");
1704     fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(3 ,"qq Resonances");
1705     fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(4 ,"Resonances");
1706     fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(5 ,"Other");
1707     fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(6 ,"#eta prime");
1708     outputContainer->Add(fhMCEtaPtOrigin) ;
1709
1710     fhMCPi0ProdVertex = new TH2F("hMCPi0ProdVertex","Selected reco pair from generated #pi^{0} #it{p}_{T} vs production vertex",200,ptmin,20+ptmin,5000,0,500) ;
1711     fhMCPi0ProdVertex->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1712     fhMCPi0ProdVertex->SetYTitle("#it{R} (cm)");
1713     outputContainer->Add(fhMCPi0ProdVertex) ;
1714
1715     fhMCEtaProdVertex = new TH2F("hMCEtaProdVertex","Selected reco pair from generated #eta #it{p}_{T} vs production vertex",200,ptmin,20+ptmin,5000,0,500) ;
1716     fhMCEtaProdVertex->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1717     fhMCEtaProdVertex->SetYTitle("#it{R} (cm)");
1718     outputContainer->Add(fhMCEtaProdVertex) ;
1719
1720     if(GetReader()->GetDataType() != AliCaloTrackReader::kMC && fAnaType==kSSCalo)
1721     {
1722       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",
1723                                           nptbins,ptmin,ptmax,200,0,2);
1724       fhMCPi0PtGenRecoFraction->SetXTitle("#it{p}^{rec}_{T} (GeV/#it{c})");
1725       fhMCPi0PtGenRecoFraction->SetYTitle("#it{E}^{#pi^{0} mother} / #it{E}^{rec}");
1726       outputContainer->Add(fhMCPi0PtGenRecoFraction) ;
1727
1728       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",
1729                                           nptbins,ptmin,ptmax,200,0,2);
1730       fhMCEtaPtGenRecoFraction->SetXTitle("#it{p}^{rec}_{T} (GeV/#it{c})");
1731       fhMCEtaPtGenRecoFraction->SetYTitle("#it{E}^{ #eta mother} / #it{E}^{rec}");
1732       outputContainer->Add(fhMCEtaPtGenRecoFraction) ;
1733
1734       fhMCPi0DecayPt = new TH1F("hMCPi0DecayPt","Number of #gamma from #pi^{0} decay  identified as #pi^{0} (#eta)",nptbins,ptmin,ptmax);
1735       fhMCPi0DecayPt->SetYTitle("#it{N}");
1736       fhMCPi0DecayPt->SetXTitle("#it{p}^{rec}_{T} (GeV/#it{c})");
1737       outputContainer->Add(fhMCPi0DecayPt) ;
1738
1739       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}",
1740                                         nptbins,ptmin,ptmax,100,0,1);
1741       fhMCPi0DecayPtFraction->SetXTitle("p^{rec}_{T} (GeV/#it{c})");
1742       fhMCPi0DecayPtFraction->SetYTitle("E^{gen} / E^{gen-mother}");
1743       outputContainer->Add(fhMCPi0DecayPtFraction) ;
1744
1745       fhMCEtaDecayPt = new TH1F("hMCEtaDecayPt","Number of #gamma from #eta decay  identified as #pi^{0} (#eta)",nptbins,ptmin,ptmax);
1746       fhMCEtaDecayPt->SetYTitle("#it{N}");
1747       fhMCEtaDecayPt->SetXTitle("#it{p}^{rec}_{T} (GeV/#it{c})");
1748       outputContainer->Add(fhMCEtaDecayPt) ;
1749
1750       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",
1751                                         nptbins,ptmin,ptmax,100,0,1);
1752       fhMCEtaDecayPtFraction->SetXTitle("#it{p}^{rec}_{T} (GeV/#it{c})");
1753       fhMCEtaDecayPtFraction->SetYTitle("#it{E}^{gen} / #it{E}^{gen-mother}");
1754       outputContainer->Add(fhMCEtaDecayPtFraction) ;
1755
1756       fhMCOtherDecayPt = new TH1F("hMCOtherDecayPt","Number of #gamma decay (not #eta or #pi^{0})  identified as #pi^{0} (#eta)",nptbins,ptmin,ptmax);
1757       fhMCOtherDecayPt->SetYTitle("#it{N}");
1758       fhMCOtherDecayPt->SetXTitle("#it{p}^{rec}_{T} (GeV/#it{c})");
1759       outputContainer->Add(fhMCOtherDecayPt) ;
1760     }
1761
1762     if(fAnaType!=kSSCalo)
1763     {
1764       fhAnglePairMCPi0  = new TH2F
1765       ("AnglePairMCPi0",
1766        "Angle between decay #gamma pair vs #it{E}_{pair}, origin #pi^{0}",nptbins,ptmin,ptmax,250,0,0.5);
1767       fhAnglePairMCPi0->SetYTitle("#alpha (rad)");
1768       fhAnglePairMCPi0->SetXTitle("#it{E}_{pair} (GeV)");
1769       outputContainer->Add(fhAnglePairMCPi0) ;
1770
1771       fhAnglePairMCEta  = new TH2F
1772       ("AnglePairMCEta",
1773        "Angle between decay #gamma pair vs #it{E}_{pair}, origin #eta",nptbins,ptmin,ptmax,250,0,0.5);
1774       fhAnglePairMCEta->SetYTitle("#alpha (rad)");
1775       fhAnglePairMCEta->SetXTitle("#it{E}_{pair} (GeV)");
1776       outputContainer->Add(fhAnglePairMCEta) ;
1777
1778       fhMassPairMCPi0  = new TH2F
1779       ("MassPairMCPi0",
1780        "#it{M} for decay #gamma pair vs #it{E}_{pair}, origin #pi^{0}",nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
1781       fhMassPairMCPi0->SetYTitle("#it{M} (GeV/#it{c}^{2})");
1782       fhMassPairMCPi0->SetXTitle("#it{E}_{pair} (GeV)");
1783       outputContainer->Add(fhMassPairMCPi0) ;
1784
1785       fhMassPairMCEta  = new TH2F
1786       ("MassPairMCEta",
1787        "#it{M} for decay #gamma pair vs #it{E}_{pair}, origin #eta",nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
1788       fhMassPairMCEta->SetYTitle("#it{M} (GeV/#it{c}^{2})");
1789       fhMassPairMCEta->SetXTitle("#it{E}_{pair} (GeV)");
1790       outputContainer->Add(fhMassPairMCEta) ;
1791     }
1792
1793     for(Int_t i = 0; i < fgkNmcTypes; i++)
1794     {
1795       fhMCE[i]  = new TH1F
1796       (Form("hE_MC%s",pname[i].Data()),
1797        Form("Identified as #pi^{0} (#eta), cluster from %s",
1798             ptype[i].Data()),
1799        nptbins,ptmin,ptmax);
1800       fhMCE[i]->SetYTitle("#it{N}");
1801       fhMCE[i]->SetXTitle("#it{E} (GeV)");
1802       outputContainer->Add(fhMCE[i]) ;
1803
1804       fhMCPt[i]  = new TH1F
1805       (Form("hPt_MC%s",pname[i].Data()),
1806        Form("Identified as #pi^{0} (#eta), cluster from %s",
1807             ptype[i].Data()),
1808        nptbins,ptmin,ptmax);
1809       fhMCPt[i]->SetYTitle("#it{N}");
1810       fhMCPt[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1811       outputContainer->Add(fhMCPt[i]) ;
1812
1813       if(fFillHighMultHistograms)
1814       {
1815         fhMCPtCentrality[i]  = new TH2F
1816         (Form("hPtCentrality_MC%s",pname[i].Data()),
1817          Form("Identified as #pi^{0} (#eta), cluster from %s",
1818               ptype[i].Data()),
1819          nptbins,ptmin,ptmax, 100,0,100);
1820         fhMCPtCentrality[i]->SetYTitle("centrality");
1821         fhMCPtCentrality[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1822         outputContainer->Add(fhMCPtCentrality[i]) ;
1823       }
1824
1825       if(fAnaType == kSSCalo)
1826       {
1827         fhMCNLocMaxPt[i] = new TH2F
1828         (Form("hNLocMaxPt_MC%s",pname[i].Data()),
1829          Form("cluster from %s, #it{p}_{T} of cluster vs NLM, accepted",ptype[i].Data()),
1830          nptbins,ptmin,ptmax,20,0,20);
1831         fhMCNLocMaxPt[i] ->SetYTitle("#it{NLM}");
1832         fhMCNLocMaxPt[i] ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1833         outputContainer->Add(fhMCNLocMaxPt[i]) ;
1834
1835         fhMCNLocMaxPtReject[i] = new TH2F
1836         (Form("hNLocMaxPtReject_MC%s",pname[i].Data()),
1837          Form("cluster from %s, #it{p}_{T} of cluster vs NLM, rejected",ptype[i].Data()),
1838          nptbins,ptmin,ptmax,20,0,20);
1839         fhMCNLocMaxPtReject[i] ->SetYTitle("#it{NLM}");
1840         fhMCNLocMaxPtReject[i] ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1841         outputContainer->Add(fhMCNLocMaxPtReject[i]) ;
1842
1843         fhMCEReject[i]  = new TH1F
1844         (Form("hEReject_MC%s",pname[i].Data()),
1845          Form("Rejected as #pi^{0} (#eta), cluster from %s",
1846               ptype[i].Data()),
1847          nptbins,ptmin,ptmax);
1848         fhMCEReject[i]->SetYTitle("#it{N}");
1849         fhMCEReject[i]->SetXTitle("#it{E} (GeV)");
1850         outputContainer->Add(fhMCEReject[i]) ;
1851
1852         fhMCPtReject[i]  = new TH1F
1853         (Form("hPtReject_MC%s",pname[i].Data()),
1854          Form("Rejected as #pi^{0} (#eta), cluster from %s",
1855               ptype[i].Data()),
1856          nptbins,ptmin,ptmax);
1857         fhMCPtReject[i]->SetYTitle("#it{N}");
1858         fhMCPtReject[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1859         outputContainer->Add(fhMCPtReject[i]) ;
1860       }
1861
1862       fhMCPtPhi[i]  = new TH2F
1863       (Form("hPtPhi_MC%s",pname[i].Data()),
1864        Form("Identified as #pi^{0} (#eta), cluster from %s",ptype[i].Data()),
1865        nptbins,ptmin,ptmax,nphibins,phimin,phimax);
1866       fhMCPtPhi[i]->SetYTitle("#phi");
1867       fhMCPtPhi[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1868       outputContainer->Add(fhMCPtPhi[i]) ;
1869
1870       fhMCPtEta[i]  = new TH2F
1871       (Form("hPtEta_MC%s",pname[i].Data()),
1872        Form("Identified as #pi^{0} (#eta), cluster from %s",
1873             ptype[i].Data()),nptbins,ptmin,ptmax,netabins,etamin,etamax);
1874       fhMCPtEta[i]->SetYTitle("#eta");
1875       fhMCPtEta[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1876       outputContainer->Add(fhMCPtEta[i]) ;
1877
1878       fhMCMassPt[i]  = new TH2F
1879       (Form("hMassPt_MC%s",pname[i].Data()),
1880        Form("all pairs #it{M}: #it{p}_{T} vs #it{M} from %s",ptype[i].Data()),
1881        nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1882       fhMCMassPt[i]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
1883       fhMCMassPt[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1884       outputContainer->Add(fhMCMassPt[i]) ;
1885
1886       fhMCSelectedMassPt[i]  = new TH2F
1887       (Form("hSelectedMassPt_MC%s",pname[i].Data()),
1888        Form("Selected #pi^{0} (#eta) pairs #it{M}: #it{p}_{T} vs #it{M} from %s",ptype[i].Data()),
1889        nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1890       fhMCSelectedMassPt[i]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
1891       fhMCSelectedMassPt[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1892       outputContainer->Add(fhMCSelectedMassPt[i]) ;
1893
1894       if(fAnaType == kSSCalo)
1895       {
1896         fhMCMassPtNoOverlap[i]  = new TH2F
1897         (Form("hMassPtNoOverlap_MC%s",pname[i].Data()),
1898          Form("all pairs #it{M}: #it{p}_{T} vs #it{M} from %s, no overlap",ptype[i].Data()),
1899          nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1900         fhMCMassPt[i]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
1901         fhMCMassPt[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1902         outputContainer->Add(fhMCMassPtNoOverlap[i]) ;
1903
1904         fhMCSelectedMassPtNoOverlap[i]  = new TH2F
1905         (Form("hSelectedMassPtNoOverlap_MC%s",pname[i].Data()),
1906          Form("Selected #pi^{0} (#eta) pairs #it{M}: #it{p}_{T} vs #it{M} from %s, no overlap",ptype[i].Data()),
1907          nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1908         fhMCSelectedMassPtNoOverlap[i]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
1909         fhMCSelectedMassPtNoOverlap[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1910         outputContainer->Add(fhMCSelectedMassPtNoOverlap[i]) ;
1911       }
1912
1913       if( fFillSelectClHisto )
1914       {
1915         fhMCPtLambda0[i]  = new TH2F(Form("hELambda0_MC%s",pname[i].Data()),
1916                                      Form("Selected pair, cluster from %s : #it{p}_{T} vs #lambda_{0}^{2}",ptype[i].Data()),
1917                                      nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1918         fhMCPtLambda0[i]->SetYTitle("#lambda_{0}^{2}");
1919         fhMCPtLambda0[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1920         outputContainer->Add(fhMCPtLambda0[i]) ;
1921
1922         fhMCPtLambda1[i]  = new TH2F(Form("hELambda1_MC%s",pname[i].Data()),
1923                                      Form("Selected pair, cluster from %s : #it{p}_{T} vs #lambda_{1}^{2}",ptype[i].Data()),
1924                                      nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1925         fhMCPtLambda1[i]->SetYTitle("#lambda_{1}^{2}");
1926         fhMCPtLambda1[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1927         outputContainer->Add(fhMCPtLambda1[i]) ;
1928
1929         if(fCalorimeter=="EMCAL" &&  GetFirstSMCoveredByTRD() >= 0)
1930         {
1931           fhMCPtLambda0NoTRD[i]  = new TH2F(Form("hELambda0NoTRD_MC%s",pname[i].Data()),
1932                                             Form("Selected pair, cluster from %s : #it{p}_{T} vs #lambda_{0}^{2}, NoTRD",ptype[i].Data()),
1933                                             nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1934           fhMCPtLambda0NoTRD[i]->SetYTitle("#lambda_{0}^{2}");
1935           fhMCPtLambda0NoTRD[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1936           outputContainer->Add(fhMCPtLambda0NoTRD[i]) ;
1937         }
1938
1939         if(!fFillOnlySimpleSSHisto)
1940         {
1941           fhMCPtDispersion[i]  = new TH2F(Form("hEDispersion_MC%s",pname[i].Data()),
1942                                           Form("Selected pair, cluster from %s : #it{p}_{T} vs dispersion^{2}",ptype[i].Data()),
1943                                           nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1944           fhMCPtDispersion[i]->SetYTitle("#it{D}^{2}");
1945           fhMCPtDispersion[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1946           outputContainer->Add(fhMCPtDispersion[i]) ;
1947
1948           fhMCPtDispEta[i]  = new TH2F (Form("hPtDispEta_MC%s",pname[i].Data()),
1949                                         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()),
1950                                         nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
1951           fhMCPtDispEta[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1952           fhMCPtDispEta[i]->SetYTitle("#sigma^{2}_{#eta #eta}");
1953           outputContainer->Add(fhMCPtDispEta[i]);
1954
1955           fhMCPtDispPhi[i]  = new TH2F (Form("hPtDispPhi_MC%s",pname[i].Data()),
1956                                         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()),
1957                                         nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
1958           fhMCPtDispPhi[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1959           fhMCPtDispPhi[i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1960           outputContainer->Add(fhMCPtDispPhi[i]);
1961
1962           fhMCPtSumEtaPhi[i]  = new TH2F (Form("hPtSumEtaPhi_MC%s",pname[i].Data()),
1963                                           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()),
1964                                           nptbins,ptmin,ptmax, 2*ssbins,-ssmax,ssmax);
1965           fhMCPtSumEtaPhi[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1966           fhMCPtSumEtaPhi[i]->SetYTitle("#delta^{2}_{#eta #phi}");
1967           outputContainer->Add(fhMCPtSumEtaPhi[i]);
1968
1969           fhMCPtDispEtaPhiDiff[i]  = new TH2F (Form("hPtDispEtaPhiDiff_MC%s",pname[i].Data()),
1970                                                Form("cluster from %s : #sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta} vs #it{p}_{T}",ptype[i].Data()),
1971                                                nptbins,ptmin,ptmax,200,-10,10);
1972           fhMCPtDispEtaPhiDiff[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1973           fhMCPtDispEtaPhiDiff[i]->SetYTitle("#sigma^{2}_{#phi #phi}-#sigma^{2}_{#eta #eta}");
1974           outputContainer->Add(fhMCPtDispEtaPhiDiff[i]);
1975
1976           fhMCPtSphericity[i]  = new TH2F (Form("hPtSphericity_MC%s",pname[i].Data()),
1977                                            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()),
1978                                            nptbins,ptmin,ptmax, 200,-1,1);
1979           fhMCPtSphericity[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1980           fhMCPtSphericity[i]->SetYTitle("#it{s} = (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi})");
1981           outputContainer->Add(fhMCPtSphericity[i]);
1982
1983           for(Int_t ie = 0; ie < 7; ie++)
1984           {
1985             fhMCDispEtaDispPhi[ie][i] = new TH2F (Form("hMCDispEtaDispPhi_EBin%d_MC%s",ie,pname[i].Data()),
1986                                                   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]),
1987                                                   ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
1988             fhMCDispEtaDispPhi[ie][i]->SetXTitle("#sigma^{2}_{#eta #eta}");
1989             fhMCDispEtaDispPhi[ie][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1990             outputContainer->Add(fhMCDispEtaDispPhi[ie][i]);
1991
1992             fhMCLambda0DispEta[ie][i] = new TH2F (Form("hMCLambda0DispEta_EBin%d_MC%s",ie,pname[i].Data()),
1993                                                   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]),
1994                                                   ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
1995             fhMCLambda0DispEta[ie][i]->SetXTitle("#lambda^{2}_{0}");
1996             fhMCLambda0DispEta[ie][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1997             outputContainer->Add(fhMCLambda0DispEta[ie][i]);
1998
1999             fhMCLambda0DispPhi[ie][i] = new TH2F (Form("hMCLambda0DispPhi_EBin%d_MC%s",ie,pname[i].Data()),
2000                                                   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]),
2001                                                   ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
2002             fhMCLambda0DispPhi[ie][i]->SetXTitle("#lambda^{2}_{0}");
2003             fhMCLambda0DispPhi[ie][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
2004             outputContainer->Add(fhMCLambda0DispPhi[ie][i]);
2005
2006           }
2007
2008           fhMCPtLambda0FracMaxCellCut[i]  = new TH2F(Form("hELambda0FracMaxCellCut_MC%s",pname[i].Data()),
2009                                                      Form("Selected pair, cluster from %s : #it{p}_{T} vs #lambda_{0}^{2}, Max cell fraction of energy < 0.5 ",ptype[i].Data()),
2010                                                      nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2011           fhMCPtLambda0FracMaxCellCut[i]->SetYTitle("#lambda_{0}^{2}");
2012           fhMCPtLambda0FracMaxCellCut[i]->SetXTitle("#it{E} (GeV)");
2013           outputContainer->Add(fhMCPtLambda0FracMaxCellCut[i]) ;
2014
2015           fhMCPtFracMaxCell[i]  = new TH2F(Form("hEFracMaxCell_MC%s",pname[i].Data()),
2016                                            Form("Selected pair, cluster from %s : #it{p}_{T} vs Max cell fraction of energy",ptype[i].Data()),
2017                                            nptbins,ptmin,ptmax,100,0,1);
2018           fhMCPtFracMaxCell[i]->SetYTitle("#it{Fraction}");
2019           fhMCPtFracMaxCell[i]->SetXTitle("#it{E} (GeV)");
2020           outputContainer->Add(fhMCPtFracMaxCell[i]) ;
2021
2022         }
2023       }
2024     }// MC particle loop
2025   }//Histos with MC
2026
2027   if(fAnaType==kSSCalo)
2028   {
2029     fhAsymmetry  = new TH2F ("hAsymmetry","#it{A} = ( #it{E}_{1} - #it{E}_{2} ) / ( #it{E}_{1} + #it{E}_{2} ) vs #it{E}",
2030                              nptbins,ptmin,ptmax, 200, -1,1);
2031     fhAsymmetry->SetXTitle("#it{E} (GeV)");
2032     fhAsymmetry->SetYTitle("#it{A} = ( #it{E}_{1} - #it{E}_{2} ) / ( #it{E}_{1} + #it{E}_{2} )");
2033     outputContainer->Add(fhAsymmetry);
2034
2035     fhSelectedAsymmetry  = new TH2F ("hSelectedAsymmetry","#it{A} = ( #it{E}_{1} - #it{E}_{2} ) / ( #it{E}_{1} + #it{E}_{2} ) vs #it{E}",
2036                                      nptbins,ptmin,ptmax, 200, -1,1);
2037     fhSelectedAsymmetry->SetXTitle("#it{E} (GeV)");
2038     fhSelectedAsymmetry->SetYTitle("#it{A} = ( #it{E}_{1} - #it{E}_{2} ) / ( #it{E}_{1} + #it{E}_{2} )");
2039     outputContainer->Add(fhSelectedAsymmetry);
2040
2041     fhSplitE  = new TH1F
2042     ("hSplitE","Selected #pi^{0} (#eta) pairs energy sum of split sub-clusters",nptbins,ptmin,ptmax);
2043     fhSplitE->SetYTitle("counts");
2044     fhSplitE->SetXTitle("#it{E} (GeV)");
2045     outputContainer->Add(fhSplitE) ;
2046
2047     fhSplitPt  = new TH1F
2048     ("hSplitPt","Selected #pi^{0} (#eta) pairs #it{p}_{T} sum of split sub-clusters",nptbins,ptmin,ptmax);
2049     fhSplitPt->SetYTitle("counts");
2050     fhSplitPt->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2051     outputContainer->Add(fhSplitPt) ;
2052
2053
2054     fhSplitPtPhi  = new TH2F
2055     ("hSplitPtPhi","Selected #pi^{0} (#eta) pairs: sum split sub-cluster #it{p}_{T} vs #phi",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
2056     fhSplitPtPhi->SetYTitle("#phi (rad)");
2057     fhSplitPtPhi->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2058     outputContainer->Add(fhSplitPtPhi) ;
2059
2060     fhSplitPtEta  = new TH2F
2061     ("hSplitPtEta","Selected #pi^{0} (#eta) pairs: sum split sub-cluster #it{p}_{T} vs #eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
2062     fhSplitPtEta->SetYTitle("#eta");
2063     fhSplitPtEta->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2064     outputContainer->Add(fhSplitPtEta) ;
2065
2066
2067     fhNLocMaxSplitPt = new TH2F("hNLocMaxSplitPt","Number of local maxima in cluster",
2068                                 nptbins,ptmin,ptmax,20,0,20);
2069     fhNLocMaxSplitPt ->SetYTitle("#it{NLM}");
2070     fhNLocMaxSplitPt ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2071     outputContainer->Add(fhNLocMaxSplitPt) ;
2072
2073
2074     fhMassSplitPt  = new TH2F
2075     ("hMassSplitPt","all pairs #it{M}: sum split sub-cluster #it{p}_{T} vs #it{M}",
2076      nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
2077     fhMassSplitPt->SetYTitle("#it{M} (GeV/#it{c}^{2})");
2078     fhMassSplitPt->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2079     outputContainer->Add(fhMassSplitPt) ;
2080
2081     fhSelectedMassSplitPt  = new TH2F
2082     ("hSelectedMassSplitPt","Selected #pi^{0} (#eta) pairs #it{M}: sum split sub-cluster #it{p}_{T} vs #it{M}",
2083      nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
2084     fhSelectedMassSplitPt->SetYTitle("#it{M} (GeV/#it{c}^{2})");
2085     fhSelectedMassSplitPt->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2086     outputContainer->Add(fhSelectedMassSplitPt) ;
2087
2088     if(IsDataMC())
2089     {
2090       fhMassSplitPtNoOverlap  = new TH2F
2091       ("hMassSplitPtNoOverlap","all pairs #it{M}: sum split sub-cluster #it{p}_{T} vs #it{M}, no overlap",
2092        nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
2093       fhMassSplitPtNoOverlap->SetYTitle("#it{M} (GeV/#it{c}^{2})");
2094       fhMassSplitPtNoOverlap->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2095       outputContainer->Add(fhMassSplitPtNoOverlap) ;
2096
2097       fhSelectedMassSplitPtNoOverlap  = new TH2F
2098       ("hSelectedMassSplitPtNoOverlap","Selected #pi^{0} (#eta) pairs #it{M}: sum split sub-cluster #it{p}_{T} vs #it{M}, no overlap",
2099        nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
2100       fhSelectedMassSplitPtNoOverlap->SetYTitle("#it{M} (GeV/#it{c}^{2})");
2101       fhSelectedMassSplitPtNoOverlap->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2102       outputContainer->Add(fhSelectedMassSplitPtNoOverlap) ;
2103
2104
2105       fhMCPi0PtRecoPtPrim  = new TH2F
2106       ("hMCPi0PtRecoPtPrim","#it{p}_{T,reco} vs #it{p}_{T,gen}",
2107        nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2108       fhMCPi0PtRecoPtPrim ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2109       fhMCPi0PtRecoPtPrim ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2110       outputContainer->Add(fhMCPi0PtRecoPtPrim ) ;
2111
2112       fhMCPi0PtRecoPtPrimNoOverlap  = new TH2F
2113       ("hMCPi0PtRecoPtPrimNoOverlap","#it{p}_{T,reco} vs #it{p}_{T,gen}, no overlap",
2114        nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2115       fhMCPi0PtRecoPtPrimNoOverlap ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2116       fhMCPi0PtRecoPtPrimNoOverlap ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2117       outputContainer->Add(fhMCPi0PtRecoPtPrimNoOverlap ) ;
2118
2119       fhMCPi0SelectedPtRecoPtPrim  = new TH2F
2120       ("hMCPi0SelectedPtRecoPtPrim","#it{p}_{T,reco} vs #it{p}_{T,gen}",
2121        nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2122       fhMCPi0SelectedPtRecoPtPrim ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2123       fhMCPi0SelectedPtRecoPtPrim ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2124       outputContainer->Add(fhMCPi0SelectedPtRecoPtPrim ) ;
2125
2126       fhMCPi0SelectedPtRecoPtPrimNoOverlap  = new TH2F
2127       ("hMCPi0SelectedPtRecoPtPrimNoOverlap","#it{p}_{T,reco} vs #it{p}_{T,gen}, no overlap",
2128        nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2129       fhMCPi0SelectedPtRecoPtPrimNoOverlap ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2130       fhMCPi0SelectedPtRecoPtPrimNoOverlap ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2131       outputContainer->Add(fhMCPi0SelectedPtRecoPtPrimNoOverlap ) ;
2132
2133
2134       fhMCPi0SplitPtRecoPtPrim  = new TH2F
2135       ("hMCPi0SplitPtRecoPtPrim","#it{p}_{T,reco} (split sum) vs #it{p}_{T,gen}",
2136        nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2137       fhMCPi0SplitPtRecoPtPrim ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2138       fhMCPi0SplitPtRecoPtPrim ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2139       outputContainer->Add(fhMCPi0SplitPtRecoPtPrim ) ;
2140
2141       fhMCPi0SplitPtRecoPtPrimNoOverlap  = new TH2F
2142       ("hMCPi0SplitPtRecoPtPrimNoOverlap","#it{p}_{T,reco} (split sum) vs #it{p}_{T,gen}, no overlap",
2143        nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2144       fhMCPi0SplitPtRecoPtPrimNoOverlap ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2145       fhMCPi0SplitPtRecoPtPrimNoOverlap ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2146       outputContainer->Add(fhMCPi0SplitPtRecoPtPrimNoOverlap ) ;
2147
2148       fhMCPi0SelectedSplitPtRecoPtPrim  = new TH2F
2149       ("hMCPi0SelectedSplitPtRecoPtPrim","#it{p}_{T,reco} (split sum) vs #it{p}_{T,gen}",
2150        nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2151       fhMCPi0SelectedSplitPtRecoPtPrim ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2152       fhMCPi0SelectedSplitPtRecoPtPrim ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2153       outputContainer->Add(fhMCPi0SelectedSplitPtRecoPtPrim ) ;
2154
2155       fhMCPi0SelectedSplitPtRecoPtPrimNoOverlap  = new TH2F
2156       ("hMCPi0SelectedSplitPtRecoPtPrimNoOverlap","#it{p}_{T,reco} (split sum) vs #it{p}_{T,gen}, no overlap",
2157        nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2158       fhMCPi0SelectedSplitPtRecoPtPrimNoOverlap ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2159       fhMCPi0SelectedSplitPtRecoPtPrimNoOverlap ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2160       outputContainer->Add(fhMCPi0SelectedSplitPtRecoPtPrimNoOverlap ) ;
2161
2162       fhMCEtaPtRecoPtPrim  = new TH2F
2163       ("hMCEtaPtRecoPtPrim","#it{p}_{T,reco} vs #it{p}_{T,gen}",
2164        nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2165       fhMCEtaPtRecoPtPrim ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2166       fhMCEtaPtRecoPtPrim ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2167       outputContainer->Add(fhMCEtaPtRecoPtPrim ) ;
2168
2169       fhMCEtaPtRecoPtPrimNoOverlap  = new TH2F
2170       ("hMCEtaPtRecoPtPrimNoOverlap","#it{p}_{T,reco} vs #it{p}_{T,gen}, no overlap",
2171        nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2172       fhMCEtaPtRecoPtPrimNoOverlap ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2173       fhMCEtaPtRecoPtPrimNoOverlap ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2174       outputContainer->Add(fhMCEtaPtRecoPtPrimNoOverlap ) ;
2175
2176       fhMCEtaSelectedPtRecoPtPrim  = new TH2F
2177       ("hMCEtaSelectedPtRecoPtPrim","#it{p}_{T,reco} vs #it{p}_{T,gen}",
2178        nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2179       fhMCEtaSelectedPtRecoPtPrim ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2180       fhMCEtaSelectedPtRecoPtPrim ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2181       outputContainer->Add(fhMCEtaSelectedPtRecoPtPrim ) ;
2182
2183       fhMCEtaSelectedPtRecoPtPrimNoOverlap  = new TH2F
2184       ("hMCEtaSelectedPtRecoPtPrimNoOverlap","#it{p}_{T,reco} vs #it{p}_{T,gen}, no overlap",
2185        nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2186       fhMCEtaSelectedPtRecoPtPrimNoOverlap ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2187       fhMCEtaSelectedPtRecoPtPrimNoOverlap ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2188       outputContainer->Add(fhMCEtaSelectedPtRecoPtPrimNoOverlap ) ;
2189
2190
2191       fhMCEtaSplitPtRecoPtPrim  = new TH2F
2192       ("hMCEtaSplitPtRecoPtPrim","#it{p}_{T,reco} (split sum) vs #it{p}_{T,gen}",
2193        nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2194       fhMCEtaSplitPtRecoPtPrim ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2195       fhMCEtaSplitPtRecoPtPrim ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2196       outputContainer->Add(fhMCEtaSplitPtRecoPtPrim ) ;
2197
2198       fhMCEtaSplitPtRecoPtPrimNoOverlap  = new TH2F
2199       ("hMCEtaSplitPtRecoPtPrimNoOverlap","#it{p}_{T,reco} (split sum) vs #it{p}_{T,gen}, no overlap",
2200        nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2201       fhMCEtaSplitPtRecoPtPrimNoOverlap ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2202       fhMCEtaSplitPtRecoPtPrimNoOverlap ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2203       outputContainer->Add(fhMCEtaSplitPtRecoPtPrimNoOverlap ) ;
2204
2205       fhMCEtaSelectedSplitPtRecoPtPrim  = new TH2F
2206       ("hMCEtaSelectedSplitPtRecoPtPrim","#it{p}_{T,reco} (split sum) vs #it{p}_{T,gen}",
2207        nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2208       fhMCEtaSelectedSplitPtRecoPtPrim ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2209       fhMCEtaSelectedSplitPtRecoPtPrim ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2210       outputContainer->Add(fhMCEtaSelectedSplitPtRecoPtPrim ) ;
2211
2212       fhMCEtaSelectedSplitPtRecoPtPrimNoOverlap  = new TH2F
2213       ("hMCEtaSelectedSplitPtRecoPtPrimNoOverlap","#it{p}_{T,reco} (split sum) vs #it{p}_{T,gen}, no overlap",
2214        nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2215       fhMCEtaSelectedSplitPtRecoPtPrimNoOverlap ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2216       fhMCEtaSelectedSplitPtRecoPtPrimNoOverlap ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2217       outputContainer->Add(fhMCEtaSelectedSplitPtRecoPtPrimNoOverlap ) ;
2218
2219
2220       for(Int_t inlm = 0; inlm < 3; inlm++)
2221       {
2222         fhMCPi0PtRecoPtPrimLocMax[inlm]  = new TH2F
2223         (Form("hMCPi0PtRecoPtPrimLocMax%d",inlm+1),Form("#it{p}_{T,reco} vs #it{p}_{T,gen}, %s",nlm[inlm].Data()),
2224          nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2225         fhMCPi0PtRecoPtPrimLocMax[inlm] ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2226         fhMCPi0PtRecoPtPrimLocMax[inlm] ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2227         outputContainer->Add(fhMCPi0PtRecoPtPrimLocMax[inlm] ) ;
2228
2229         fhMCPi0SelectedPtRecoPtPrimLocMax[inlm]  = new TH2F
2230         (Form("hMCPi0SelectedPtRecoPtPrimLocMax%d",inlm+1),Form("#it{p}_{T,reco} vs #it{p}_{T,gen}, %s",nlm[inlm].Data()),
2231          nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2232         fhMCPi0SelectedPtRecoPtPrimLocMax[inlm] ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2233         fhMCPi0SelectedPtRecoPtPrimLocMax[inlm] ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2234         outputContainer->Add(fhMCPi0SelectedPtRecoPtPrimLocMax[inlm] ) ;
2235
2236         fhMCPi0SplitPtRecoPtPrimLocMax[inlm]  = new TH2F
2237         (Form("hMCPi0SplitPtRecoPtPrimLocMax%d",inlm+1),Form("#it{p}_{T,reco} (split sum) vs #it{p}_{T,gen}, %s",nlm[inlm].Data()),
2238          nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2239         fhMCPi0SplitPtRecoPtPrimLocMax[inlm] ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2240         fhMCPi0SplitPtRecoPtPrimLocMax[inlm] ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2241         outputContainer->Add(fhMCPi0SplitPtRecoPtPrimLocMax[inlm] ) ;
2242
2243         fhMCPi0SelectedSplitPtRecoPtPrimLocMax[inlm]  = new TH2F
2244         (Form("hMCPi0SelectedSplitPtRecoPtPrimLocMax%d",inlm+1),Form("#it{p}_{T,reco} (split sum) vs #it{p}_{T,gen}, %s",nlm[inlm].Data()),
2245          nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2246         fhMCPi0SelectedSplitPtRecoPtPrimLocMax[inlm] ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2247         fhMCPi0SelectedSplitPtRecoPtPrimLocMax[inlm] ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2248         outputContainer->Add(fhMCPi0SelectedSplitPtRecoPtPrimLocMax[inlm] ) ;
2249
2250         fhMCEtaPtRecoPtPrimLocMax[inlm]  = new TH2F
2251         (Form("hMCEtaPtRecoPtPrimLocMax%d",inlm+1),Form("#it{p}_{T,reco} vs #it{p}_{T,gen}, %s",nlm[inlm].Data()),
2252          nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2253         fhMCEtaPtRecoPtPrimLocMax[inlm] ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2254         fhMCEtaPtRecoPtPrimLocMax[inlm] ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2255         outputContainer->Add(fhMCEtaPtRecoPtPrimLocMax[inlm] ) ;
2256
2257         fhMCEtaSelectedPtRecoPtPrimLocMax[inlm]  = new TH2F
2258         (Form("hMCEtaSelectedPtRecoPtPrimLocMax%d",inlm+1),Form("#it{p}_{T,reco} vs #it{p}_{T,gen}, %s",nlm[inlm].Data()),
2259          nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2260         fhMCEtaSelectedPtRecoPtPrimLocMax[inlm] ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2261         fhMCEtaSelectedPtRecoPtPrimLocMax[inlm] ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2262         outputContainer->Add(fhMCEtaSelectedPtRecoPtPrimLocMax[inlm] ) ;
2263
2264         fhMCEtaSplitPtRecoPtPrimLocMax[inlm]  = new TH2F
2265         (Form("hMCEtaSplitPtRecoPtPrimLocMax%d",inlm+1),Form("#it{p}_{T,reco} (split sum) vs #it{p}_{T,gen}, %s",nlm[inlm].Data()),
2266          nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2267         fhMCEtaSplitPtRecoPtPrimLocMax[inlm] ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2268         fhMCEtaSplitPtRecoPtPrimLocMax[inlm] ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2269         outputContainer->Add(fhMCEtaSplitPtRecoPtPrimLocMax[inlm] ) ;
2270
2271         fhMCEtaSelectedSplitPtRecoPtPrimLocMax[inlm]  = new TH2F
2272         (Form("hMCEtaSelectedSplitPtRecoPtPrimLocMax%d",inlm+1),Form("#it{p}_{T,reco} (split sum) vs #it{p}_{T,gen}, %s",nlm[inlm].Data()),
2273          nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2274         fhMCEtaSelectedSplitPtRecoPtPrimLocMax[inlm] ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
2275         fhMCEtaSelectedSplitPtRecoPtPrimLocMax[inlm] ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
2276         outputContainer->Add(fhMCEtaSelectedSplitPtRecoPtPrimLocMax[inlm] ) ;
2277
2278       }
2279
2280       for(Int_t i = 0; i < fgkNmcTypes; i++)
2281       {
2282         fhMCPtAsymmetry[i]  = new TH2F (Form("hEAsymmetry_MC%s",pname[i].Data()),
2283                                         Form("cluster from %s : #it{A} = ( #it{E}_{1} - #it{E}_{2} ) / ( #it{E}_{1} + #it{E}_{2} ) vs #it{E}",ptype[i].Data()),
2284                                         nptbins,ptmin,ptmax, 200,-1,1);
2285         fhMCPtAsymmetry[i]->SetXTitle("#it{E} (GeV)");
2286         fhMCPtAsymmetry[i]->SetYTitle("#it{A} = ( #it{E}_{1} - #it{E}_{2} ) / ( #it{E}_{1} + #it{E}_{2} )");
2287         outputContainer->Add(fhMCPtAsymmetry[i]);
2288
2289         fhMCSplitE[i]  = new TH1F
2290         (Form("hSplitE_MC%s",pname[i].Data()),
2291          Form("cluster from %s, energy sum of split sub-clusters",ptype[i].Data()),
2292          nptbins,ptmin,ptmax);
2293         fhMCSplitE[i]->SetYTitle("counts");
2294         fhMCSplitE[i]->SetXTitle("#it{E} (GeV)");
2295         outputContainer->Add(fhMCSplitE[i]) ;
2296
2297         fhMCSplitPt[i]  = new TH1F
2298         (Form("hSplitPt_MC%s",pname[i].Data()),
2299          Form("cluster from %s, #it{p}_{T} sum of split sub-clusters",ptype[i].Data()),
2300          nptbins,ptmin,ptmax);
2301         fhMCSplitPt[i]->SetYTitle("counts");
2302         fhMCSplitPt[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2303         outputContainer->Add(fhMCSplitPt[i]) ;
2304
2305
2306         fhMCSplitPtPhi[i]  = new TH2F
2307         (Form("hSplitPtPhi_MC%s",pname[i].Data()),
2308          Form("Identified as #pi^{0} (#eta), cluster from %s",ptype[i].Data()),
2309          nptbins,ptmin,ptmax,nphibins,phimin,phimax);
2310         fhMCSplitPtPhi[i]->SetYTitle("#phi");
2311         fhMCSplitPtPhi[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2312         outputContainer->Add(fhMCSplitPtPhi[i]) ;
2313
2314         fhMCSplitPtEta[i]  = new TH2F
2315         (Form("hSplitPtEta_MC%s",pname[i].Data()),
2316          Form("Identified as #pi^{0} (#eta), cluster from %s",
2317               ptype[i].Data()),nptbins,ptmin,ptmax,netabins,etamin,etamax);
2318         fhMCSplitPtEta[i]->SetYTitle("#eta");
2319         fhMCSplitPtEta[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2320         outputContainer->Add(fhMCSplitPtEta[i]) ;
2321
2322
2323         fhMCNLocMaxSplitPt[i] = new TH2F
2324         (Form("hNLocMaxSplitPt_MC%s",pname[i].Data()),
2325          Form("cluster from %s, #it{p}_{T} sum of split sub-clusters, for NLM",ptype[i].Data()),
2326          nptbins,ptmin,ptmax,20,0,20);
2327         fhMCNLocMaxSplitPt[i] ->SetYTitle("#it{NLM}");
2328         fhMCNLocMaxSplitPt[i] ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2329         outputContainer->Add(fhMCNLocMaxSplitPt[i]) ;
2330
2331         fhMCMassSplitPt[i]  = new TH2F
2332         (Form("hMassSplitPt_MC%s",pname[i].Data()),
2333          Form("all pairs #it{M}: split #it{p}_{T} vs #it{M} from %s",ptype[i].Data()),
2334          nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
2335         fhMCMassSplitPt[i]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
2336         fhMCMassSplitPt[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2337         outputContainer->Add(fhMCMassSplitPt[i]) ;
2338
2339         fhMCSelectedMassSplitPt[i]  = new TH2F
2340         (Form("hSelectedMassSplitPt_MC%s",pname[i].Data()),
2341          Form("Selected #pi^{0} (#eta) pairs #it{M}: split #it{p}_{T} vs #it{M} from %s",ptype[i].Data()),
2342          nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
2343         fhMCSelectedMassSplitPt[i]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
2344         fhMCSelectedMassSplitPt[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2345         outputContainer->Add(fhMCSelectedMassSplitPt[i]) ;
2346
2347         fhMCMassSplitPtNoOverlap[i]  = new TH2F
2348         (Form("hMassSplitPtNoOverlap_MC%s",pname[i].Data()),
2349          Form("all pairs #it{M}: split #it{p}_{T} vs #it{M} from %s, no overlap",ptype[i].Data()),
2350          nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
2351         fhMCMassSplitPtNoOverlap[i]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
2352         fhMCMassSplitPtNoOverlap[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2353         outputContainer->Add(fhMCMassSplitPtNoOverlap[i]) ;
2354
2355         fhMCSelectedMassSplitPtNoOverlap[i]  = new TH2F
2356         (Form("hSelectedMassSplitPtNoOverlap_MC%s",pname[i].Data()),
2357          Form("Selected #pi^{0} (#eta) pairs #it{M}: split #it{p}_{T} vs #it{M} from %s, no overlap",ptype[i].Data()),
2358          nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
2359         fhMCSelectedMassSplitPtNoOverlap[i]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
2360         fhMCSelectedMassSplitPtNoOverlap[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2361         outputContainer->Add(fhMCSelectedMassSplitPtNoOverlap[i]) ;
2362       }
2363     }
2364   }
2365
2366   if(fAnaType==kSSCalo && fFillSelectClHisto && !fFillOnlySimpleSSHisto )
2367   {
2368     for(Int_t i = 0; i< 3; i++)
2369     {
2370       fhPtAsymmetryLocMax[i]  = new TH2F(Form("hEAsymmetryLocMax%d",i+1),
2371                                          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()),
2372                                          nptbins,ptmin,ptmax,200, -1,1);
2373       fhPtAsymmetryLocMax[i]->SetYTitle("#it{A} = ( #it{E}_{1} - #it{E}_{2} ) / ( #it{E}_{1} + #it{E}_{2} )");
2374       fhPtAsymmetryLocMax[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2375       outputContainer->Add(fhPtAsymmetryLocMax[i]) ;
2376     }
2377
2378     for(Int_t ie = 0; ie < 7; ie++)
2379     {
2380
2381       fhAsymmetryLambda0[ie] = new TH2F (Form("hAsymmetryLambda0_EBin%d",ie),
2382                                          Form("#lambda_{0}^{2} vs A for %d < #it{E} < %d GeV",bin[ie],bin[ie+1]),
2383                                          ssbins,ssmin,ssmax , 200,-1,1);
2384       fhAsymmetryLambda0[ie]->SetXTitle("#lambda_{0}^{2}");
2385       fhAsymmetryLambda0[ie]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
2386       outputContainer->Add(fhAsymmetryLambda0[ie]);
2387
2388       fhAsymmetryDispEta[ie] = new TH2F (Form("hAsymmetryDispEta_EBin%d",ie),
2389                                          Form("#sigma^{2}_{#eta #eta} vs #it{A} for %d < #it{E} < %d GeV",bin[ie],bin[ie+1]),
2390                                          ssbins,ssmin,ssmax , 200,-1,1);
2391       fhAsymmetryDispEta[ie]->SetXTitle("#sigma^{2}_{#eta #eta}");
2392       fhAsymmetryDispEta[ie]->SetYTitle("#it{A} = ( #it{E}_{1} - #it{E}_{2} ) / ( #it{E}_{1} + #it{E}_{2} )");
2393       outputContainer->Add(fhAsymmetryDispEta[ie]);
2394
2395       fhAsymmetryDispPhi[ie] = new TH2F (Form("hAsymmetryDispPhi_EBin%d",ie),
2396                                          Form("#sigma^{2}_{#phi #phi} vs #it{A} for %d < #it{E} < %d GeV",bin[ie],bin[ie+1]),
2397                                          ssbins,ssmin,ssmax , 200,-1,1);
2398       fhAsymmetryDispPhi[ie]->SetXTitle("#sigma^{2}_{#phi #phi}");
2399       fhAsymmetryDispPhi[ie]->SetYTitle("#it{A} = ( #it{E}_{1} - #it{E}_{2} ) / ( #it{E}_{1} + #it{E}_{2} )");
2400       outputContainer->Add(fhAsymmetryDispPhi[ie]);
2401     }
2402
2403
2404     if(IsDataMC())
2405     {
2406       for(Int_t i = 0; i < fgkNmcTypes; i++)
2407       {
2408         for(Int_t ie = 0; ie < 7; ie++)
2409         {
2410           fhMCAsymmetryLambda0[ie][i] = new TH2F (Form("hMCAsymmetryLambda0_EBin%d_MC%s",ie,pname[i].Data()),
2411                                                   Form("cluster from %s : #lambda_{0}^{2} vs A for %d < #it{E} < %d GeV",ptype[i].Data(),bin[ie],bin[ie+1]),
2412                                                   ssbins,ssmin,ssmax , 200,-1,1);
2413           fhMCAsymmetryLambda0[ie][i]->SetXTitle("#lambda_{0}^{2}");
2414           fhMCAsymmetryLambda0[ie][i]->SetYTitle("#it{A} = ( #it{E}_{1} - #it{E}_{2} ) / ( #it{E}_{1} + #it{E}_{2} )");
2415           outputContainer->Add(fhMCAsymmetryLambda0[ie][i]);
2416
2417           fhMCAsymmetryDispEta[ie][i] = new TH2F (Form("hMCAsymmetryDispEta_EBin%d_MC%s",ie,pname[i].Data()),
2418                                                   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]),
2419                                                   ssbins,ssmin,ssmax , 200,-1,1);
2420           fhMCAsymmetryDispEta[ie][i]->SetXTitle("#sigma^{2}_{#eta #eta}");
2421           fhMCAsymmetryDispEta[ie][i]->SetYTitle("#it{A} = ( #it{E}_{1} - #it{E}_{2} ) / ( #it{E}_{1} + #it{E}_{2} )");
2422           outputContainer->Add(fhMCAsymmetryDispEta[ie][i]);
2423
2424           fhMCAsymmetryDispPhi[ie][i] = new TH2F (Form("hMCAsymmetryDispPhi_EBin%d_MC%s",ie,pname[i].Data()),
2425                                                   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]),
2426                                                   ssbins,ssmin,ssmax , 200,-1,1);
2427           fhMCAsymmetryDispPhi[ie][i]->SetXTitle("#sigma^{2}_{#phi #phi}");
2428           fhMCAsymmetryDispPhi[ie][i]->SetYTitle("#it{A} = ( #it{E}_{1} - #it{E}_{2} ) / ( #it{E}_{1} + #it{E}_{2} )");
2429           outputContainer->Add(fhMCAsymmetryDispPhi[ie][i]);
2430         }
2431       }
2432     }
2433   }
2434
2435   if(fFillPileUpHistograms)
2436   {
2437
2438     TString pileUpName[] = {"SPD","EMCAL","SPDOrEMCAL","SPDAndEMCAL","SPDAndNotEMCAL","EMCALAndNotSPD","NotSPDAndNotEMCAL"} ;
2439
2440     for(Int_t i = 0 ; i < 7 ; i++)
2441     {
2442       fhPtPileUp[i]  = new TH1F(Form("hPtPileUp%s",pileUpName[i].Data()),
2443                                 Form("Selected #pi^{0} (#eta) #it{p}_{T} distribution, %s Pile-Up event",pileUpName[i].Data()), nptbins,ptmin,ptmax);
2444       fhPtPileUp[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2445       outputContainer->Add(fhPtPileUp[i]);
2446
2447       fhPtCellTimePileUp[i]  = new TH2F(Form("hPtCellTimePileUp%s",pileUpName[i].Data()),
2448                                         Form("Pt vs cell time in cluster, %s Pile-Up event",pileUpName[i].Data()),
2449                                         nptbins,ptmin,ptmax,ntimptbins,timemin,timemax);
2450       fhPtCellTimePileUp[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2451       fhPtCellTimePileUp[i]->SetYTitle("#it{t}_{cell} (ns)");
2452       outputContainer->Add(fhPtCellTimePileUp[i]);
2453
2454       fhPtTimeDiffPileUp[i]  = new TH2F(Form("hPtTimeDiffPileUp%s",pileUpName[i].Data()),
2455                                         Form("Pt vs t_{max}-t_{cell} in cluster, %s Pile-Up event",pileUpName[i].Data()),
2456                                         nptbins,ptmin,ptmax,400,-200,200);
2457       fhPtTimeDiffPileUp[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2458       fhPtTimeDiffPileUp[i]->SetYTitle("#it{t}_{max}-#it{t}_{cell} (ns)");
2459       outputContainer->Add(fhPtTimeDiffPileUp[i]);
2460
2461     }
2462
2463     fhTimePtNoCut  = new TH2F ("hTimePt_NoCut","#it{t} of cluster vs #it{E} of clusters, no cut", nptbins,ptmin,ptmax, ntimptbins,timemin,timemax);
2464     fhTimePtNoCut->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2465     fhTimePtNoCut->SetYTitle("#it{t} (ns)");
2466     outputContainer->Add(fhTimePtNoCut);
2467
2468     fhTimePtSPD  = new TH2F ("hTimePt_SPD","#it{t} of cluster vs #it{E} of clusters, SPD cut", nptbins,ptmin,ptmax, ntimptbins,timemin,timemax);
2469     fhTimePtSPD->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2470     fhTimePtSPD->SetYTitle("#it{t} (ns)");
2471     outputContainer->Add(fhTimePtSPD);
2472
2473     fhTimePtSPDMulti  = new TH2F ("hTimePt_SPDMulti","time of cluster vs #it{E} of clusters, SPD multi cut", nptbins,ptmin,ptmax, ntimptbins,timemin,timemax);
2474     fhTimePtSPDMulti->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2475     fhTimePtSPDMulti->SetYTitle("#it{t} (ns)");
2476     outputContainer->Add(fhTimePtSPDMulti);
2477
2478     fhTimeNPileUpVertSPD  = new TH2F ("hTime_NPileUpVertSPD","#it{t} of cluster vs #it{N} pile-up SPD vertex", ntimptbins,timemin,timemax,50,0,50);
2479     fhTimeNPileUpVertSPD->SetYTitle("# vertex ");
2480     fhTimeNPileUpVertSPD->SetXTitle("#it{t} (ns)");
2481     outputContainer->Add(fhTimeNPileUpVertSPD);
2482
2483     fhTimeNPileUpVertTrack  = new TH2F ("hTime_NPileUpVertTracks","#it{t} of cluster vs #it{N} pile-up Tracks vertex", ntimptbins,timemin,timemax, 50,0,50 );
2484     fhTimeNPileUpVertTrack->SetYTitle("# vertex ");
2485     fhTimeNPileUpVertTrack->SetXTitle("#it{t} (ns)");
2486     outputContainer->Add(fhTimeNPileUpVertTrack);
2487
2488     fhTimeNPileUpVertContributors  = new TH2F ("hTime_NPileUpVertContributors","#it{t} of cluster vs #it{N} constributors to pile-up SPD vertex", ntimptbins,timemin,timemax,50,0,50);
2489     fhTimeNPileUpVertContributors->SetYTitle("# vertex ");
2490     fhTimeNPileUpVertContributors->SetXTitle("#it{t} (ns)");
2491     outputContainer->Add(fhTimeNPileUpVertContributors);
2492
2493     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);
2494     fhTimePileUpMainVertexZDistance->SetYTitle("distance #it{Z} (cm) ");
2495     fhTimePileUpMainVertexZDistance->SetXTitle("#it{t} (ns)");
2496     outputContainer->Add(fhTimePileUpMainVertexZDistance);
2497
2498     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);
2499     fhTimePileUpMainVertexZDiamond->SetYTitle("diamond distance #it{Z} (cm) ");
2500     fhTimePileUpMainVertexZDiamond->SetXTitle("#it{t} (ns)");
2501     outputContainer->Add(fhTimePileUpMainVertexZDiamond);
2502
2503                 fhPtNPileUpSPDVtx  = new TH2F ("hPt_NPileUpVertSPD","#it{p}_{T} of cluster vs #it{N} pile-up SPD vertex",
2504                                                                                                                                          nptbins,ptmin,ptmax,20,0,20);
2505                 fhPtNPileUpSPDVtx->SetYTitle("# vertex ");
2506                 fhPtNPileUpSPDVtx->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2507                 outputContainer->Add(fhPtNPileUpSPDVtx);
2508
2509                 fhPtNPileUpTrkVtx  = new TH2F ("hPt_NPileUpVertTracks","#it{p}_{T} of cluster vs #it{N} pile-up Tracks vertex",
2510                                                                                                                                          nptbins,ptmin,ptmax, 20,0,20 );
2511                 fhPtNPileUpTrkVtx->SetYTitle("# vertex ");
2512                 fhPtNPileUpTrkVtx->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2513                 outputContainer->Add(fhPtNPileUpTrkVtx);
2514
2515                 fhPtNPileUpSPDVtxTimeCut  = new TH2F ("hPt_NPileUpVertSPD_TimeCut","#it{p}_{T} of cluster vs N pile-up SPD vertex, |tof| < 25 ns",
2516                                           nptbins,ptmin,ptmax,20,0,20);
2517                 fhPtNPileUpSPDVtxTimeCut->SetYTitle("# vertex ");
2518                 fhPtNPileUpSPDVtxTimeCut->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2519                 outputContainer->Add(fhPtNPileUpSPDVtxTimeCut);
2520
2521                 fhPtNPileUpTrkVtxTimeCut  = new TH2F ("hPt_NPileUpVertTracks_TimeCut","#it{p}_{T} of cluster vs N pile-up Tracks vertex, |tof| < 25 ns",
2522                                                                                                                                                                         nptbins,ptmin,ptmax, 20,0,20 );
2523                 fhPtNPileUpTrkVtxTimeCut->SetYTitle("# vertex ");
2524                 fhPtNPileUpTrkVtxTimeCut->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2525                 outputContainer->Add(fhPtNPileUpTrkVtxTimeCut);
2526
2527     fhPtNPileUpSPDVtxTimeCut2  = new TH2F ("hPt_NPileUpVertSPD_TimeCut2","#it{p}_{T} of cluster vs N pile-up SPD vertex, -25 < tof < 75 ns",
2528                                            nptbins,ptmin,ptmax,20,0,20);
2529                 fhPtNPileUpSPDVtxTimeCut2->SetYTitle("# vertex ");
2530                 fhPtNPileUpSPDVtxTimeCut2->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2531                 outputContainer->Add(fhPtNPileUpSPDVtxTimeCut2);
2532
2533                 fhPtNPileUpTrkVtxTimeCut2  = new TH2F ("hPt_NPileUpVertTracks_TimeCut2","#it{p}_{T} of cluster vs N pile-up Tracks vertex, -25 < tof < 75 ns",
2534                                            nptbins,ptmin,ptmax, 20,0,20 );
2535                 fhPtNPileUpTrkVtxTimeCut2->SetYTitle("# vertex ");
2536                 fhPtNPileUpTrkVtxTimeCut2->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2537                 outputContainer->Add(fhPtNPileUpTrkVtxTimeCut2);
2538
2539   }
2540
2541   //Keep neutral meson selection histograms if requiered
2542   //Setting done in AliNeutralMesonSelection
2543
2544   if(fAnaType!=kSSCalo && GetNeutralMesonSelection())
2545   {
2546     TList * nmsHistos = GetNeutralMesonSelection()->GetCreateOutputObjects() ;
2547
2548     if(GetNeutralMesonSelection()->AreNeutralMesonSelectionHistosKept())
2549       for(Int_t i = 0; i < nmsHistos->GetEntries(); i++) outputContainer->Add(nmsHistos->At(i)) ;
2550
2551     delete nmsHistos;
2552   }
2553
2554   return outputContainer ;
2555
2556 }
2557
2558 //_____________________________________________
2559 Int_t AliAnaPi0EbE::GetMCIndex(const Int_t tag)
2560 {
2561   // Assign mc index depending on MC bit set
2562
2563   if       ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0)  )
2564   {
2565     return kmcPi0 ;
2566   }//pi0
2567   else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta)  )
2568   {
2569     return kmcEta ;
2570   }//eta
2571   else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt) ||
2572              GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation) ||
2573              GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCISR))
2574   {
2575     return kmcPhoton ;
2576   }//direct photon
2577   else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) &&
2578              GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay) )
2579   {
2580     return kmcPi0Decay ;
2581   }//decay photon from pi0
2582   else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) &&
2583              GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay) )
2584   {
2585     return kmcEtaDecay ;
2586   }//decay photon from eta
2587   else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) &&
2588              GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay) )
2589   {
2590     return kmcOtherDecay ;
2591   }//decay photon from other than eta or pi0
2592   else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron))
2593   {
2594     return kmcElectron ;
2595   }//electron
2596   else
2597   {
2598     return kmcHadron ;
2599   }//other particles
2600
2601 }
2602
2603 //__________________________________________________________________
2604 void AliAnaPi0EbE::HasPairSameMCMother(AliAODPWG4Particle * photon1,
2605                                        AliAODPWG4Particle * photon2,
2606                                        Int_t & label, Int_t & tag)
2607 {
2608   // Check the labels of pare in case mother was same pi0 or eta
2609   // Set the new AOD accordingly
2610
2611   Int_t  label1 = photon1->GetLabel();
2612   Int_t  label2 = photon2->GetLabel();
2613
2614   if(label1 < 0 || label2 < 0 ) return ;
2615
2616   //Int_t tag1 = GetMCAnalysisUtils()->CheckOrigin(label1, GetReader());
2617   //Int_t tag2 = GetMCAnalysisUtils()->CheckOrigin(label2, GetReader());
2618   Int_t tag1 = photon1->GetTag();
2619   Int_t tag2 = photon2->GetTag();
2620
2621   if(GetDebug() > 0) printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - Origin of: photon1 %d; photon2 %d \n",tag1, tag2);
2622   if( (GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCPi0Decay) &&
2623        GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCPi0Decay)    ) ||
2624       (GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCEtaDecay) &&
2625        GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCEtaDecay)    )
2626      )
2627   {
2628
2629     //Check if pi0/eta mother is the same
2630     if(GetReader()->ReadStack())
2631     {
2632       if(label1>=0)
2633       {
2634         TParticle * mother1 = GetMCStack()->Particle(label1);//photon in kine tree
2635         label1 = mother1->GetFirstMother();
2636         //mother1 = GetMCStack()->Particle(label1);//pi0
2637       }
2638       if(label2>=0)
2639       {
2640         TParticle * mother2 = GetMCStack()->Particle(label2);//photon in kine tree
2641         label2 = mother2->GetFirstMother();
2642         //mother2 = GetMCStack()->Particle(label2);//pi0
2643       }
2644     } // STACK
2645     else if(GetReader()->ReadAODMCParticles())
2646     {//&& (input > -1)){
2647       if(label1>=0)
2648       {
2649         AliAODMCParticle * mother1 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles())->At(label1);//photon in kine tree
2650         label1 = mother1->GetMother();
2651         //mother1 = GetMCStack()->Particle(label1);//pi0
2652       }
2653       if(label2>=0)
2654       {
2655         AliAODMCParticle * mother2 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles())->At(label2);//photon in kine tree
2656         label2 = mother2->GetMother();
2657         //mother2 = GetMCStack()->Particle(label2);//pi0
2658       }
2659     }// AOD
2660
2661     //printf("mother1 %d, mother2 %d\n",label1,label2);
2662     if( label1 == label2 && label1>=0 )
2663     {
2664       label = label1;
2665
2666       TLorentzVector mom1 = *(photon1->Momentum());
2667       TLorentzVector mom2 = *(photon2->Momentum());
2668
2669       Double_t angle = mom2.Angle(mom1.Vect());
2670       Double_t mass  = (mom1+mom2).M();
2671       Double_t epair = (mom1+mom2).E();
2672
2673       if(GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCPi0Decay))
2674       {
2675         fhMassPairMCPi0 ->Fill(epair,mass);
2676         fhAnglePairMCPi0->Fill(epair,angle);
2677         GetMCAnalysisUtils()->SetTagBit(tag,AliMCAnalysisUtils::kMCPi0);
2678       }
2679       else
2680       {
2681         fhMassPairMCEta ->Fill(epair,mass);
2682         fhAnglePairMCEta->Fill(epair,angle);
2683         GetMCAnalysisUtils()->SetTagBit(tag,AliMCAnalysisUtils::kMCEta);
2684       }
2685
2686     } // same label
2687   } // both from eta or pi0 decay
2688
2689 }
2690
2691 //____________________________________________________________________________
2692 void AliAnaPi0EbE::Init()
2693 {
2694   //Init
2695   //Do some checks
2696   if(fCalorimeter == "PHOS" && !GetReader()->IsPHOSSwitchedOn() && NewOutputAOD()){
2697     printf("AliAnaPi0EbE::Init() - !!STOP: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!\n");
2698     abort();
2699   }
2700   else  if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn() && NewOutputAOD()){
2701     printf("AliAnaPi0EbE::Init() - !!STOP: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!\n");
2702     abort();
2703   }
2704
2705 }
2706
2707 //____________________________________________________________________________
2708 void AliAnaPi0EbE::InitParameters()
2709 {
2710   //Initialize the parameters of the analysis.
2711   AddToHistogramsName("AnaPi0EbE_");
2712
2713   fInputAODGammaConvName = "PhotonsCTS" ;
2714   fAnaType = kIMCalo ;
2715   fCalorimeter = "EMCAL" ;
2716   fMinDist  = 2.;
2717   fMinDist2 = 4.;
2718   fMinDist3 = 5.;
2719
2720   fNLMECutMin[0] = 10.;
2721   fNLMECutMin[1] = 6. ;
2722   fNLMECutMin[2] = 6. ;
2723 }
2724
2725 //__________________________________________________________________
2726 void  AliAnaPi0EbE::MakeAnalysisFillAOD()
2727 {
2728   //Do analysis and fill aods
2729
2730   switch(fAnaType)
2731   {
2732     case kIMCalo:
2733       MakeInvMassInCalorimeter();
2734       break;
2735
2736     case kSSCalo:
2737       MakeShowerShapeIdentification();
2738       break;
2739
2740     case kIMCaloTracks:
2741       MakeInvMassInCalorimeterAndCTS();
2742       break;
2743
2744   }
2745 }
2746
2747 //____________________________________________
2748 void  AliAnaPi0EbE::MakeInvMassInCalorimeter()
2749 {
2750   //Do analysis and fill aods
2751   //Search for the photon decay in calorimeters
2752   //Read photon list from AOD, produced in class AliAnaPhoton
2753   //Check if 2 photons have the mass of the pi0.
2754
2755   TLorentzVector mom1;
2756   TLorentzVector mom2;
2757   TLorentzVector mom ;
2758
2759   Int_t tag   = 0;
2760   Int_t label = 0;
2761
2762   if(!GetInputAODBranch())
2763   {
2764     AliFatal(Form("No input calo photons in AOD with name branch < %s >, STOP \n",GetInputAODName().Data()));
2765     return; // coverity
2766   }
2767
2768   //Get shower shape information of clusters
2769   TObjArray *clusters = 0;
2770   if     (fCalorimeter=="EMCAL") clusters = GetEMCALClusters();
2771   else if(fCalorimeter=="PHOS")  clusters = GetPHOSClusters() ;
2772
2773   for(Int_t iphoton = 0; iphoton < GetInputAODBranch()->GetEntriesFast()-1; iphoton++)
2774   {
2775     AliAODPWG4Particle * photon1 =  (AliAODPWG4Particle*) (GetInputAODBranch()->At(iphoton));
2776
2777     // Vertex cut in case of mixed events
2778     Int_t evtIndex1 = 0 ;
2779     if(GetMixedEvent())
2780     {
2781       evtIndex1 = GetMixedEvent()->EventIndexForCaloCluster(photon1->GetCaloLabel(0)) ;
2782       if(TMath::Abs(GetVertex(evtIndex1)[2]) > GetZvertexCut()) continue ;  //vertex cut
2783     }
2784
2785     mom1 = *(photon1->Momentum());
2786
2787     //Get original cluster, to recover some information
2788     Int_t iclus = -1;
2789     AliVCluster *cluster1 = FindCluster(clusters,photon1->GetCaloLabel(0),iclus);
2790
2791     if(!cluster1)
2792     {
2793       printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - First cluster not found\n");
2794       return;
2795     }
2796
2797     for(Int_t jphoton = iphoton+1; jphoton < GetInputAODBranch()->GetEntriesFast(); jphoton++)
2798     {
2799       AliAODPWG4Particle * photon2 =  (AliAODPWG4Particle*) (GetInputAODBranch()->At(jphoton));
2800
2801       // Do analysis only when one of the decays is isolated
2802       // Run AliAnaParticleIsolation before
2803       if(fSelectIsolatedDecay)
2804       {
2805         Bool_t isolated1 = ((AliAODPWG4ParticleCorrelation*) photon1)->IsIsolated();
2806         Bool_t isolated2 = ((AliAODPWG4ParticleCorrelation*) photon2)->IsIsolated();
2807         if(!isolated1 && !isolated2) continue;
2808       }
2809
2810       // Vertex cut in case of mixed events
2811       Int_t evtIndex2 = 0 ;
2812       if(GetMixedEvent())
2813       {
2814         evtIndex2 = GetMixedEvent()->EventIndexForCaloCluster(photon2->GetCaloLabel(0)) ;
2815
2816         if(evtIndex1 == evtIndex2)
2817           continue ;
2818
2819         if(TMath::Abs(GetVertex(evtIndex2)[2]) > GetZvertexCut()) continue ;  //vertex cut
2820       }
2821
2822       mom2 = *(photon2->Momentum());
2823
2824       //Get original cluster, to recover some information
2825       Int_t iclus2 = -1;
2826       AliVCluster *cluster2 = FindCluster(clusters,photon2->GetCaloLabel(0),iclus2,iclus+1);
2827       // start new loop from iclus1+1 to gain some time
2828
2829       if(!cluster2)
2830       {
2831         printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - Second cluster not found\n");
2832         return;
2833       }
2834
2835       Float_t e1    = photon1->E();
2836       Float_t e2    = photon2->E();
2837
2838       //Select clusters with good time window difference
2839       Float_t tof1  = cluster1->GetTOF()*1e9;;
2840       Float_t tof2  = cluster2->GetTOF()*1e9;;
2841       Double_t t12diff = tof1-tof2;
2842       fhEPairDiffTime->Fill(e1+e2,    t12diff);
2843       if(TMath::Abs(t12diff) > GetPairTimeCut()) continue;
2844
2845       //Play with the MC stack if available
2846       Int_t mcIndex = kmcHadron;
2847       if(IsDataMC())
2848       {
2849         HasPairSameMCMother(photon1, photon2, label, tag) ;
2850         mcIndex = GetMCIndex(tag);
2851       }
2852
2853       // Check the invariant mass for different selection on the local maxima
2854       // Name of AOD method TO BE FIXED
2855       Int_t nMaxima1 = photon1->GetFiducialArea();
2856       Int_t nMaxima2 = photon2->GetFiducialArea();
2857
2858       mom = mom1+mom2;
2859
2860       Double_t mass  = mom.M();
2861       Double_t epair = mom.E();
2862       Float_t ptpair = mom.Pt();
2863
2864       if(fFillAllNLMHistograms)
2865       {
2866         if(nMaxima1==nMaxima2)
2867         {
2868           if     (nMaxima1==1) fhMassPairLocMax[0]->Fill(epair,mass);
2869           else if(nMaxima1==2) fhMassPairLocMax[1]->Fill(epair,mass);
2870           else                 fhMassPairLocMax[2]->Fill(epair,mass);
2871         }
2872         else if(nMaxima1==1 || nMaxima2==1)
2873         {
2874           if  (nMaxima1==2 || nMaxima2==2) fhMassPairLocMax[3]->Fill(epair,mass);
2875           else                             fhMassPairLocMax[4]->Fill(epair,mass);
2876         }
2877         else
2878           fhMassPairLocMax[5]->Fill(epair,mass);
2879
2880         // combinations with SS axis cut and NLM cut
2881         if(nMaxima1 == 1 && cluster2->GetM02() > 0.3) fhMassPairLocMax[6]->Fill(epair,mass);
2882         if(nMaxima2 == 1 && cluster1->GetM02() > 0.3) fhMassPairLocMax[6]->Fill(epair,mass);
2883         if(nMaxima1 >  1 && cluster2->GetM02() < 0.3 && cluster2->GetM02()> 0.1 ) fhMassPairLocMax[7]->Fill(epair,mass);
2884         if(nMaxima2 >  1 && cluster1->GetM02() < 0.3 && cluster1->GetM02()> 0.1 ) fhMassPairLocMax[7]->Fill(epair,mass);
2885       }
2886
2887       //
2888       // Skip events with too few or too many  NLM
2889       //
2890       if((nMaxima1 < fNLMCutMin || nMaxima1 > fNLMCutMax) || (nMaxima2 < fNLMCutMin || nMaxima2 > fNLMCutMax)) continue ;
2891
2892       if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - NLM of out of range: cluster1 %d, cluster2 %d \n",nMaxima1, nMaxima2);
2893
2894       //Mass of all pairs
2895       fhMass  ->Fill( epair,mass);
2896       fhMassPt->Fill(ptpair,mass);
2897       if(IsDataMC()) fhMCMassPt[mcIndex]->Fill(ptpair,mass);
2898
2899       //
2900       // Select good pair (good phi, pt cuts, aperture and invariant mass)
2901       //
2902       if(!GetNeutralMesonSelection()->SelectPair(mom1, mom2,fCalorimeter)) continue;
2903
2904       if(GetDebug()>1)
2905         printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - Selected gamma pair: pt %f, phi %f, eta%f \n",
2906                mom.Pt(), mom.Phi()*TMath::RadToDeg(), mom.Eta());
2907
2908       //
2909       // Tag both photons as decay if not done before
2910       // set the corresponding bit for pi0 or eta or "side" case
2911       //
2912       Int_t bit1 = photon1->GetBtag(); // temporary
2913       if( bit1 < 0 ) bit1 = 0 ; // temporary
2914       if( !GetNeutralMesonSelection()->CheckDecayBit(bit1) )
2915       {
2916         if( GetDebug() > 1 )
2917           printf("AliAnaPi0EbE::MakeInvMassInCalorimeter - pT1 %2.2f; bit requested %d; decay bit1: In %d, ",
2918                  mom1.Pt(), GetNeutralMesonSelection()->GetDecayBit(), bit1);
2919
2920         photon1->SetTagged(kTRUE); // temporary
2921         GetNeutralMesonSelection()->SetDecayBit(bit1);
2922         photon1->SetBtag(bit1); // temporary
2923
2924         if( GetDebug() > 1 )
2925           printf("Out %d \n", bit1);
2926
2927         fhPtDecay->Fill(photon1->Pt());
2928
2929         //Fill some histograms about shower shape
2930         if(fFillSelectClHisto && cluster1 && GetReader()->GetDataType()!=AliCaloTrackReader::kMC)
2931           FillSelectedClusterHistograms(cluster1, mom1.Pt(), nMaxima1, photon1->GetTag());
2932
2933         if(IsDataMC())
2934         {
2935           Int_t mcIndex1 = GetMCIndex(photon1->GetTag());
2936           fhMCPtDecay[mcIndex1]->Fill(photon1->Pt());
2937         }
2938       }
2939
2940       Int_t bit2 = photon2->GetBtag(); // temporary
2941       if( bit2 < 0 ) bit2 = 0 ; // temporary
2942       if( !GetNeutralMesonSelection()->CheckDecayBit(bit2) )
2943       {
2944         if( GetDebug() > 1 )
2945           printf("AliAnaPi0EbE::MakeInvMassInCalorimeter - pT2 %2.2f; bit requested %d; decay bit2: In %d, ",
2946                  mom2.Pt(), GetNeutralMesonSelection()->GetDecayBit(), bit2);
2947
2948         photon2->SetTagged(kTRUE); // temporary
2949         GetNeutralMesonSelection()->SetDecayBit(bit2);
2950         photon2->SetBtag(bit2); // temporary
2951
2952         if( GetDebug() > 1 )
2953           printf("Out %d \n", bit2);
2954
2955         fhPtDecay->Fill(photon2->Pt());
2956
2957         //Fill some histograms about shower shape
2958         if(fFillSelectClHisto && cluster2 && GetReader()->GetDataType()!=AliCaloTrackReader::kMC)
2959           FillSelectedClusterHistograms(cluster2, mom2.Pt(), nMaxima2, photon2->GetTag());
2960
2961         if(IsDataMC())
2962         {
2963           Int_t mcIndex2 = GetMCIndex(photon2->GetTag());
2964           fhMCPtDecay[mcIndex2]->Fill(photon2->Pt());
2965         }
2966       }
2967
2968       //Mass of selected pairs
2969       fhSelectedMass  ->Fill( epair,mass);
2970       fhSelectedMassPt->Fill(ptpair,mass);
2971       if(IsDataMC())fhMCSelectedMassPt[mcIndex]->Fill(ptpair,mass);
2972
2973       // Fill histograms to undertand pile-up before other cuts applied
2974       // Remember to relax time cuts in the reader
2975       FillPileUpHistograms(ptpair,((cluster1->GetTOF()+cluster2->GetTOF())*1e9)/2,cluster1);
2976
2977       //Create AOD for analysis
2978       AliAODPWG4Particle pi0 = AliAODPWG4Particle(mom);
2979
2980       if     ( (GetNeutralMesonSelection()->GetParticle()).Contains("Pi0") ) pi0.SetIdentifiedParticleType(AliCaloPID::kPi0);
2981       else if( (GetNeutralMesonSelection()->GetParticle()).Contains("Eta") ) pi0.SetIdentifiedParticleType(AliCaloPID::kEta);
2982       else
2983       {
2984         printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - Particle type declared in AliNeutralMeson not correct, do not add \n");
2985         return ;
2986       }
2987       pi0.SetDetector(photon1->GetDetector());
2988
2989       // MC
2990       pi0.SetLabel(label);
2991       pi0.SetTag(tag);
2992
2993       //Set the indeces of the original caloclusters
2994       pi0.SetCaloLabel(photon1->GetCaloLabel(0), photon2->GetCaloLabel(0));
2995       //pi0.SetInputFileIndex(input);
2996
2997       AddAODParticle(pi0);
2998
2999     }//2n photon loop
3000
3001   }//1st photon loop
3002
3003   if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - End fill AODs \n");
3004
3005 }
3006
3007 //__________________________________________________
3008 void  AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS()
3009 {
3010   //Do analysis and fill aods
3011   //Search for the photon decay in calorimeters
3012   //Read photon list from AOD, produced in class AliAnaPhoton and AliGammaConversion
3013   //Check if 2 photons have the mass of the pi0.
3014
3015   TLorentzVector mom1;
3016   TLorentzVector mom2;
3017   TLorentzVector mom ;
3018   Int_t tag   = 0;
3019   Int_t label = 0;
3020   Int_t evtIndex = 0;
3021
3022   // Check calorimeter input
3023   if(!GetInputAODBranch())
3024   {
3025     printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - No input calo photons in AOD branch with name < %s > , STOP\n",GetInputAODName().Data());
3026     abort();
3027   }
3028
3029   // Get the array with conversion photons
3030   TClonesArray * inputAODGammaConv = (TClonesArray *) GetReader()->GetOutputEvent()->FindListObject(fInputAODGammaConvName);
3031   if(!inputAODGammaConv)
3032   {
3033     inputAODGammaConv = (TClonesArray *) GetReader()->GetInputEvent()->FindListObject(fInputAODGammaConvName);
3034
3035     if(!inputAODGammaConv)
3036     {
3037       printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - No input gamma conversions in AOD branch with name < %s >\n",fInputAODGammaConvName.Data());
3038
3039       return;
3040     }
3041   }
3042
3043   //Get shower shape information of clusters
3044   TObjArray *clusters = 0;
3045   if     (fCalorimeter=="EMCAL") clusters = GetEMCALClusters();
3046   else if(fCalorimeter=="PHOS")  clusters = GetPHOSClusters() ;
3047
3048   Int_t nCTS  = inputAODGammaConv->GetEntriesFast();
3049   Int_t nCalo = GetInputAODBranch()->GetEntriesFast();
3050   if(nCTS<=0 || nCalo <=0)
3051   {
3052     if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - nCalo %d, nCTS %d, cannot loop\n",nCalo,nCTS);
3053     return;
3054   }
3055
3056   if(GetDebug() > 1)
3057     printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - Number of conversion photons %d\n",nCTS);
3058
3059   // Do the loop, first calo, second CTS
3060   for(Int_t iphoton = 0; iphoton < GetInputAODBranch()->GetEntriesFast(); iphoton++)
3061   {
3062     AliAODPWG4Particle * photon1 =  (AliAODPWG4Particle*) (GetInputAODBranch()->At(iphoton));
3063     mom1 = *(photon1->Momentum());
3064
3065     // Do analysis only when one of the decays is isolated
3066     // Run AliAnaParticleIsolation before
3067     if(fSelectIsolatedDecay)
3068     {
3069       Bool_t isolated1 = ((AliAODPWG4ParticleCorrelation*) photon1)->IsIsolated();
3070       if(!isolated1) continue;
3071     }
3072
3073     //Get original cluster, to recover some information
3074     Int_t iclus = -1;
3075     AliVCluster *cluster = FindCluster(clusters,photon1->GetCaloLabel(0),iclus);
3076
3077     for(Int_t jphoton = 0; jphoton < nCTS; jphoton++)
3078     {
3079       AliAODPWG4Particle * photon2 =  (AliAODPWG4Particle*) (inputAODGammaConv->At(jphoton));
3080
3081       if(GetMixedEvent())
3082         evtIndex = GetMixedEvent()->EventIndexForCaloCluster(photon2->GetCaloLabel(0)) ;
3083       if(TMath::Abs(GetVertex(evtIndex)[2]) > GetZvertexCut()) continue ;  //vertex cut
3084
3085       mom2 = *(photon2->Momentum());
3086
3087       mom = mom1+mom2;
3088
3089       Double_t mass  = mom.M();
3090       Double_t epair = mom.E();
3091       Float_t ptpair = mom.Pt();
3092
3093       Int_t nMaxima = photon1->GetFiducialArea();
3094       if(fFillAllNLMHistograms)
3095       {
3096         if     (nMaxima==1) fhMassPairLocMax[0]->Fill(epair,mass);
3097         else if(nMaxima==2) fhMassPairLocMax[1]->Fill(epair,mass);
3098         else                fhMassPairLocMax[2]->Fill(epair,mass);
3099       }
3100
3101       if(nMaxima < fNLMCutMin || nMaxima > fNLMCutMax) continue ;
3102       if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - NLM %d of out of range \n",nMaxima);
3103
3104       //Play with the MC stack if available
3105       Int_t mcIndex = kmcHadron;
3106       if(IsDataMC())
3107       {
3108         Int_t   label2 = photon2->GetLabel();
3109         if(label2 >= 0 )photon2->SetTag(GetMCAnalysisUtils()->CheckOrigin(label2, GetReader()));
3110
3111         HasPairSameMCMother(photon1, photon2, label, tag) ;
3112         mcIndex = GetMCIndex(tag);
3113       }
3114
3115       //Mass of selected pairs
3116       fhMass  ->Fill( epair,mass);
3117       fhMassPt->Fill(ptpair,mass);
3118       if(IsDataMC()) fhMCMassPt[mcIndex]->Fill(ptpair,mass);
3119
3120       //
3121       // Select good pair (good phi, pt cuts, aperture and invariant mass)
3122       //
3123       if(!GetNeutralMesonSelection()->SelectPair(mom1, mom2,fCalorimeter)) continue ;
3124
3125       if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - Selected gamma pair: pt %f, phi %f, eta%f\n",
3126                                 mom.Pt(), mom.Phi()*TMath::RadToDeg(), mom.Eta());
3127
3128       //
3129       // Tag both photons as decay if not done before
3130       // set the corresponding bit for pi0 or eta or "side" case
3131       //
3132       Int_t bit1 = photon1->GetBtag(); // temporary
3133       if( bit1 < 0 ) bit1 = 0 ; // temporary
3134       if( !GetNeutralMesonSelection()->CheckDecayBit(bit1) )
3135       {
3136         photon1->SetTagged(kTRUE); // temporary
3137         GetNeutralMesonSelection()->SetDecayBit(bit1);
3138         photon1->SetBtag(bit1); // temporary
3139         fhPtDecay->Fill(photon1->Pt());
3140
3141         //Fill some histograms about shower shape
3142         if(fFillSelectClHisto && cluster && GetReader()->GetDataType()!=AliCaloTrackReader::kMC)
3143           FillSelectedClusterHistograms(cluster, mom1.Pt(), nMaxima, photon1->GetTag());
3144
3145         if(IsDataMC())
3146         {
3147           Int_t mcIndex1 = GetMCIndex(photon1->GetTag());
3148           fhMCPtDecay[mcIndex1]->Fill(photon1->Pt());
3149         }
3150       }
3151
3152       Int_t bit2 = photon2->GetBtag(); // temporary
3153       if( bit2 < 0 ) bit2 = 0 ; // temporary
3154       if( !GetNeutralMesonSelection()->CheckDecayBit(bit2) )
3155       {
3156         photon2->SetTagged(kTRUE); // temporary
3157         GetNeutralMesonSelection()->SetDecayBit(bit2);
3158         photon2->SetBtag(bit2); // temporary
3159       }
3160
3161       //Mass of selected pairs
3162       fhSelectedMass  ->Fill( epair,mass);
3163       fhSelectedMassPt->Fill(ptpair,mass);
3164       if(IsDataMC()) fhMCSelectedMassPt[mcIndex]->Fill(ptpair,mass);
3165
3166       // Fill histograms to undertand pile-up before other cuts applied
3167       // Remember to relax time cuts in the reader
3168       if(cluster) FillPileUpHistograms(mom.Pt(),cluster->GetTOF()*1e9,cluster);
3169
3170       //Create AOD for analysis
3171
3172       AliAODPWG4Particle pi0 = AliAODPWG4Particle(mom);
3173
3174       if     ( (GetNeutralMesonSelection()->GetParticle()).Contains("Pi0") ) pi0.SetIdentifiedParticleType(AliCaloPID::kPi0);
3175       else if( (GetNeutralMesonSelection()->GetParticle()).Contains("Eta") ) pi0.SetIdentifiedParticleType(AliCaloPID::kEta);
3176       else
3177       {
3178         printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - Particle type declared in AliNeutralMeson not correct, do not add \n");
3179         return ;
3180       }
3181       pi0.SetDetector(photon1->GetDetector());
3182
3183       // MC
3184       pi0.SetLabel(label);
3185       pi0.SetTag(tag);
3186
3187       //Set the indeces of the original tracks or caloclusters
3188       pi0.SetCaloLabel (photon1->GetCaloLabel(0) , -1);
3189       pi0.SetTrackLabel(photon2->GetTrackLabel(0), photon2->GetTrackLabel(1));
3190       //pi0.SetInputFileIndex(input);
3191
3192       AddAODParticle(pi0);
3193
3194     }//2n photon loop
3195
3196   }//1st photon loop
3197
3198   if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - End fill AODs \n");
3199
3200 }
3201
3202 //_________________________________________________
3203 void  AliAnaPi0EbE::MakeShowerShapeIdentification()
3204 {
3205   //Search for pi0 in fCalorimeter with shower shape analysis
3206
3207   TObjArray * pl        = 0x0;
3208   AliVCaloCells * cells = 0x0;
3209   //Select the Calorimeter of the photon
3210   if      (fCalorimeter == "EMCAL" )
3211   {
3212     pl    = GetEMCALClusters();
3213     cells = GetEMCALCells();
3214   }
3215   else if (fCalorimeter == "PHOS")
3216   {
3217     AliFatal("kSSCalo case not implememted for PHOS");
3218     return; // for coverity
3219
3220     //pl    = GetPHOSClusters();
3221     //cells = GetPHOSCells();
3222   }
3223
3224   if(!pl)
3225   {
3226     Info("MakeShowerShapeIdentification","TObjArray with %s clusters is NULL!\n",fCalorimeter.Data());
3227     return;
3228   }
3229
3230   TLorentzVector mom ;
3231   for(Int_t icalo = 0; icalo < pl->GetEntriesFast(); icalo++)
3232   {
3233     AliVCluster * calo = (AliVCluster*) (pl->At(icalo));
3234
3235     Int_t evtIndex = 0 ;
3236     if (GetMixedEvent())
3237     {
3238       evtIndex=GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ;
3239
3240     if(TMath::Abs(GetVertex(evtIndex)[2]) > GetZvertexCut()) continue ;  //vertex cut
3241     }
3242
3243     //Get Momentum vector,
3244     Double_t vertex[]={0,0,0};
3245     if(GetReader()->GetDataType() != AliCaloTrackReader::kMC)
3246     {
3247       calo->GetMomentum(mom,GetVertex(evtIndex)) ;
3248     }//Assume that come from vertex in straight line
3249     else
3250     {
3251       calo->GetMomentum(mom,vertex) ;
3252     }
3253
3254     //If too small or big pt, skip it
3255     if(mom.E() < GetMinEnergy() || mom.E() > GetMaxEnergy() ) continue ;
3256
3257     //Check acceptance selection
3258     if(IsFiducialCutOn())
3259     {
3260       Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter) ;
3261       if(! in ) continue ;
3262     }
3263
3264     if(GetDebug() > 1)
3265       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());
3266
3267     //Play with the MC stack if available
3268     //Check origin of the candidates
3269     Int_t tag   = 0 ;
3270     if(IsDataMC())
3271     {
3272       tag = GetMCAnalysisUtils()->CheckOrigin(calo->GetLabels(),calo->GetNLabels(),GetReader());
3273       //GetMCAnalysisUtils()->CheckMultipleOrigin(calo->GetLabels(),calo->GetNLabels(), GetReader(), aodpi0.GetInputFileIndex(), tag);
3274       if(GetDebug() > 0) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - Origin of candidate %d\n",tag);
3275     }
3276
3277     //Int_t nMaxima = GetCaloUtils()->GetNumberOfLocalMaxima(calo, cells); // NLM
3278
3279     //Check Distance to Bad channel, set bit.
3280     Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
3281     if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;
3282     if(distBad < fMinDist){ //In bad channel (PHOS cristal size 2.2x2.2 cm)
3283       //FillRejectedClusterHistograms(mom,tag,nMaxima);
3284       continue ;
3285     }
3286
3287     if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - FillAOD: Bad channel cut passed %4.2f\n",distBad);
3288
3289     //If too low number of cells, skip it
3290     if ( calo->GetNCells() < GetCaloPID()->GetClusterSplittingMinNCells())
3291     {
3292       //FillRejectedClusterHistograms(mom,tag,nMaxima);
3293       continue ;
3294     }
3295
3296     if(GetDebug() > 1)
3297       printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - FillAOD: N cells cut passed %d > %d\n",
3298              calo->GetNCells(), GetCaloPID()->GetClusterSplittingMinNCells());
3299
3300     //.......................................
3301     // TOF cut, BE CAREFUL WITH THIS CUT
3302     Double_t tof = calo->GetTOF()*1e9;
3303     if(tof < fTimeCutMin || tof > fTimeCutMax)
3304     {
3305       //FillRejectedClusterHistograms(mom,tag,nMaxima);
3306       continue ;
3307     }
3308
3309     //Check PID
3310     //PID selection or bit setting
3311     Int_t    nMaxima  = 0;
3312     Double_t mass     = 0, angle    = 0;
3313     Int_t    absId1   =-1, absId2   =-1;
3314     Float_t  distbad1 =-1, distbad2 =-1;
3315     Bool_t   fidcut1  = 0, fidcut2  = 0;
3316     TLorentzVector    l1, l2;
3317
3318     Int_t idPartType = GetCaloPID()->GetIdentifiedParticleTypeFromClusterSplitting(calo,cells,GetCaloUtils(),
3319                                                                                    GetVertex(evtIndex),nMaxima,
3320                                                                                    mass,angle,l1,l2,absId1,absId2,
3321                                                                                    distbad1,distbad2,fidcut1,fidcut2) ;
3322
3323
3324     if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - PDG of identified particle %d\n",idPartType);
3325
3326     // Skip events where one of the new clusters (lowest energy) is close to an EMCal border or a bad channel
3327     if( (fCheckSplitDistToBad) &&
3328        (!fidcut2 || !fidcut1 || distbad1 < fMinDist || distbad2 < fMinDist))
3329     {
3330       if(GetDebug() > 1)
3331         Info("MakeShowerShapeIdentification", "Dist to bad channel cl %f, cl1 %f, cl2 %f; fid cl1 %d, cl2 %d \n",
3332                calo->GetDistanceToBadChannel(),distbad1,distbad2, fidcut1,fidcut2);
3333
3334       //FillRejectedClusterHistograms(mom,tag,nMaxima);
3335       continue ;
3336     }
3337
3338     //Skip events with too few or too many  NLM
3339     if(nMaxima < fNLMCutMin || nMaxima > fNLMCutMax)
3340     {
3341       //FillRejectedClusterHistograms(mom,tag,nMaxima);
3342       continue ;
3343     }
3344
3345     if(GetDebug() > 1)
3346       printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - NLM %d accepted \n",nMaxima);
3347
3348     //Skip matched clusters with tracks
3349     if(fRejectTrackMatch && IsTrackMatched(calo, GetReader()->GetInputEvent()))
3350     {
3351       FillRejectedClusterHistograms(mom,tag,nMaxima);
3352       continue ;
3353     }
3354
3355     Float_t l0 = calo->GetM02();
3356     Float_t e1 = l1.Energy();
3357     Float_t e2 = l2.Energy();
3358     TLorentzVector l12 = l1+l2;
3359     Float_t ptSplit = l12.Pt();
3360     Float_t  eSplit = e1+e2;
3361
3362     //mass of all clusters
3363     fhMass       ->Fill(mom.E() ,mass);
3364     fhMassPt     ->Fill(mom.Pt(),mass);
3365     fhMassSplitPt->Fill(ptSplit ,mass);
3366     fhPtLambda0NoSplitCut->Fill(mom.Pt(),l0);
3367
3368     // Asymmetry of all clusters
3369     Float_t asy =-10;
3370
3371     if(e1+e2 > 0) asy = (e1-e2) / (e1+e2);
3372     fhAsymmetry->Fill(mom.E(),asy);
3373
3374     // Divide NLM in 3 cases, 1 local maxima, 2 local maxima, more than 2 local maxima
3375     Int_t indexMax = -1;
3376     if     (nMaxima==1) indexMax = 0 ;
3377     else if(nMaxima==2) indexMax = 1 ;
3378     else                indexMax = 2 ;
3379     fhMassPtLocMax[indexMax]->Fill(mom.Pt(),mass);
3380
3381     Int_t   mcIndex   =-1;
3382     Int_t   noverlaps = 0;
3383     Float_t ptprim    = 0;
3384     if(IsDataMC())
3385     {
3386       mcIndex = GetMCIndex(tag);
3387
3388       Bool_t ok      = kFALSE;
3389       Int_t  mcLabel = calo->GetLabel();
3390
3391       TLorentzVector primary = GetMCAnalysisUtils()->GetMother(mcLabel,GetReader(),ok);
3392
3393       Int_t mesonLabel = -1;
3394
3395       if(mcIndex == kmcPi0 || mcIndex == kmcEta)
3396       {
3397         if(mcIndex == kmcPi0)
3398         {
3399           TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(mcLabel,111,GetReader(),ok,mesonLabel);
3400           if(grandmom.E() > 0 && ok) ptprim =  grandmom.Pt();
3401         }
3402         else
3403         {
3404           TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(mcLabel,221,GetReader(),ok,mesonLabel);
3405           if(grandmom.E() > 0 && ok) ptprim =  grandmom.Pt();
3406         }
3407       }
3408
3409       const UInt_t nlabels = calo->GetNLabels();
3410       Int_t overpdg[nlabels];
3411       noverlaps = GetMCAnalysisUtils()->GetNOverlaps(calo->GetLabels(), nlabels,tag,mesonLabel,GetReader(),overpdg);
3412
3413       fhMCMassPt     [mcIndex]->Fill(mom.Pt(),mass);
3414       fhMCMassSplitPt[mcIndex]->Fill(ptSplit ,mass);
3415       if(mcIndex==kmcPi0)
3416       {
3417         fhMCPi0PtRecoPtPrim                     ->Fill(mom.Pt(),ptprim);
3418         fhMCPi0SplitPtRecoPtPrim                ->Fill(ptSplit ,ptprim);
3419         fhMCPi0PtRecoPtPrimLocMax     [indexMax]->Fill(mom.Pt(),ptprim);
3420         fhMCPi0SplitPtRecoPtPrimLocMax[indexMax]->Fill(ptSplit ,ptprim);
3421
3422       }
3423       else if(mcIndex==kmcEta)
3424       {
3425         fhMCEtaPtRecoPtPrim                     ->Fill(mom.Pt(),ptprim);
3426         fhMCEtaSplitPtRecoPtPrim                ->Fill(ptSplit ,ptprim);
3427         fhMCEtaPtRecoPtPrimLocMax     [indexMax]->Fill(mom.Pt(),ptprim);
3428         fhMCEtaSplitPtRecoPtPrimLocMax[indexMax]->Fill(ptSplit ,ptprim);
3429       }
3430
3431       if(noverlaps==0)
3432       {
3433         if(mcIndex==kmcPi0)
3434         {
3435           fhMCPi0PtRecoPtPrimNoOverlap     ->Fill(mom.Pt(),ptprim);
3436           fhMCPi0SplitPtRecoPtPrimNoOverlap->Fill(ptSplit ,ptprim);
3437         }
3438         else if(mcIndex==kmcEta)
3439         {
3440           fhMCEtaPtRecoPtPrimNoOverlap     ->Fill(mom.Pt(),ptprim);
3441           fhMCEtaSplitPtRecoPtPrimNoOverlap->Fill(ptSplit ,ptprim);
3442         }
3443
3444         fhMassNoOverlap       ->Fill(mom.E() ,mass);
3445         fhMassPtNoOverlap     ->Fill(mom.Pt(),mass);
3446         fhMassSplitPtNoOverlap->Fill(ptSplit ,mass);
3447
3448         fhMCMassPtNoOverlap     [mcIndex]->Fill(mom.Pt(),mass);
3449         fhMCMassSplitPtNoOverlap[mcIndex]->Fill(ptSplit ,mass);
3450       }
3451
3452       fhMCPtAsymmetry[mcIndex]->Fill(mom.Pt(),asy);
3453     }
3454
3455     // If cluster does not pass pid, not pi0/eta, skip it.
3456     if     (GetOutputAODName().Contains("Pi0") && idPartType != AliCaloPID::kPi0)
3457     {
3458       if(GetDebug() > 1) Info("MakeShowerShapeIdentification","Cluster is not Pi0\n");
3459       FillRejectedClusterHistograms(mom,tag,nMaxima);
3460       continue ;
3461     }
3462
3463     else if(GetOutputAODName().Contains("Eta") && idPartType != AliCaloPID::kEta)
3464     {
3465       if(GetDebug() > 1) Info("MakeShowerShapeIdentification","Cluster is not Eta\n");
3466       FillRejectedClusterHistograms(mom,tag,nMaxima);
3467       continue ;
3468     }
3469
3470     if(GetDebug() > 1)
3471       Info("MakeShowerShapeIdentification","Pi0/Eta selection cuts passed: pT %3.2f, pdg %d\n",
3472              mom.Pt(), idPartType);
3473
3474     //Mass and asymmetry of selected pairs
3475     fhSelectedAsymmetry  ->Fill(mom.E() ,asy );
3476     fhSelectedMass       ->Fill(mom.E() ,mass);
3477     fhSelectedMassPt     ->Fill(mom.Pt(),mass);
3478     fhSelectedMassSplitPt->Fill(ptSplit ,mass);
3479     if(fFillAllNLMHistograms) fhSelectedMassPtLocMax[indexMax]->Fill(mom.Pt(),mass);
3480
3481     Int_t   nSM  = GetModuleNumber(calo);
3482     if(nSM < GetCaloUtils()->GetNumberOfSuperModulesUsed() && nSM >=0 && fFillAllNLMHistograms)
3483     {
3484       fhSelectedMassPtLocMaxSM   [indexMax][nSM]->Fill(mom.Pt(),mass);
3485       fhSelectedLambda0PtLocMaxSM[indexMax][nSM]->Fill(mom.Pt(),l0  );
3486     }
3487
3488     if(IsDataMC())
3489     {
3490       if(mcIndex==kmcPi0)
3491       {
3492         fhMCPi0SelectedPtRecoPtPrim                     ->Fill(mom.Pt(),ptprim);
3493         fhMCPi0SelectedSplitPtRecoPtPrim                ->Fill(ptSplit ,ptprim);
3494         fhMCPi0SelectedPtRecoPtPrimLocMax     [indexMax]->Fill(mom.Pt(),ptprim);
3495         fhMCPi0SelectedSplitPtRecoPtPrimLocMax[indexMax]->Fill(ptSplit ,ptprim);
3496       }
3497       else if(mcIndex==kmcEta)
3498       {
3499         fhMCEtaSelectedPtRecoPtPrim                     ->Fill(mom.Pt(),ptprim);
3500         fhMCEtaSelectedSplitPtRecoPtPrim                ->Fill(ptSplit ,ptprim);
3501         fhMCEtaSelectedPtRecoPtPrimLocMax     [indexMax]->Fill(mom.Pt(),ptprim);
3502         fhMCEtaSelectedSplitPtRecoPtPrimLocMax[indexMax]->Fill(ptSplit ,ptprim);
3503       }
3504
3505       if(noverlaps==0)
3506       {
3507         fhSelectedMassNoOverlap       ->Fill(mom.E() ,mass);
3508         fhSelectedMassPtNoOverlap     ->Fill(mom.Pt(),mass);
3509         fhSelectedMassSplitPtNoOverlap->Fill(ptSplit ,mass);
3510
3511         if(mcIndex==kmcPi0)
3512         {
3513           fhMCPi0SelectedPtRecoPtPrimNoOverlap     ->Fill(mom.Pt(),ptprim);
3514           fhMCPi0SelectedSplitPtRecoPtPrimNoOverlap->Fill(ptSplit ,ptprim);
3515         }
3516         else if(mcIndex==kmcEta)
3517         {
3518           fhMCEtaSelectedPtRecoPtPrimNoOverlap     ->Fill(mom.Pt(),ptprim);
3519           fhMCEtaSelectedSplitPtRecoPtPrimNoOverlap->Fill(ptSplit ,ptprim);
3520         }
3521       }
3522     }
3523
3524     fhSplitE        ->Fill( eSplit);
3525     fhSplitPt       ->Fill(ptSplit);
3526     Float_t phi = mom.Phi();
3527     if(phi<0) phi+=TMath::TwoPi();
3528     fhSplitPtPhi    ->Fill(ptSplit,phi);
3529     fhSplitPtEta    ->Fill(ptSplit,mom.Eta());
3530     fhNLocMaxSplitPt->Fill(ptSplit ,nMaxima);
3531
3532     //Check split-clusters with good time window difference
3533     Double_t tof1  = cells->GetCellTime(absId1);
3534     GetCaloUtils()->RecalibrateCellTime(tof1, fCalorimeter, absId1,GetReader()->GetInputEvent()->GetBunchCrossNumber());
3535     tof1*=1.e9;
3536
3537     Double_t tof2  = cells->GetCellTime(absId2);
3538     GetCaloUtils()->RecalibrateCellTime(tof2, fCalorimeter, absId2,GetReader()->GetInputEvent()->GetBunchCrossNumber());
3539     tof2*=1.e9;
3540
3541     Double_t t12diff = tof1-tof2;
3542     fhEPairDiffTime->Fill(e1+e2,    t12diff);
3543
3544     if(IsDataMC())
3545     {
3546       fhMCSplitE        [mcIndex]->Fill( eSplit);
3547       fhMCSplitPt       [mcIndex]->Fill(ptSplit);
3548       fhMCSplitPtPhi    [mcIndex]->Fill(ptSplit,phi);
3549       fhMCSplitPtEta    [mcIndex]->Fill(ptSplit,mom.Eta());
3550       fhMCNLocMaxSplitPt[mcIndex]->Fill(ptSplit ,nMaxima);
3551       fhMCNLocMaxPt     [mcIndex]->Fill(mom.Pt(),nMaxima);
3552
3553       fhMCSelectedMassPt     [mcIndex]->Fill(mom.Pt(),mass);
3554       fhMCSelectedMassSplitPt[mcIndex]->Fill(ptSplit,mass);
3555       fhMCSelectedMassPtLocMax[mcIndex][indexMax]->Fill(mom.Pt(),mass);
3556
3557       if(noverlaps==0)
3558       {
3559         fhMCSelectedMassPtNoOverlap     [mcIndex]->Fill(mom.Pt(),mass);
3560         fhMCSelectedMassSplitPtNoOverlap[mcIndex]->Fill(ptSplit,mass);
3561       }
3562     }
3563
3564     // Remove clusters with NLM=x depeding on a minimim energy cut
3565     if(nMaxima == 1 && fNLMECutMin[0] > mom.E()) continue;
3566     if(nMaxima == 2 && fNLMECutMin[1] > mom.E()) continue;
3567     if(nMaxima >  2 && fNLMECutMin[2] > mom.E()) continue;
3568
3569     //Fill some histograms about shower shape
3570     if(fFillSelectClHisto && GetReader()->GetDataType()!=AliCaloTrackReader::kMC)
3571     {
3572       FillSelectedClusterHistograms(calo, mom.Pt(), nMaxima, tag, asy);
3573     }
3574
3575     // Fill histograms to undertand pile-up before other cuts applied
3576     // Remember to relax time cuts in the reader
3577     Double_t tofcluster   = calo->GetTOF()*1e9;
3578
3579     FillPileUpHistograms(mom.Pt(),tofcluster,calo);
3580
3581     if(fFillEMCALBCHistograms && fCalorimeter=="EMCAL")
3582       FillEMCALBCHistograms(mom.E(), mom.Eta(), mom.Phi(), tofcluster);
3583
3584     //-----------------------
3585     //Create AOD for analysis
3586
3587     AliAODPWG4Particle aodpi0 = AliAODPWG4Particle(mom);
3588     aodpi0.SetLabel(calo->GetLabel());
3589
3590     //Set the indeces of the original caloclusters
3591     aodpi0.SetCaloLabel(calo->GetID(),-1);
3592     aodpi0.SetDetector(fCalorimeter);
3593
3594     if     (distBad > fMinDist3) aodpi0.SetDistToBad(2) ;
3595     else if(distBad > fMinDist2) aodpi0.SetDistToBad(1) ;
3596     else                         aodpi0.SetDistToBad(0) ;
3597
3598     // Check if cluster is pi0 via cluster splitting
3599     aodpi0.SetIdentifiedParticleType(idPartType);
3600
3601     // Add number of local maxima to AOD, method name in AOD to be FIXED
3602     aodpi0.SetFiducialArea(nMaxima);
3603
3604     aodpi0.SetTag(tag);
3605
3606     //Add AOD with pi0 object to aod branch
3607     AddAODParticle(aodpi0);
3608
3609   }//loop
3610
3611   if(GetDebug() > 1) Info("MakeShowerShapeIdentification","End fill AODs \n");
3612
3613 }
3614
3615 //______________________________________________
3616 void  AliAnaPi0EbE::MakeAnalysisFillHistograms()
3617 {
3618   //Do analysis and fill histograms
3619
3620   if(!GetOutputAODBranch())
3621   {
3622     AliFatal(Form("No output pi0 in AOD branch with name < %s >,STOP \n",GetOutputAODName().Data()));
3623   }
3624
3625   //Loop on stored AOD pi0
3626   Int_t naod = GetOutputAODBranch()->GetEntriesFast();
3627   if(GetDebug() > 0) Info("MakeAnalysisFillHistograms","aod branch entries %d\n", naod);
3628
3629   Float_t cen = GetEventCentrality();
3630   Float_t ep  = GetEventPlaneAngle();
3631
3632   for(Int_t iaod = 0; iaod < naod ; iaod++)
3633   {
3634     AliAODPWG4Particle* pi0 =  (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
3635     Int_t pdg = pi0->GetIdentifiedParticleType();
3636
3637     if( ( pdg != AliCaloPID::kPi0 && pdg != AliCaloPID::kEta ) ) continue;
3638
3639     //Fill pi0 histograms
3640     Float_t ener  = pi0->E();
3641     Float_t pt    = pi0->Pt();
3642     Float_t phi   = pi0->Phi();
3643     if(phi < 0) phi+=TMath::TwoPi();
3644     Float_t eta = pi0->Eta();
3645
3646     fhPt     ->Fill(pt  );
3647     fhE      ->Fill(ener);
3648
3649     fhPtEta  ->Fill(pt  ,eta);
3650     fhPtPhi  ->Fill(pt  ,phi);
3651     fhEtaPhi ->Fill(eta ,phi);
3652
3653     if(fFillHighMultHistograms)
3654     {
3655       fhPtCentrality ->Fill(pt,cen) ;
3656       fhPtEventPlane ->Fill(pt,ep ) ;
3657     }
3658
3659     if(IsDataMC())
3660     {
3661       Int_t tag     = pi0->GetTag();
3662       Int_t label   = pi0->GetLabel();
3663       Int_t mcIndex = GetMCIndex(tag);
3664
3665       fhMCE    [mcIndex] ->Fill(ener);
3666       fhMCPt   [mcIndex] ->Fill(pt);
3667       fhMCPtPhi[mcIndex] ->Fill(pt,phi);
3668       fhMCPtEta[mcIndex] ->Fill(pt,eta);
3669
3670       if(fFillHighMultHistograms) fhMCPtCentrality[mcIndex]->Fill(pt,cen);
3671
3672       if((mcIndex==kmcPi0Decay || mcIndex==kmcEtaDecay ||
3673           mcIndex==kmcPi0      || mcIndex==kmcEta         ) &&
3674          fAnaType==kSSCalo)
3675       {
3676         Float_t efracMC   = 0;
3677         Int_t   momlabel  = -1;
3678         Bool_t  ok        = kFALSE;
3679
3680         TLorentzVector mom   = GetMCAnalysisUtils()->GetMother(label,GetReader(),ok);
3681         if(!ok) continue;
3682
3683         if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0))
3684         {
3685           TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(label,111,GetReader(),ok,momlabel);
3686           if(grandmom.E() > 0 && ok)
3687           {
3688             efracMC =  grandmom.E()/ener;
3689             fhMCPi0PtGenRecoFraction ->Fill(pt,efracMC);
3690           }
3691         }
3692         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay))
3693         {
3694           fhMCPi0DecayPt->Fill(pt);
3695           TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(label,111,GetReader(),ok,momlabel);
3696           if(grandmom.E() > 0 && ok)
3697           {
3698             efracMC =  mom.E()/grandmom.E();
3699             fhMCPi0DecayPtFraction ->Fill(pt,efracMC);
3700           }
3701         }
3702         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta))
3703         {
3704           TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(label,221,GetReader(),ok,momlabel);
3705           if(grandmom.E() > 0 && ok)
3706           {
3707             efracMC =  grandmom.E()/ener;
3708             fhMCEtaPtGenRecoFraction ->Fill(pt,efracMC);
3709           }
3710         }
3711         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay))
3712         {
3713           fhMCEtaDecayPt->Fill(pt);
3714           TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(label,221,GetReader(),ok,momlabel);
3715           if(grandmom.E() > 0 && ok)
3716           {
3717             efracMC =  mom.E()/grandmom.E();
3718             fhMCEtaDecayPtFraction ->Fill(pt,efracMC);
3719           }
3720         }
3721         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay))
3722         {
3723           fhMCOtherDecayPt->Fill(pt);
3724         }
3725
3726       }
3727
3728       if( mcIndex==kmcPi0 || mcIndex==kmcEta )
3729       {
3730         Float_t prodR     = -1;
3731         Int_t   momindex  = -1;
3732         Int_t   mompdg    = -1;
3733         Int_t   momstatus = -1;
3734
3735         if(GetReader()->ReadStack())
3736         {
3737           TParticle* ancestor = GetMCStack()->Particle(label);
3738           momindex  = ancestor->GetFirstMother();
3739           if(momindex < 0) return;
3740           TParticle* mother = GetMCStack()->Particle(momindex);
3741           mompdg    = TMath::Abs(mother->GetPdgCode());
3742           momstatus = mother->GetStatusCode();
3743           prodR = mother->R();
3744         }
3745         else
3746         {
3747           TClonesArray * mcparticles = GetReader()->GetAODMCParticles();
3748           AliAODMCParticle* ancestor = (AliAODMCParticle *) mcparticles->At(label);
3749           momindex  = ancestor->GetMother();
3750           if(momindex < 0) return;
3751           AliAODMCParticle* mother = (AliAODMCParticle *) mcparticles->At(momindex);
3752           mompdg    = TMath::Abs(mother->GetPdgCode());
3753           momstatus = mother->GetStatus();
3754           prodR = TMath::Sqrt(mother->Xv()*mother->Xv()+mother->Yv()*mother->Yv());
3755         }
3756
3757         if( mcIndex==kmcPi0 )
3758         {
3759           fhMCPi0ProdVertex->Fill(pt,prodR);
3760
3761           if     (momstatus  == 21) fhMCPi0PtOrigin->Fill(pt,0.5);//parton
3762           else if(mompdg     < 22 ) fhMCPi0PtOrigin->Fill(pt,1.5);//quark
3763           else if(mompdg     > 2100  && mompdg   < 2210) fhMCPi0PtOrigin->Fill(pt,2.5);// resonances
3764           else if(mompdg    == 221) fhMCPi0PtOrigin->Fill(pt,8.5);//eta
3765           else if(mompdg    == 331) fhMCPi0PtOrigin->Fill(pt,9.5);//eta prime
3766           else if(mompdg    == 213) fhMCPi0PtOrigin->Fill(pt,4.5);//rho
3767           else if(mompdg    == 223) fhMCPi0PtOrigin->Fill(pt,5.5);//omega
3768           else if(mompdg    >= 310   && mompdg    <= 323) fhMCPi0PtOrigin->Fill(pt,6.5);//k0S, k+-,k*
3769           else if(mompdg    == 130) fhMCPi0PtOrigin->Fill(pt,6.5);//k0L
3770           else if(momstatus == 11 || momstatus  == 12 ) fhMCPi0PtOrigin->Fill(pt,3.5);//resonances
3771           else                      fhMCPi0PtOrigin->Fill(pt,7.5);//other?
3772         }
3773         else if (mcIndex==kmcEta )
3774         {
3775           fhMCEtaProdVertex->Fill(pt,prodR);
3776
3777           if     (momstatus  == 21) fhMCEtaPtOrigin->Fill(pt,0.5);//parton
3778           else if(mompdg     < 22 ) fhMCEtaPtOrigin->Fill(pt,1.5);//quark
3779           else if(mompdg     > 2100  && mompdg   < 2210) fhMCEtaPtOrigin->Fill(pt,2.5);// resonances
3780           else if(mompdg    == 221) fhMCEtaPtOrigin->Fill(pt,8.5);//eta
3781           else if(mompdg    == 331) fhMCEtaPtOrigin->Fill(pt,9.5);//eta prime
3782           else if(mompdg    == 213) fhMCEtaPtOrigin->Fill(pt,4.5);//rho
3783           else if(mompdg    == 223) fhMCEtaPtOrigin->Fill(pt,5.5);//omega
3784           else if(mompdg    >= 310   && mompdg    <= 323) fhMCEtaPtOrigin->Fill(pt,6.5);//k0S, k+-,k*
3785           else if(mompdg    == 130) fhMCEtaPtOrigin->Fill(pt,6.5);//k0L
3786           else if(momstatus == 11 || momstatus  == 12 ) fhMCEtaPtOrigin->Fill(pt,3.5);//resonances
3787           else                      fhMCEtaPtOrigin->Fill(pt,7.5);//other?
3788         }
3789       }
3790
3791     }//Histograms with MC
3792
3793   }// aod loop
3794
3795 }
3796
3797 //__________________________________________________________________
3798 void AliAnaPi0EbE::Print(const Option_t * opt) const
3799 {
3800   //Print some relevant parameters set for the analysis
3801   if(! opt)
3802     return;
3803
3804   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
3805   AliAnaCaloTrackCorrBaseClass::Print("");
3806   printf("Analysis Type = %d \n",  fAnaType) ;
3807   if(fAnaType == kSSCalo)
3808   {
3809     printf("Calorimeter            =     %s\n", fCalorimeter.Data()) ;
3810     printf("Min Distance to Bad Channel   = %2.1f\n",fMinDist);
3811     printf("Min Distance to Bad Channel 2 = %2.1f\n",fMinDist2);
3812     printf("Min Distance to Bad Channel 3 = %2.1f\n",fMinDist3);
3813   }
3814   printf("    \n") ;
3815
3816 }
3817
3818