PWGGA/CaloTrackCorrelations/AliAnaPi0EbE.cxx

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