PWGGA/CaloTrackCorrelations/AliAnaPi0EbE.cxx

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