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