PWGGA/CaloTrackCorrelations/AliAnaPi0EbE.cxx

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