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