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   {
2618     return kmcOtherDecay ;
2619   }//decay photon from other than eta or pi0
2620   else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron))
2621   {
2622     return kmcElectron ;
2623   }//electron
2624   else
2625   {
2626     return kmcHadron ;
2627   }//other particles
2628
2629 }
2630
2631 //__________________________________________________________________
2632 void AliAnaPi0EbE::HasPairSameMCMother(Int_t label1 , Int_t label2,
2633                                        Int_t tag1   , Int_t tag2,
2634                                        Int_t & label, Int_t & tag)
2635 {
2636   // Check the labels of pair in case mother was same pi0 or eta
2637   // Set the new AOD accordingly
2638
2639
2640   if(label1 < 0 || label2 < 0 || label1 == label2) return ;
2641
2642   AliDebug(0,Form("Origin of: photon1 %d; photon2 %d",tag1, tag2));
2643
2644   if( (GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCPi0Decay) &&
2645        GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCPi0Decay)    ) ||
2646       (GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCEtaDecay) &&
2647        GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCEtaDecay)    )
2648      )
2649   {
2650     Int_t pdg1    = -1;//, pdg2    = -1;
2651     Int_t ndaugh1 = -1, ndaugh2 = -1;
2652     //Check if pi0/eta mother is the same
2653     if(GetReader()->ReadStack())
2654     {
2655       if(label1>=0)
2656       {
2657         TParticle * mother1 = GetMCStack()->Particle(label1);//photon in kine tree
2658         label1 = mother1->GetFirstMother();
2659         mother1 = GetMCStack()->Particle(label1);//pi0
2660         pdg1=mother1->GetPdgCode();
2661         ndaugh1 = mother1->GetNDaughters();
2662       }
2663       if(label2>=0)
2664       {
2665         TParticle * mother2 = GetMCStack()->Particle(label2);//photon in kine tree
2666         label2 = mother2->GetFirstMother();
2667         mother2 = GetMCStack()->Particle(label2);//pi0
2668         //pdg2=mother2->GetPdgCode();
2669         ndaugh2 = mother2->GetNDaughters();
2670       }
2671     } // STACK
2672     else if(GetReader()->ReadAODMCParticles())
2673     {//&& (input > -1)){
2674       if(label1>=0)
2675       {
2676         AliAODMCParticle * mother1 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles())->At(label1);//photon in kine tree
2677         label1 = mother1->GetMother();
2678         mother1 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles())->At(label1);//pi0
2679         pdg1=mother1->GetPdgCode();
2680         ndaugh1 = mother1->GetNDaughters();
2681       }
2682       if(label2>=0)
2683       {
2684         AliAODMCParticle * mother2 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles())->At(label2);//photon in kine tree
2685         label2 = mother2->GetMother();
2686         mother2 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles())->At(label2);//pi0
2687         //pdg2=mother2->GetPdgCode();
2688         ndaugh2 = mother2->GetNDaughters();
2689       }
2690     }// AOD
2691
2692     //printf("mother1 %d, mother2 %d\n",label1,label2);
2693     if( label1 == label2 && label1>=0  && ndaugh1==ndaugh2 && ndaugh1==2)
2694     {
2695       label = label1;
2696
2697       Double_t angle = fMomentum2.Angle(fMomentum1.Vect());
2698       Double_t mass  = (fMomentum1+fMomentum2).M();
2699       Double_t epair = (fMomentum1+fMomentum2).E();
2700
2701       if(pdg1==111)
2702       {
2703         //printf("Real pi0 pair: pt %f, mass %f\n",(fMomentum1+fMomentum2).Pt(),mass);
2704         fhMassPairMCPi0 ->Fill(epair,mass);
2705         fhAnglePairMCPi0->Fill(epair,angle);
2706         GetMCAnalysisUtils()->SetTagBit(tag,AliMCAnalysisUtils::kMCPi0);
2707 //        printf(" Lab1 %d (%d), lab2 %d (%d), pdg1 %d, pdg2 %d, Is In calo %d, %d, Is lost %d, %d\n",
2708 //               label1,photon1->GetLabel(),label2,photon2->GetLabel(), pdg1, pdg2,
2709 //               GetMCAnalysisUtils()->CheckTagBit(photon1->GetTag(),AliMCAnalysisUtils::kMCDecayPairInCalo),
2710 //               GetMCAnalysisUtils()->CheckTagBit(photon2->GetTag(),AliMCAnalysisUtils::kMCDecayPairInCalo),
2711 //               GetMCAnalysisUtils()->CheckTagBit(photon1->GetTag(),AliMCAnalysisUtils::kMCDecayPairLost),
2712 //               GetMCAnalysisUtils()->CheckTagBit(photon2->GetTag(),AliMCAnalysisUtils::kMCDecayPairLost));
2713       }
2714       else  if(pdg1==221)
2715       {
2716         //printf("Real eta pair\n");
2717         fhMassPairMCEta ->Fill(epair,mass);
2718         fhAnglePairMCEta->Fill(epair,angle);
2719         GetMCAnalysisUtils()->SetTagBit(tag,AliMCAnalysisUtils::kMCEta);
2720       }
2721
2722     } // same label
2723   } // both from eta or pi0 decay
2724
2725 }
2726
2727 //____________________________________________________________________________
2728 void AliAnaPi0EbE::Init()
2729 {
2730   //Init
2731   //Do some checks
2732
2733   if       ( GetCalorimeter() == kPHOS  && !GetReader()->IsPHOSSwitchedOn()  && NewOutputAOD() )
2734     AliFatal("STOP: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!");
2735   else  if ( GetCalorimeter() == kEMCAL && !GetReader()->IsEMCALSwitchedOn() && NewOutputAOD() )
2736     AliFatal("STOP: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!");
2737
2738 }
2739
2740 //____________________________________________________________________________
2741 void AliAnaPi0EbE::InitParameters()
2742 {
2743   //Initialize the parameters of the analysis.
2744   AddToHistogramsName("AnaPi0EbE_");
2745
2746   fInputAODGammaConvName = "PhotonsCTS" ;
2747   fAnaType = kIMCalo ;
2748   fMinDist  = 2.;
2749   fMinDist2 = 4.;
2750   fMinDist3 = 5.;
2751
2752   fNLMECutMin[0] = 10.;
2753   fNLMECutMin[1] = 6. ;
2754   fNLMECutMin[2] = 6. ;
2755
2756   fIsoCandMinPt = 8;
2757   fR            = 0.4;
2758
2759 }
2760
2761 //__________________________________________________________________
2762 void  AliAnaPi0EbE::MakeAnalysisFillAOD()
2763 {
2764   //Do analysis and fill aods
2765
2766   AliDebug(1,Form("Start analysis type %d",fAnaType));
2767
2768   switch(fAnaType)
2769   {
2770     case kIMCalo:
2771       MakeInvMassInCalorimeter();
2772       break;
2773
2774     case kSSCalo:
2775       MakeShowerShapeIdentification();
2776       break;
2777
2778     case kIMCaloTracks:
2779       MakeInvMassInCalorimeterAndCTS();
2780       break;
2781
2782   }
2783
2784   AliDebug(1,"End");
2785
2786 }
2787
2788 //____________________________________________
2789 void  AliAnaPi0EbE::MakeInvMassInCalorimeter()
2790 {
2791   //Do analysis and fill aods
2792   //Search for the photon decay in calorimeters
2793   //Read photon list from AOD, produced in class AliAnaPhoton
2794   //Check if 2 photons have the mass of the pi0.
2795
2796   if(!GetInputAODBranch())
2797   {
2798     AliFatal(Form("No input calo photons in AOD with name branch < %s >, STOP",GetInputAODName().Data()));
2799     return; // coverity
2800   }
2801
2802   //Get shower shape information of clusters
2803   TObjArray *clusters = 0;
2804   if     (GetCalorimeter()==kEMCAL) clusters = GetEMCALClusters();
2805   else if(GetCalorimeter()==kPHOS)  clusters = GetPHOSClusters() ;
2806
2807   Int_t nphoton = GetInputAODBranch()->GetEntriesFast();
2808   for(Int_t iphoton = 0; iphoton < nphoton-1; iphoton++)
2809   {
2810     AliAODPWG4Particle * photon1 =  (AliAODPWG4Particle*) (GetInputAODBranch()->At(iphoton));
2811
2812     // Vertex cut in case of mixed events
2813     Int_t evtIndex1 = 0 ;
2814     if(GetMixedEvent())
2815     {
2816       evtIndex1 = GetMixedEvent()->EventIndexForCaloCluster(photon1->GetCaloLabel(0)) ;
2817       if(TMath::Abs(GetVertex(evtIndex1)[2]) > GetZvertexCut()) continue ;  //vertex cut
2818     }
2819
2820     fMomentum1 = *(photon1->Momentum());
2821
2822     //Get original cluster, to recover some information
2823     Int_t iclus = -1;
2824     AliVCluster *cluster1 = FindCluster(clusters,photon1->GetCaloLabel(0),iclus);
2825
2826     if(!cluster1)
2827     {
2828       AliWarning("First cluster not found");
2829       return;
2830     }
2831
2832     for(Int_t jphoton = iphoton+1; jphoton < nphoton; jphoton++)
2833     {
2834       AliAODPWG4Particle * photon2 =  (AliAODPWG4Particle*) (GetInputAODBranch()->At(jphoton));
2835
2836       // Do analysis only when one of the decays is isolated
2837       // Run AliAnaParticleIsolation before
2838       if(fSelectIsolatedDecay)
2839       {
2840         Bool_t isolated1 = ((AliAODPWG4ParticleCorrelation*) photon1)->IsIsolated();
2841         Bool_t isolated2 = ((AliAODPWG4ParticleCorrelation*) photon2)->IsIsolated();
2842         if(!isolated1 && !isolated2) continue;
2843       }
2844
2845       // Vertex cut in case of mixed events
2846       Int_t evtIndex2 = 0 ;
2847       if(GetMixedEvent())
2848       {
2849         evtIndex2 = GetMixedEvent()->EventIndexForCaloCluster(photon2->GetCaloLabel(0)) ;
2850
2851         if(evtIndex1 == evtIndex2)
2852           continue ;
2853
2854         if(TMath::Abs(GetVertex(evtIndex2)[2]) > GetZvertexCut()) continue ;  //vertex cut
2855       }
2856
2857       fMomentum2 = *(photon2->Momentum());
2858
2859       //Get original cluster, to recover some information
2860       Int_t iclus2 = -1;
2861       AliVCluster *cluster2 = FindCluster(clusters,photon2->GetCaloLabel(0),iclus2,iclus+1);
2862       // start new loop from iclus1+1 to gain some time
2863
2864       if(!cluster2)
2865       {
2866         AliWarning("Second cluster not found");
2867         return;
2868       }
2869
2870       Float_t e1    = photon1->E();
2871       Float_t e2    = photon2->E();
2872
2873       //Select clusters with good time window difference
2874       Float_t tof1  = cluster1->GetTOF()*1e9;;
2875       Float_t tof2  = cluster2->GetTOF()*1e9;;
2876       Double_t t12diff = tof1-tof2;
2877       fhEPairDiffTime->Fill(e1+e2,    t12diff);
2878       if(TMath::Abs(t12diff) > GetPairTimeCut()) continue;
2879
2880       //Play with the MC stack if available
2881       Int_t mcIndex = kmcHadron;
2882       Int_t tag     = 0;
2883       Int_t label   =-1;
2884       if(IsDataMC())
2885       {
2886         HasPairSameMCMother(photon1->GetLabel(), photon2->GetLabel(),
2887                             photon1->GetTag()  , photon2->GetTag(),
2888                             label, tag) ;
2889         mcIndex = GetMCIndex(tag);
2890       }
2891
2892       // Check the invariant mass for different selection on the local maxima
2893       Int_t nMaxima1 = photon1->GetNLM();
2894       Int_t nMaxima2 = photon2->GetNLM();
2895
2896       fMomentum = fMomentum1+fMomentum2;
2897
2898       Double_t mass  = fMomentum.M();
2899       Double_t epair = fMomentum.E();
2900       Float_t ptpair = fMomentum.Pt();
2901
2902       if(fFillAllNLMHistograms)
2903       {
2904         if(nMaxima1==nMaxima2)
2905         {
2906           if     (nMaxima1==1) fhMassPairLocMax[0]->Fill(epair,mass);
2907           else if(nMaxima1==2) fhMassPairLocMax[1]->Fill(epair,mass);
2908           else                 fhMassPairLocMax[2]->Fill(epair,mass);
2909         }
2910         else if(nMaxima1==1 || nMaxima2==1)
2911         {
2912           if  (nMaxima1==2 || nMaxima2==2) fhMassPairLocMax[3]->Fill(epair,mass);
2913           else                             fhMassPairLocMax[4]->Fill(epair,mass);
2914         }
2915         else
2916           fhMassPairLocMax[5]->Fill(epair,mass);
2917
2918         // combinations with SS axis cut and NLM cut
2919         if(nMaxima1 == 1 && cluster2->GetM02() > 0.3) fhMassPairLocMax[6]->Fill(epair,mass);
2920         if(nMaxima2 == 1 && cluster1->GetM02() > 0.3) fhMassPairLocMax[6]->Fill(epair,mass);
2921         if(nMaxima1 >  1 && cluster2->GetM02() < 0.3 && cluster2->GetM02()> 0.1 ) fhMassPairLocMax[7]->Fill(epair,mass);
2922         if(nMaxima2 >  1 && cluster1->GetM02() < 0.3 && cluster1->GetM02()> 0.1 ) fhMassPairLocMax[7]->Fill(epair,mass);
2923       }
2924
2925       //
2926       // Skip events with too few or too many  NLM
2927       //
2928       if((nMaxima1 < fNLMCutMin || nMaxima1 > fNLMCutMax) || (nMaxima2 < fNLMCutMin || nMaxima2 > fNLMCutMax))
2929       {
2930         AliDebug(1,Form("NLM out of range: cluster1 %d, cluster2 %d",nMaxima1, nMaxima2));
2931         continue ;
2932       }
2933
2934       //Mass of all pairs
2935       fhMass  ->Fill( epair,mass);
2936       fhMassPt->Fill(ptpair,mass);
2937       if(IsDataMC() && mcIndex < 2) fhMCMassPt[mcIndex]->Fill(ptpair,mass);
2938
2939       if(fSelectPairInIsoCone && fMomentum1.Pt() > fIsoCandMinPt)
2940       {
2941         //Double_t angleOp  = fMomentum1.Angle(fMomentum2.Vect());
2942         Double_t radius   = GetIsolationCut()->Radius(fMomentum1.Eta(),fMomentum1.Phi(),fMomentum2.Eta(),fMomentum2.Phi());
2943
2944         //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);
2945
2946         if(radius < fR) fhMassPtIsoRCut->Fill(ptpair,mass);
2947       }
2948
2949       //
2950       // Select good pair (good phi, pt cuts, aperture and invariant mass)
2951       //
2952       if(!GetNeutralMesonSelection()->SelectPair(fMomentum1, fMomentum2,GetCalorimeter())) continue;
2953
2954       AliDebug(1,Form("Selected gamma pair: pt %f, phi %f, eta%f",
2955                       fMomentum.Pt(), fMomentum.Phi()*TMath::RadToDeg(), fMomentum.Eta()));
2956
2957       //
2958       // Tag both photons as decay if not done before
2959       // set the corresponding bit for pi0 or eta or "side" case
2960       //
2961       Int_t bit1 = photon1->DecayTag();
2962       if( bit1 < 0 ) bit1 = 0 ;
2963       if( !GetNeutralMesonSelection()->CheckDecayBit(bit1) )
2964       {
2965         AliDebug(1,Form("pT1 %2.2f; bit requested %d; decay bit1: In %d",fMomentum1.Pt(), GetNeutralMesonSelection()->GetDecayBit(), bit1));
2966
2967         GetNeutralMesonSelection()->SetDecayBit(bit1);
2968         photon1->SetDecayTag(bit1);
2969
2970         AliDebug(1,Form("\t Out %d", bit1));
2971
2972         fhPtDecay->Fill(photon1->Pt());
2973
2974         //Fill some histograms about shower shape
2975         if(fFillSelectClHisto && cluster1 && GetReader()->GetDataType()!=AliCaloTrackReader::kMC)
2976           FillSelectedClusterHistograms(cluster1, fMomentum1.Pt(), nMaxima1, photon1->GetTag());
2977
2978         if(IsDataMC())
2979         {
2980           Int_t mcIndex1 = GetMCIndex(photon1->GetTag());
2981           fhMCPtDecay[mcIndex1]->Fill(photon1->Pt());
2982           if(GetMCAnalysisUtils()->CheckTagBit(photon1->GetTag(),AliMCAnalysisUtils::kMCDecayPairLost))
2983           {
2984             if     ( mcIndex1 == kmcPi0Decay ) fhMCPtDecayLostPairPi0->Fill(photon1->Pt());
2985             else if( mcIndex1 == kmcEtaDecay ) fhMCPtDecayLostPairEta->Fill(photon1->Pt());
2986           }
2987         }
2988       }
2989
2990       Int_t bit2 = photon2->DecayTag();
2991       if( bit2 < 0 ) bit2 = 0 ;
2992       if( !GetNeutralMesonSelection()->CheckDecayBit(bit2) )
2993       {
2994         AliDebug(1,Form("pT2 %2.2f; bit requested %d; decay bit2: In %d",fMomentum2.Pt(), GetNeutralMesonSelection()->GetDecayBit(), bit2));
2995
2996         GetNeutralMesonSelection()->SetDecayBit(bit2);
2997         photon2->SetDecayTag(bit2);
2998
2999         AliDebug(1,Form("\t Out %d", bit2));
3000
3001         fhPtDecay->Fill(photon2->Pt());
3002
3003         //Fill some histograms about shower shape
3004         if(fFillSelectClHisto && cluster2 && GetReader()->GetDataType()!=AliCaloTrackReader::kMC)
3005           FillSelectedClusterHistograms(cluster2, fMomentum2.Pt(), nMaxima2, photon2->GetTag());
3006
3007         if(IsDataMC())
3008         {
3009           Int_t mcIndex2 = GetMCIndex(photon2->GetTag());
3010           fhMCPtDecay[mcIndex2]->Fill(photon2->Pt());
3011           if(GetMCAnalysisUtils()->CheckTagBit(photon2->GetTag(),AliMCAnalysisUtils::kMCDecayPairLost))
3012           {
3013             if     ( mcIndex2 == kmcPi0Decay ) fhMCPtDecayLostPairPi0->Fill(photon2->Pt());
3014             else if( mcIndex2 == kmcEtaDecay ) fhMCPtDecayLostPairEta->Fill(photon2->Pt());
3015           }
3016         }
3017       }
3018
3019       //Mass of selected pairs
3020       fhSelectedMass  ->Fill( epair,mass);
3021       fhSelectedMassPt->Fill(ptpair,mass);
3022       if(IsDataMC()  && mcIndex < 2) fhMCSelectedMassPt[mcIndex]->Fill(ptpair,mass);
3023
3024       // Fill histograms to undertand pile-up before other cuts applied
3025       // Remember to relax time cuts in the reader
3026       FillPileUpHistograms(ptpair,((cluster1->GetTOF()+cluster2->GetTOF())*1e9)/2,cluster1);
3027
3028       //Create AOD for analysis
3029       AliAODPWG4Particle pi0 = AliAODPWG4Particle(fMomentum);
3030
3031       if     ( (GetNeutralMesonSelection()->GetParticle()).Contains("Pi0") ) pi0.SetIdentifiedParticleType(AliCaloPID::kPi0);
3032       else if( (GetNeutralMesonSelection()->GetParticle()).Contains("Eta") ) pi0.SetIdentifiedParticleType(AliCaloPID::kEta);
3033       else
3034       {
3035         AliWarning("Particle type declared in AliNeutralMeson not correct, do not add");
3036         return ;
3037       }
3038       pi0.SetDetectorTag(photon1->GetDetectorTag());
3039
3040       // MC
3041       pi0.SetLabel(label);
3042       pi0.SetTag(tag);
3043
3044       //Set the indeces of the original caloclusters
3045       pi0.SetCaloLabel(photon1->GetCaloLabel(0), photon2->GetCaloLabel(0));
3046       //pi0.SetInputFileIndex(input);
3047
3048       AddAODParticle(pi0);
3049
3050     }//2n photon loop
3051
3052   }//1st photon loop
3053
3054   AliDebug(1,"End fill AODs");
3055
3056 }
3057
3058 //__________________________________________________
3059 void  AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS()
3060 {
3061   //Do analysis and fill aods
3062   //Search for the photon decay in calorimeters
3063   //Read photon list from AOD, produced in class AliAnaPhoton and AliGammaConversion
3064   //Check if 2 photons have the mass of the pi0.
3065
3066   // Check calorimeter input
3067   if(!GetInputAODBranch())
3068   {
3069     AliFatal(Form("No input calo photons in AOD branch with name < %s > , STOP",GetInputAODName().Data()));
3070   }
3071
3072   // Get the array with conversion photons
3073   TClonesArray * inputAODGammaConv = (TClonesArray *) GetReader()->GetOutputEvent()->FindListObject(fInputAODGammaConvName);
3074   if(!inputAODGammaConv)
3075   {
3076     inputAODGammaConv = (TClonesArray *) GetReader()->GetInputEvent()->FindListObject(fInputAODGammaConvName);
3077
3078     if(!inputAODGammaConv)
3079     {
3080       AliFatal(Form("No input gamma conversions in AOD branch with name < %s >",fInputAODGammaConvName.Data()));
3081       return; // coverity
3082     }
3083   }
3084
3085   //Get shower shape information of clusters
3086   TObjArray *clusters = 0;
3087   if     (GetCalorimeter()==kEMCAL) clusters = GetEMCALClusters();
3088   else if(GetCalorimeter()==kPHOS)  clusters = GetPHOSClusters() ;
3089
3090   Int_t nCTS  = inputAODGammaConv->GetEntriesFast();
3091   Int_t nCalo = GetInputAODBranch()->GetEntriesFast();
3092   if(nCTS<=0 || nCalo <=0)
3093   {
3094     AliDebug(1,Form("nCalo %d, nCTS %d, cannot loop",nCalo,nCTS));
3095     return;
3096   }
3097
3098   AliDebug(1,Form("Number of conversion photons %d and number of clusters %d",nCTS,nCalo));
3099
3100   // Do the loop, first calo, second CTS
3101   for(Int_t iphoton = 0; iphoton < GetInputAODBranch()->GetEntriesFast(); iphoton++)
3102   {
3103     AliAODPWG4Particle * photon1 =  (AliAODPWG4Particle*) (GetInputAODBranch()->At(iphoton));
3104     fMomentum1 = *(photon1->Momentum());
3105
3106     // Do analysis only when one of the decays is isolated
3107     // Run AliAnaParticleIsolation before
3108     if(fSelectIsolatedDecay)
3109     {
3110       Bool_t isolated1 = ((AliAODPWG4ParticleCorrelation*) photon1)->IsIsolated();
3111       if(!isolated1) continue;
3112     }
3113
3114     //Get original cluster, to recover some information
3115     Int_t iclus = -1;
3116     AliVCluster *cluster = FindCluster(clusters,photon1->GetCaloLabel(0),iclus);
3117
3118     for(Int_t jphoton = 0; jphoton < nCTS; jphoton++)
3119     {
3120       AliAODPWG4Particle * photon2 =  (AliAODPWG4Particle*) (inputAODGammaConv->At(jphoton));
3121
3122       Int_t evtIndex = 0;
3123       if(GetMixedEvent())
3124       {
3125         evtIndex = GetMixedEvent()->EventIndexForCaloCluster(photon2->GetCaloLabel(0)) ;
3126         if(TMath::Abs(GetVertex(evtIndex)[2]) > GetZvertexCut()) continue ;  //vertex cut
3127       }
3128
3129       fMomentum2 = *(photon2->Momentum());
3130
3131       fMomentum = fMomentum1+fMomentum2;
3132
3133       Double_t mass  = fMomentum.M();
3134       Double_t epair = fMomentum.E();
3135       Float_t ptpair = fMomentum.Pt();
3136
3137       Int_t nMaxima = photon1->GetNLM();
3138       if(fFillAllNLMHistograms)
3139       {
3140         if     (nMaxima==1) fhMassPairLocMax[0]->Fill(epair,mass);
3141         else if(nMaxima==2) fhMassPairLocMax[1]->Fill(epair,mass);
3142         else                fhMassPairLocMax[2]->Fill(epair,mass);
3143       }
3144
3145       if(nMaxima < fNLMCutMin || nMaxima > fNLMCutMax)
3146       {
3147         AliDebug(1,Form("NLM %d out of range",nMaxima));
3148         continue ;
3149       }
3150
3151       //Play with the MC stack if available
3152       Int_t mcIndex = kmcHadron;
3153       Int_t tag     = 0;
3154       Int_t label   =-1;
3155       if(IsDataMC())
3156       {
3157         Int_t   label2 = photon2->GetLabel();
3158         if(label2 >= 0 )photon2->SetTag(GetMCAnalysisUtils()->CheckOrigin(label2, GetReader(),kCTS));
3159
3160         HasPairSameMCMother(photon1->GetLabel(), photon2->GetLabel(),
3161                             photon1->GetTag()  , photon2->GetTag(),
3162                             label, tag) ;
3163         mcIndex = GetMCIndex(tag);
3164       }
3165
3166       //Mass of selected pairs
3167       fhMass  ->Fill( epair,mass);
3168       fhMassPt->Fill(ptpair,mass);
3169       if(IsDataMC() && mcIndex < 2 ) fhMCMassPt[mcIndex]->Fill(ptpair,mass);
3170
3171       //
3172       // Select good pair (good phi, pt cuts, aperture and invariant mass)
3173       //
3174       if(!GetNeutralMesonSelection()->SelectPair(fMomentum1, fMomentum2,GetCalorimeter())) continue ;
3175
3176       AliDebug(1,Form("Selected gamma pair: pt %f, phi %f, eta%f",fMomentum.Pt(), fMomentum.Phi()*TMath::RadToDeg(), fMomentum.Eta()));
3177
3178       //
3179       // Tag both photons as decay if not done before
3180       // set the corresponding bit for pi0 or eta or "side" case
3181       //
3182       Int_t bit1 = photon1->DecayTag();
3183       if( bit1 < 0 ) bit1 = 0 ;
3184       if( !GetNeutralMesonSelection()->CheckDecayBit(bit1) )
3185       {
3186         GetNeutralMesonSelection()->SetDecayBit(bit1);
3187         photon1->SetDecayTag(bit1);
3188
3189         fhPtDecay->Fill(photon1->Pt());
3190
3191         //Fill some histograms about shower shape
3192         if(fFillSelectClHisto && cluster && GetReader()->GetDataType()!=AliCaloTrackReader::kMC)
3193           FillSelectedClusterHistograms(cluster, fMomentum1.Pt(), nMaxima, photon1->GetTag());
3194
3195         if(IsDataMC())
3196         {
3197           Int_t mcIndex1 = GetMCIndex(photon1->GetTag());
3198           fhMCPtDecay[mcIndex1]->Fill(photon1->Pt());
3199           if(GetMCAnalysisUtils()->CheckTagBit(photon1->GetTag(),AliMCAnalysisUtils::kMCDecayPairLost))
3200           {
3201             if     ( mcIndex1 == kmcPi0Decay ) fhMCPtDecayLostPairPi0->Fill(photon1->Pt());
3202             else if( mcIndex1 == kmcEtaDecay ) fhMCPtDecayLostPairEta->Fill(photon1->Pt());
3203           }
3204         }
3205       }
3206
3207       Int_t bit2 = photon2->DecayTag();
3208       if( bit2 < 0 ) bit2 = 0 ;
3209       if( !GetNeutralMesonSelection()->CheckDecayBit(bit2) )
3210       {
3211         GetNeutralMesonSelection()->SetDecayBit(bit2);
3212         photon2->SetDecayTag(bit2);
3213       }
3214
3215       //Mass of selected pairs
3216       fhSelectedMass  ->Fill( epair,mass);
3217       fhSelectedMassPt->Fill(ptpair,mass);
3218       if(IsDataMC()  && mcIndex < 2) fhMCSelectedMassPt[mcIndex]->Fill(ptpair,mass);
3219
3220       // Fill histograms to undertand pile-up before other cuts applied
3221       // Remember to relax time cuts in the reader
3222       if(cluster) FillPileUpHistograms(fMomentum.Pt(),cluster->GetTOF()*1e9,cluster);
3223
3224       //Create AOD for analysis
3225
3226       AliAODPWG4Particle pi0 = AliAODPWG4Particle(fMomentum);
3227
3228       if     ( (GetNeutralMesonSelection()->GetParticle()).Contains("Pi0") ) pi0.SetIdentifiedParticleType(AliCaloPID::kPi0);
3229       else if( (GetNeutralMesonSelection()->GetParticle()).Contains("Eta") ) pi0.SetIdentifiedParticleType(AliCaloPID::kEta);
3230       else
3231       {
3232         AliWarning("Particle type declared in AliNeutralMeson not correct, do not add");
3233         return ;
3234       }
3235       pi0.SetDetectorTag(photon1->GetDetectorTag());
3236
3237       // MC
3238       pi0.SetLabel(label);
3239       pi0.SetTag(tag);
3240
3241       //Set the indeces of the original tracks or caloclusters
3242       pi0.SetCaloLabel (photon1->GetCaloLabel(0) , -1);
3243       pi0.SetTrackLabel(photon2->GetTrackLabel(0), photon2->GetTrackLabel(1));
3244       //pi0.SetInputFileIndex(input);
3245
3246       AddAODParticle(pi0);
3247
3248     }//2n photon loop
3249
3250   }//1st photon loop
3251
3252   AliDebug(1,"End fill AODs");
3253
3254 }
3255
3256 //_________________________________________________
3257 void  AliAnaPi0EbE::MakeShowerShapeIdentification()
3258 {
3259   //Search for pi0 in GetCalorimeter() with shower shape analysis
3260
3261   TObjArray * pl        = 0x0;
3262   AliVCaloCells * cells = 0x0;
3263   //Select the Calorimeter of the photon
3264   if      (GetCalorimeter() == kEMCAL )
3265   {
3266     pl    = GetEMCALClusters();
3267     cells = GetEMCALCells();
3268   }
3269   else if (GetCalorimeter() == kPHOS)
3270   {
3271     AliFatal("kSSCalo case not implememted for PHOS");
3272     return; // for coverity
3273
3274     //pl    = GetPHOSClusters();
3275     //cells = GetPHOSCells();
3276   }
3277
3278   if(!pl)
3279   {
3280     AliInfo(Form("TObjArray with %s clusters is NULL!",GetCalorimeterString().Data()));
3281     return;
3282   }
3283
3284   for(Int_t icalo = 0; icalo < pl->GetEntriesFast(); icalo++)
3285   {
3286     AliVCluster * calo = (AliVCluster*) (pl->At(icalo));
3287
3288     Int_t evtIndex = 0 ;
3289     if (GetMixedEvent())
3290     {
3291       evtIndex=GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ;
3292       if(TMath::Abs(GetVertex(evtIndex)[2]) > GetZvertexCut()) continue ;  //vertex cut
3293     }
3294
3295     //Get Momentum vector,
3296     Double_t vertex[]={0,0,0};
3297     if(GetReader()->GetDataType() != AliCaloTrackReader::kMC)
3298     {
3299       calo->GetMomentum(fMomentum,GetVertex(evtIndex)) ;
3300     }//Assume that come from vertex in straight line
3301     else
3302     {
3303       calo->GetMomentum(fMomentum,vertex) ;
3304     }
3305
3306     //If too small or big pt, skip it
3307     if(fMomentum.E() < GetMinEnergy() || fMomentum.E() > GetMaxEnergy() ) continue ;
3308
3309     //Check acceptance selection
3310     if(IsFiducialCutOn())
3311     {
3312       Bool_t in = GetFiducialCut()->IsInFiducialCut(fMomentum.Eta(),fMomentum.Phi(),GetCalorimeter()) ;
3313       if(! in ) continue ;
3314     }
3315
3316     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()));
3317
3318     //Play with the MC stack if available
3319     //Check origin of the candidates
3320     Int_t tag   = 0 ;
3321     if(IsDataMC())
3322     {
3323       tag = GetMCAnalysisUtils()->CheckOrigin(calo->GetLabels(),calo->GetNLabels(),GetReader(),GetCalorimeter());
3324       //GetMCAnalysisUtils()->CheckMultipleOrigin(calo->GetLabels(),calo->GetNLabels(), GetReader(), aodpi0.GetInputFileIndex(), tag);
3325       AliDebug(1,Form("Origin of candidate %d",tag));
3326     }
3327
3328     //Int_t nMaxima = GetCaloUtils()->GetNumberOfLocalMaxima(calo, cells); // NLM
3329
3330     //Check Distance to Bad channel, set bit.
3331     Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
3332     if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;
3333     if(distBad < fMinDist){ //In bad channel (PHOS cristal size 2.2x2.2 cm)
3334       //FillRejectedClusterHistograms(tag,nMaxima);
3335       continue ;
3336     }
3337
3338     AliDebug(1,Form("Bad channel cut passed %4.2f",distBad));
3339
3340     //If too low number of cells, skip it
3341     if ( calo->GetNCells() < GetCaloPID()->GetClusterSplittingMinNCells())
3342     {
3343       //FillRejectedClusterHistograms(tag,nMaxima);
3344       continue ;
3345     }
3346
3347     AliDebug(1,Form("N cells cut passed %d > %d",calo->GetNCells(), GetCaloPID()->GetClusterSplittingMinNCells()));
3348
3349     //.......................................
3350     // TOF cut, BE CAREFUL WITH THIS CUT
3351     Double_t tof = calo->GetTOF()*1e9;
3352     if(tof < fTimeCutMin || tof > fTimeCutMax)
3353     {
3354       //FillRejectedClusterHistograms(tag,nMaxima);
3355       continue ;
3356     }
3357
3358     //Check PID
3359     //PID selection or bit setting
3360     Int_t    nMaxima  = 0;
3361     Double_t mass     = 0, angle    = 0;
3362     Int_t    absId1   =-1, absId2   =-1;
3363     Float_t  distbad1 =-1, distbad2 =-1;
3364     Bool_t   fidcut1  = 0, fidcut2  = 0;
3365
3366     Int_t idPartType = GetCaloPID()->GetIdentifiedParticleTypeFromClusterSplitting(calo,cells,GetCaloUtils(),
3367                                                                                    GetVertex(evtIndex),nMaxima,
3368                                                                                    mass,angle,
3369                                                                                    fMomentum1,fMomentum2,
3370                                                                                    absId1,absId2,
3371                                                                                    distbad1,distbad2,
3372                                                                                    fidcut1,fidcut2) ;
3373
3374
3375     AliDebug(1,Form("PDG of identified particle %d",idPartType));
3376
3377     // Skip events where one of the new clusters (lowest energy) is close to an EMCal border or a bad channel
3378     if( (fCheckSplitDistToBad) &&
3379        (!fidcut2 || !fidcut1 || distbad1 < fMinDist || distbad2 < fMinDist))
3380     {
3381       AliDebug(1,Form("Dist to bad channel cl %f, cl1 %f, cl2 %f; fid cl1 %d, cl2 %d",
3382                       calo->GetDistanceToBadChannel(),distbad1,distbad2, fidcut1,fidcut2));
3383
3384       //FillRejectedClusterHistograms(tag,nMaxima);
3385       continue ;
3386     }
3387
3388     //Skip events with too few or too many  NLM
3389     if(nMaxima < fNLMCutMin || nMaxima > fNLMCutMax)
3390     {
3391       //FillRejectedClusterHistograms(tag,nMaxima);
3392       continue ;
3393     }
3394
3395     AliDebug(1,Form("NLM %d accepted",nMaxima));
3396
3397     //Skip matched clusters with tracks
3398     if(fRejectTrackMatch && IsTrackMatched(calo, GetReader()->GetInputEvent()))
3399     {
3400       FillRejectedClusterHistograms(tag,nMaxima);
3401       continue ;
3402     }
3403
3404     Float_t l0 = calo->GetM02();
3405     Float_t e1 = fMomentum1.Energy();
3406     Float_t e2 = fMomentum2.Energy();
3407     fMomentum12 = fMomentum1+fMomentum2;
3408     Float_t ptSplit = fMomentum12.Pt();
3409     Float_t  eSplit = e1+e2;
3410
3411     //mass of all clusters
3412     fhMass       ->Fill(fMomentum.E() ,mass);
3413     fhMassPt     ->Fill(fMomentum.Pt(),mass);
3414     fhMassSplitPt->Fill(ptSplit ,mass);
3415     fhPtLambda0NoSplitCut->Fill(fMomentum.Pt(),l0);
3416
3417     // Asymmetry of all clusters
3418     Float_t asy =-10;
3419
3420     if(e1+e2 > 0) asy = (e1-e2) / (e1+e2);
3421     fhAsymmetry->Fill(fMomentum.E(),asy);
3422
3423     // Divide NLM in 3 cases, 1 local maxima, 2 local maxima, more than 2 local maxima
3424     Int_t indexMax = -1;
3425     if     (nMaxima==1) indexMax = 0 ;
3426     else if(nMaxima==2) indexMax = 1 ;
3427     else                indexMax = 2 ;
3428     fhMassPtLocMax[indexMax]->Fill(fMomentum.Pt(),mass);
3429
3430     Int_t   mcIndex   =-1;
3431     Int_t   noverlaps = 0;
3432     Float_t ptprim    = 0;
3433     if(IsDataMC())
3434     {
3435       mcIndex = GetMCIndex(tag);
3436
3437       Bool_t ok      = kFALSE;
3438       Int_t  mcLabel = calo->GetLabel();
3439
3440       fPrimaryMom = GetMCAnalysisUtils()->GetMother(mcLabel,GetReader(),ok);
3441
3442       Int_t mesonLabel = -1;
3443
3444       if(mcIndex == kmcPi0 || mcIndex == kmcEta)
3445       {
3446         if(mcIndex == kmcPi0)
3447         {
3448           fGrandMotherMom = GetMCAnalysisUtils()->GetMotherWithPDG(mcLabel,111,GetReader(),ok,mesonLabel);
3449           if(fGrandMotherMom.E() > 0 && ok) ptprim =  fGrandMotherMom.Pt();
3450         }
3451         else
3452         {
3453           fGrandMotherMom = GetMCAnalysisUtils()->GetMotherWithPDG(mcLabel,221,GetReader(),ok,mesonLabel);
3454           if(fGrandMotherMom.E() > 0 && ok) ptprim = fGrandMotherMom.Pt();
3455         }
3456       }
3457
3458       const UInt_t nlabels = calo->GetNLabels();
3459       Int_t overpdg[nlabels];
3460       noverlaps = GetMCAnalysisUtils()->GetNOverlaps(calo->GetLabels(), nlabels,tag,mesonLabel,GetReader(),overpdg);
3461
3462       fhMCMassPt     [mcIndex]->Fill(fMomentum.Pt(),mass);
3463       fhMCMassSplitPt[mcIndex]->Fill(ptSplit ,mass);
3464       if(mcIndex==kmcPi0)
3465       {
3466         fhMCPi0PtRecoPtPrim                     ->Fill(fMomentum.Pt(),ptprim);
3467         fhMCPi0SplitPtRecoPtPrim                ->Fill(ptSplit ,ptprim);
3468         fhMCPi0PtRecoPtPrimLocMax     [indexMax]->Fill(fMomentum.Pt(),ptprim);
3469         fhMCPi0SplitPtRecoPtPrimLocMax[indexMax]->Fill(ptSplit ,ptprim);
3470
3471       }
3472       else if(mcIndex==kmcEta)
3473       {
3474         fhMCEtaPtRecoPtPrim                     ->Fill(fMomentum.Pt(),ptprim);
3475         fhMCEtaSplitPtRecoPtPrim                ->Fill(ptSplit ,ptprim);
3476         fhMCEtaPtRecoPtPrimLocMax     [indexMax]->Fill(fMomentum.Pt(),ptprim);
3477         fhMCEtaSplitPtRecoPtPrimLocMax[indexMax]->Fill(ptSplit ,ptprim);
3478       }
3479
3480       if(noverlaps==0)
3481       {
3482         if(mcIndex==kmcPi0)
3483         {
3484           fhMCPi0PtRecoPtPrimNoOverlap     ->Fill(fMomentum.Pt(),ptprim);
3485           fhMCPi0SplitPtRecoPtPrimNoOverlap->Fill(ptSplit ,ptprim);
3486         }
3487         else if(mcIndex==kmcEta)
3488         {
3489           fhMCEtaPtRecoPtPrimNoOverlap     ->Fill(fMomentum.Pt(),ptprim);
3490           fhMCEtaSplitPtRecoPtPrimNoOverlap->Fill(ptSplit ,ptprim);
3491         }
3492
3493         fhMassNoOverlap       ->Fill(fMomentum.E() ,mass);
3494         fhMassPtNoOverlap     ->Fill(fMomentum.Pt(),mass);
3495         fhMassSplitPtNoOverlap->Fill(ptSplit ,mass);
3496
3497         fhMCMassPtNoOverlap     [mcIndex]->Fill(fMomentum.Pt(),mass);
3498         fhMCMassSplitPtNoOverlap[mcIndex]->Fill(ptSplit ,mass);
3499       }
3500
3501       fhMCPtAsymmetry[mcIndex]->Fill(fMomentum.Pt(),asy);
3502     }
3503
3504     // If cluster does not pass pid, not pi0/eta, skip it.
3505     if     (GetOutputAODName().Contains("Pi0") && idPartType != AliCaloPID::kPi0)
3506     {
3507       AliDebug(1,"Cluster is not Pi0");
3508       FillRejectedClusterHistograms(tag,nMaxima);
3509       continue ;
3510     }
3511
3512     else if(GetOutputAODName().Contains("Eta") && idPartType != AliCaloPID::kEta)
3513     {
3514       AliDebug(1,"Cluster is not Eta");
3515       FillRejectedClusterHistograms(tag,nMaxima);
3516       continue ;
3517     }
3518
3519     AliDebug(1,Form("Pi0/Eta selection cuts passed: pT %3.2f, pdg %d",fMomentum.Pt(), idPartType));
3520
3521     //Mass and asymmetry of selected pairs
3522     fhSelectedAsymmetry  ->Fill(fMomentum.E() ,asy );
3523     fhSelectedMass       ->Fill(fMomentum.E() ,mass);
3524     fhSelectedMassPt     ->Fill(fMomentum.Pt(),mass);
3525     fhSelectedMassSplitPt->Fill(ptSplit ,mass);
3526     fhSelectedMassPtLocMax[indexMax]->Fill(fMomentum.Pt(),mass);
3527
3528     Int_t   nSM  = GetModuleNumber(calo);
3529     if(nSM < GetCaloUtils()->GetNumberOfSuperModulesUsed() && nSM >=0 && fFillAllNLMHistograms)
3530     {
3531       fhSelectedMassPtLocMaxSM   [indexMax][nSM]->Fill(fMomentum.Pt(),mass);
3532       fhSelectedLambda0PtLocMaxSM[indexMax][nSM]->Fill(fMomentum.Pt(),l0  );
3533     }
3534
3535     if(IsDataMC())
3536     {
3537       if(mcIndex==kmcPi0)
3538       {
3539         fhMCPi0SelectedPtRecoPtPrim                     ->Fill(fMomentum.Pt(),ptprim);
3540         fhMCPi0SelectedSplitPtRecoPtPrim                ->Fill(ptSplit ,ptprim);
3541         fhMCPi0SelectedPtRecoPtPrimLocMax     [indexMax]->Fill(fMomentum.Pt(),ptprim);
3542         fhMCPi0SelectedSplitPtRecoPtPrimLocMax[indexMax]->Fill(ptSplit ,ptprim);
3543       }
3544       else if(mcIndex==kmcEta)
3545       {
3546         fhMCEtaSelectedPtRecoPtPrim                     ->Fill(fMomentum.Pt(),ptprim);
3547         fhMCEtaSelectedSplitPtRecoPtPrim                ->Fill(ptSplit ,ptprim);
3548         fhMCEtaSelectedPtRecoPtPrimLocMax     [indexMax]->Fill(fMomentum.Pt(),ptprim);
3549         fhMCEtaSelectedSplitPtRecoPtPrimLocMax[indexMax]->Fill(ptSplit ,ptprim);
3550       }
3551
3552       if(noverlaps==0)
3553       {
3554         fhSelectedMassNoOverlap       ->Fill(fMomentum.E() ,mass);
3555         fhSelectedMassPtNoOverlap     ->Fill(fMomentum.Pt(),mass);
3556         fhSelectedMassSplitPtNoOverlap->Fill(ptSplit ,mass);
3557
3558         if(mcIndex==kmcPi0)
3559         {
3560           fhMCPi0SelectedPtRecoPtPrimNoOverlap     ->Fill(fMomentum.Pt(),ptprim);
3561           fhMCPi0SelectedSplitPtRecoPtPrimNoOverlap->Fill(ptSplit ,ptprim);
3562         }
3563         else if(mcIndex==kmcEta)
3564         {
3565           fhMCEtaSelectedPtRecoPtPrimNoOverlap     ->Fill(fMomentum.Pt(),ptprim);
3566           fhMCEtaSelectedSplitPtRecoPtPrimNoOverlap->Fill(ptSplit ,ptprim);
3567         }
3568       }
3569     }
3570
3571     fhSplitE        ->Fill( eSplit);
3572     fhSplitPt       ->Fill(ptSplit);
3573     Float_t phi = fMomentum.Phi();
3574     if(phi<0) phi+=TMath::TwoPi();
3575     fhSplitPtPhi    ->Fill(ptSplit,phi);
3576     fhSplitPtEta    ->Fill(ptSplit,fMomentum.Eta());
3577     fhNLocMaxSplitPt->Fill(ptSplit ,nMaxima);
3578
3579     //Check split-clusters with good time window difference
3580     Double_t tof1  = cells->GetCellTime(absId1);
3581     GetCaloUtils()->RecalibrateCellTime(tof1, GetCalorimeter(), absId1,GetReader()->GetInputEvent()->GetBunchCrossNumber());
3582     tof1*=1.e9;
3583
3584     Double_t tof2  = cells->GetCellTime(absId2);
3585     GetCaloUtils()->RecalibrateCellTime(tof2, GetCalorimeter(), absId2,GetReader()->GetInputEvent()->GetBunchCrossNumber());
3586     tof2*=1.e9;
3587
3588     Double_t t12diff = tof1-tof2;
3589     fhEPairDiffTime->Fill(e1+e2,    t12diff);
3590
3591     if(IsDataMC())
3592     {
3593       fhMCSplitE        [mcIndex]->Fill( eSplit);
3594       fhMCSplitPt       [mcIndex]->Fill(ptSplit);
3595       fhMCSplitPtPhi    [mcIndex]->Fill(ptSplit,phi);
3596       fhMCSplitPtEta    [mcIndex]->Fill(ptSplit,fMomentum.Eta());
3597       fhMCNLocMaxSplitPt[mcIndex]->Fill(ptSplit ,nMaxima);
3598       fhMCNLocMaxPt     [mcIndex]->Fill(fMomentum.Pt(),nMaxima);
3599
3600       fhMCSelectedMassPt     [mcIndex]->Fill(fMomentum.Pt(),mass);
3601       fhMCSelectedMassSplitPt[mcIndex]->Fill(ptSplit,mass);
3602       fhMCSelectedMassPtLocMax[mcIndex][indexMax]->Fill(fMomentum.Pt(),mass);
3603
3604       if(noverlaps==0)
3605       {
3606         fhMCSelectedMassPtNoOverlap     [mcIndex]->Fill(fMomentum.Pt(),mass);
3607         fhMCSelectedMassSplitPtNoOverlap[mcIndex]->Fill(ptSplit,mass);
3608       }
3609     }
3610
3611     // Remove clusters with NLM=x depeding on a minimim energy cut
3612     if(nMaxima == 1 && fNLMECutMin[0] > fMomentum.E()) continue;
3613     if(nMaxima == 2 && fNLMECutMin[1] > fMomentum.E()) continue;
3614     if(nMaxima >  2 && fNLMECutMin[2] > fMomentum.E()) continue;
3615
3616     //Fill some histograms about shower shape
3617     if(fFillSelectClHisto && GetReader()->GetDataType()!=AliCaloTrackReader::kMC)
3618     {
3619       FillSelectedClusterHistograms(calo, fMomentum.Pt(), nMaxima, tag, asy);
3620     }
3621
3622     // Fill histograms to undertand pile-up before other cuts applied
3623     // Remember to relax time cuts in the reader
3624     Double_t tofcluster   = calo->GetTOF()*1e9;
3625
3626     FillPileUpHistograms(fMomentum.Pt(),tofcluster,calo);
3627
3628     if(fFillEMCALBCHistograms && GetCalorimeter()==kEMCAL)
3629       FillEMCALBCHistograms(fMomentum.E(), fMomentum.Eta(), fMomentum.Phi(), tofcluster);
3630
3631     //-----------------------
3632     //Create AOD for analysis
3633
3634     AliAODPWG4Particle aodpi0 = AliAODPWG4Particle(fMomentum);
3635     aodpi0.SetLabel(calo->GetLabel());
3636
3637     //Set the indeces of the original caloclusters
3638     aodpi0.SetCaloLabel(calo->GetID(),-1);
3639     aodpi0.SetDetectorTag(GetCalorimeter());
3640
3641     if     (distBad > fMinDist3) aodpi0.SetDistToBad(2) ;
3642     else if(distBad > fMinDist2) aodpi0.SetDistToBad(1) ;
3643     else                         aodpi0.SetDistToBad(0) ;
3644
3645     // Check if cluster is pi0 via cluster splitting
3646     aodpi0.SetIdentifiedParticleType(idPartType);
3647
3648     aodpi0.SetM02(l0);
3649     aodpi0.SetNLM(nMaxima);
3650     aodpi0.SetTime(tofcluster);
3651     aodpi0.SetNCells(calo->GetNCells());
3652     aodpi0.SetSModNumber(nSM);
3653
3654     aodpi0.SetTag(tag);
3655
3656     //Add AOD with pi0 object to aod branch
3657     AddAODParticle(aodpi0);
3658
3659   }//loop
3660
3661   AliDebug(1,"End fill AODs");
3662
3663 }
3664
3665 //______________________________________________
3666 void  AliAnaPi0EbE::MakeAnalysisFillHistograms()
3667 {
3668   //Do analysis and fill histograms
3669
3670   if(!GetOutputAODBranch())
3671   {
3672     AliFatal(Form("No output pi0 in AOD branch with name < %s >,STOP",GetOutputAODName().Data()));
3673     return;
3674   }
3675
3676   //Loop on stored AOD pi0
3677   Int_t naod = GetOutputAODBranch()->GetEntriesFast();
3678
3679   AliDebug(1,Form("AOD branch entries %d", naod));
3680
3681   Float_t cen = GetEventCentrality();
3682   Float_t ep  = GetEventPlaneAngle();
3683
3684   for(Int_t iaod = 0; iaod < naod ; iaod++)
3685   {
3686     AliAODPWG4Particle* pi0 =  (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
3687     Int_t pdg = pi0->GetIdentifiedParticleType();
3688
3689     if( ( pdg != AliCaloPID::kPi0 && pdg != AliCaloPID::kEta ) ) continue;
3690
3691     //Fill pi0 histograms
3692     Float_t ener  = pi0->E();
3693     Float_t pt    = pi0->Pt();
3694     Float_t phi   = pi0->Phi();
3695     if(phi < 0) phi+=TMath::TwoPi();
3696     Float_t eta = pi0->Eta();
3697
3698     fhPt     ->Fill(pt  );
3699     fhE      ->Fill(ener);
3700
3701     fhPtEta  ->Fill(pt  ,eta);
3702     fhPtPhi  ->Fill(pt  ,phi);
3703     fhEtaPhi ->Fill(eta ,phi);
3704
3705     if(IsHighMultiplicityAnalysisOn())
3706     {
3707       fhPtCentrality ->Fill(pt,cen) ;
3708       fhPtEventPlane ->Fill(pt,ep ) ;
3709     }
3710
3711     if(IsDataMC())
3712     {
3713       Int_t tag     = pi0->GetTag();
3714       Int_t label   = pi0->GetLabel();
3715       Int_t mcIndex = GetMCIndex(tag);
3716
3717       if(fAnaType!=kSSCalo && mcIndex > 1) continue;
3718
3719       fhMCE    [mcIndex] ->Fill(ener);
3720       fhMCPt   [mcIndex] ->Fill(pt);
3721       fhMCPtPhi[mcIndex] ->Fill(pt,phi);
3722       fhMCPtEta[mcIndex] ->Fill(pt,eta);
3723
3724       if(IsHighMultiplicityAnalysisOn()) fhMCPtCentrality[mcIndex]->Fill(pt,cen);
3725
3726       if((mcIndex==kmcPi0Decay || mcIndex==kmcEtaDecay ||
3727           mcIndex==kmcPi0      || mcIndex==kmcEta         ) &&
3728          fAnaType==kSSCalo)
3729       {
3730         Float_t efracMC   = 0;
3731         Int_t   momlabel  = -1;
3732         Bool_t  ok        = kFALSE;
3733
3734         fPrimaryMom = GetMCAnalysisUtils()->GetMother(label,GetReader(),ok);
3735         if(!ok) continue;
3736
3737         if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0))
3738         {
3739           fGrandMotherMom = GetMCAnalysisUtils()->GetMotherWithPDG(label,111,GetReader(),ok,momlabel);
3740           if(fGrandMotherMom.E() > 0 && ok)
3741           {
3742             efracMC =  fGrandMotherMom.E()/ener;
3743             fhMCPi0PtGenRecoFraction ->Fill(pt,efracMC);
3744           }
3745         }
3746         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay))
3747         {
3748           fhMCPi0DecayPt->Fill(pt);
3749           fGrandMotherMom = GetMCAnalysisUtils()->GetMotherWithPDG(label,111,GetReader(),ok,momlabel);
3750           if(fGrandMotherMom.E() > 0 && ok)
3751           {
3752             efracMC =  fPrimaryMom.E()/fGrandMotherMom.E();
3753             fhMCPi0DecayPtFraction ->Fill(pt,efracMC);
3754           }
3755         }
3756         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta))
3757         {
3758           fGrandMotherMom = GetMCAnalysisUtils()->GetMotherWithPDG(label,221,GetReader(),ok,momlabel);
3759           if(fGrandMotherMom.E() > 0 && ok)
3760           {
3761             efracMC =  fGrandMotherMom.E()/ener;
3762             fhMCEtaPtGenRecoFraction ->Fill(pt,efracMC);
3763           }
3764         }
3765         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay))
3766         {
3767           fhMCEtaDecayPt->Fill(pt);
3768           fGrandMotherMom = GetMCAnalysisUtils()->GetMotherWithPDG(label,221,GetReader(),ok,momlabel);
3769           if(fGrandMotherMom.E() > 0 && ok)
3770           {
3771             efracMC =  fPrimaryMom.E()/fGrandMotherMom.E();
3772             fhMCEtaDecayPtFraction ->Fill(pt,efracMC);
3773           }
3774         }
3775         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton))
3776         {
3777           fhMCOtherDecayPt->Fill(pt);
3778         }
3779
3780       }
3781
3782       if( mcIndex==kmcPi0 || mcIndex==kmcEta )
3783       {
3784         Float_t prodR     = -1;
3785         Int_t   momindex  = -1;
3786         Int_t   mompdg    = -1;
3787         Int_t   momstatus = -1;
3788
3789         if(GetReader()->ReadStack())
3790         {
3791           TParticle* ancestor = GetMCStack()->Particle(label);
3792           momindex  = ancestor->GetFirstMother();
3793           if(momindex < 0) return;
3794           TParticle* mother = GetMCStack()->Particle(momindex);
3795           mompdg    = TMath::Abs(mother->GetPdgCode());
3796           momstatus = mother->GetStatusCode();
3797           prodR = mother->R();
3798         }
3799         else
3800         {
3801           TClonesArray * mcparticles = GetReader()->GetAODMCParticles();
3802           AliAODMCParticle* ancestor = (AliAODMCParticle *) mcparticles->At(label);
3803           momindex  = ancestor->GetMother();
3804           if(momindex < 0) return;
3805           AliAODMCParticle* mother = (AliAODMCParticle *) mcparticles->At(momindex);
3806           mompdg    = TMath::Abs(mother->GetPdgCode());
3807           momstatus = mother->GetStatus();
3808           prodR = TMath::Sqrt(mother->Xv()*mother->Xv()+mother->Yv()*mother->Yv());
3809         }
3810
3811         if( mcIndex==kmcPi0 )
3812         {
3813           fhMCPi0ProdVertex->Fill(pt,prodR);
3814
3815           if     (momstatus  == 21) fhMCPi0PtOrigin->Fill(pt,0.5);//parton
3816           else if(mompdg     < 22 ) fhMCPi0PtOrigin->Fill(pt,1.5);//quark
3817           else if(mompdg     > 2100  && mompdg   < 2210) fhMCPi0PtOrigin->Fill(pt,2.5);// resonances
3818           else if(mompdg    == 221) fhMCPi0PtOrigin->Fill(pt,8.5);//eta
3819           else if(mompdg    == 331) fhMCPi0PtOrigin->Fill(pt,9.5);//eta prime
3820           else if(mompdg    == 213) fhMCPi0PtOrigin->Fill(pt,4.5);//rho
3821           else if(mompdg    == 223) fhMCPi0PtOrigin->Fill(pt,5.5);//omega
3822           else if(mompdg    >= 310   && mompdg    <= 323) fhMCPi0PtOrigin->Fill(pt,6.5);//k0S, k+-,k*
3823           else if(mompdg    == 130) fhMCPi0PtOrigin->Fill(pt,6.5);//k0L
3824           else if(momstatus == 11 || momstatus  == 12 ) fhMCPi0PtOrigin->Fill(pt,3.5);//resonances
3825           else                      fhMCPi0PtOrigin->Fill(pt,7.5);//other?
3826         }
3827         else if (mcIndex==kmcEta )
3828         {
3829           fhMCEtaProdVertex->Fill(pt,prodR);
3830
3831           if     (momstatus  == 21) fhMCEtaPtOrigin->Fill(pt,0.5);//parton
3832           else if(mompdg     < 22 ) fhMCEtaPtOrigin->Fill(pt,1.5);//quark
3833           else if(mompdg     > 2100  && mompdg   < 2210) fhMCEtaPtOrigin->Fill(pt,2.5);// resonances
3834           else if(mompdg    == 221) fhMCEtaPtOrigin->Fill(pt,8.5);//eta
3835           else if(mompdg    == 331) fhMCEtaPtOrigin->Fill(pt,9.5);//eta prime
3836           else if(mompdg    == 213) fhMCEtaPtOrigin->Fill(pt,4.5);//rho
3837           else if(mompdg    == 223) fhMCEtaPtOrigin->Fill(pt,5.5);//omega
3838           else if(mompdg    >= 310   && mompdg    <= 323) fhMCEtaPtOrigin->Fill(pt,6.5);//k0S, k+-,k*
3839           else if(mompdg    == 130) fhMCEtaPtOrigin->Fill(pt,6.5);//k0L
3840           else if(momstatus == 11 || momstatus  == 12 ) fhMCEtaPtOrigin->Fill(pt,3.5);//resonances
3841           else                      fhMCEtaPtOrigin->Fill(pt,7.5);//other?
3842         }
3843       }
3844
3845     }//Histograms with MC
3846
3847   }// aod loop
3848
3849   AliDebug(1,"End");
3850
3851 }
3852
3853 //__________________________________________________________________
3854 void AliAnaPi0EbE::Print(const Option_t * opt) const
3855 {
3856   //Print some relevant parameters set for the analysis
3857   if(! opt)
3858     return;
3859
3860   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
3861   AliAnaCaloTrackCorrBaseClass::Print("");
3862   printf("Analysis Type = %d \n",  fAnaType) ;
3863   if(fAnaType == kSSCalo)
3864   {
3865     printf("Calorimeter            =     %s\n", GetCalorimeterString().Data()) ;
3866     printf("Min Distance to Bad Channel   = %2.1f\n",fMinDist);
3867     printf("Min Distance to Bad Channel 2 = %2.1f\n",fMinDist2);
3868     printf("Min Distance to Bad Channel 3 = %2.1f\n",fMinDist3);
3869   }
3870   printf("    \n") ;
3871
3872 }
3873
3874