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