PWGGA/CaloTrackCorrelations/AliAnaPi0EbE.cxx

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