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