PWGGA/CaloTrackCorrelations/AliAnaPi0EbE.cxx

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