PWGGA/CaloTrackCorrelations/AliAnaParticleHadronCorrelation.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 particle - hadron correlations
  18 // Particle (for example direct gamma) must be found
  19 // in a previous analysis
  20 //
  21 //-- Author: Gustavo Conesa (LNF-INFN) (LPSC-IN2P3-CNRS)
  22 //           Yaxian Mao (LPSC-IN2P3-CNRS) and (CNWU) first usable implementation.
  23 //           Xiangrong Zhu (CNWU), implementtion of own mixing.
  24 //////////////////////////////////////////////////////////////////////////////
  25
  26
  27 // --- ROOT system ---
  28 //#include "TClonesArray.h"
  29 #include <TClass.h>
  30 #include <TMath.h>
  31 #include <TH2F.h>
  32 #include <TDatabasePDG.h>
  33
  34 //---- ANALYSIS system ----
  35 #include "AliNeutralMesonSelection.h"
  36 #include "AliAnaParticleHadronCorrelation.h"
  37 #include "AliCaloTrackReader.h"
  38 #include "AliAODPWG4ParticleCorrelation.h"
  39 #include "AliFiducialCut.h"
  40 #include "AliVTrack.h"
  41 #include "AliVCluster.h"
  42 #include "AliMCAnalysisUtils.h"
  43 #include "TParticle.h"
  44 #include "AliStack.h"
  45 #include "AliAODMCParticle.h"
  46 #include "AliMixedEvent.h"
  47 #include "AliAnalysisManager.h"
  48 #include "AliInputEventHandler.h"
  49 #include "AliEventplane.h"
  50
  51 ClassImp(AliAnaParticleHadronCorrelation)
  52
  53
  54 //___________________________________________________________________
  55   AliAnaParticleHadronCorrelation::AliAnaParticleHadronCorrelation():
  56     AliAnaCaloTrackCorrBaseClass(),
  57     fFillAODWithReferences(0),      fCheckLeadingWithNeutralClusters(0),
  58     fMaxAssocPt(1000.),             fMinAssocPt(0.),
  59     fDeltaPhiMaxCut(0.),            fDeltaPhiMinCut(0.),
  60     fSelectIsolated(0),             fMakeSeveralUE(0),
  61     fUeDeltaPhiMaxCut(0.),          fUeDeltaPhiMinCut(0.),
  62     fPi0AODBranchName(""),          fNeutralCorr(0),
  63     fPi0Trigger(0),                 fDecayTrigger(0),
  64     fNDecayBits(0),                 fDecayBits(),
  65     fMakeAbsoluteLeading(0),        fMakeNearSideLeading(0),
  66     fLeadingTriggerIndex(-1),       fHMPIDCorrelation(0),  fFillBradHisto(0),
  67     fNAssocPtBins(0),               fAssocPtBinLimit(),
  68     fCorrelVzBin(0),
  69     fListMixTrackEvents(),          fListMixCaloEvents(),
  70     fUseMixStoredInReader(0),       fFillNeutralEventMixPool(0),
  71     fM02MaxCut(0),                  fM02MinCut(0),
  72     fSelectLeadingHadronAngle(0),   fFillLeadHadOppositeHisto(0),
  73     fMinLeadHadPhi(0),              fMaxLeadHadPhi(0),
  74     fMinLeadHadPt(0),               fMaxLeadHadPt(0),
  75     fFillEtaGapsHisto(1),           fFillMomImbalancePtAssocBinsHisto(0),
  76     fMCGenTypeMin(0),               fMCGenTypeMax(0),
  77     fTrackVector(),                 fMomentum(),
  78     fDecayMom1(),                   fDecayMom2(),
  79     //Histograms
  80     fhPtTriggerInput(0),            fhPtTriggerSSCut(0),
  81     fhPtTriggerIsoCut(0),           fhPtTriggerFidCut(0),
  82     fhPtTrigger(0),                 fhPtTriggerVtxBC0(0),
  83     fhPtTriggerVzBin(0),            fhPtTriggerBin(0),
  84     fhPhiTrigger(0),                fhEtaTrigger(0),
  85     fhPtTriggerMC(),
  86     fhPtDecayTrigger(),             fhPtDecayTriggerMC(),
  87     fhPtTriggerCentrality(0),       fhPtTriggerEventPlane(0),
  88     fhTriggerEventPlaneCentrality(0),
  89     fhPtTriggerMixed(0),            fhPtTriggerMixedVzBin(0), fhPtTriggerMixedBin(0),
  90     fhPhiTriggerMixed(0),           fhEtaTriggerMixed(0),
  91     fhPtLeadingOppositeHadron(0),   fhPtDiffPhiLeadingOppositeHadron(0), fhPtDiffEtaLeadingOppositeHadron(0),
  92     fhPtNoLeadingOppositeHadron(0), fhEtaPhiNoLeadingOppositeHadron(0),
  93     fhDeltaPhiDeltaEtaCharged(0),
  94     fhPhiCharged(0),                fhEtaCharged(0),
  95     fhDeltaPhiCharged(0),           fhDeltaEtaCharged(0),
  96     fhDeltaPhiChargedPt(0),         fhDeltaPhiUeChargedPt(0),
  97     fhUePart(0),
  98     fhXECharged(0),                 fhXECharged_Cone2(0),      fhXEUeCharged(0),
  99     fhXEPosCharged(0),              fhXENegCharged(0),
 100     fhPtHbpXECharged(0),            fhPtHbpXECharged_Cone2(0), fhPtHbpXEUeCharged(0),
 101     fhZTCharged(0),                 fhZTUeCharged(0),
 102     fhZTPosCharged(0),              fhZTNegCharged(0),
 103     fhPtHbpZTCharged(0),            fhPtHbpZTUeCharged(0),
 104     fhXEChargedMC(),                fhDeltaPhiChargedMC(),
 105     fhDeltaPhiDeltaEtaChargedPtA3GeV(0),
 106     fhDeltaPhiChargedPtA3GeV(0),    fhDeltaEtaChargedPtA3GeV(0),
 107     //Pile-Up
 108     fhDeltaPhiChargedPileUp(),      fhDeltaEtaChargedPileUp(),
 109     fhDeltaPhiChargedPtA3GeVPileUp(), fhDeltaEtaChargedPtA3GeVPileUp(),
 110     fhXEChargedPileUp(),            fhXEUeChargedPileUp(),
 111     fhZTChargedPileUp(),            fhZTUeChargedPileUp(),
 112     fhPtTrigChargedPileUp(),
 113     fhDeltaPhiChargedOtherBC(),     fhDeltaPhiChargedPtA3GeVOtherBC(),
 114     fhXEChargedOtherBC(),           fhXEUeChargedOtherBC(),
 115     fhZTChargedOtherBC(),           fhZTUeChargedOtherBC(),
 116     fhPtTrigChargedOtherBC(),
 117     fhDeltaPhiChargedBC0(),         fhDeltaPhiChargedPtA3GeVBC0(),
 118     fhXEChargedBC0(),               fhXEUeChargedBC0(),
 119     fhZTChargedBC0(),               fhZTUeChargedBC0(),
 120     fhPtTrigChargedBC0(),
 121     fhDeltaPhiChargedVtxBC0(),      fhDeltaPhiChargedPtA3GeVVtxBC0(),
 122     fhXEChargedVtxBC0(),            fhXEUeChargedVtxBC0(),
 123     fhZTChargedVtxBC0(),            fhZTUeChargedVtxBC0(),
 124     fhPtTrigChargedVtxBC0(),
 125     fhDeltaPhiUeLeftCharged(0),
 126     fhDeltaPhiUeLeftUpCharged(0),   fhDeltaPhiUeRightUpCharged(0),
 127     fhDeltaPhiUeLeftDownCharged(0), fhDeltaPhiUeRightDownCharged(0),
 128     fhXEUeLeftCharged(0),
 129     fhXEUeLeftUpCharged(0),         fhXEUeRightUpCharged(0),
 130     fhXEUeLeftDownCharged(0),       fhXEUeRightDownCharged(0),
 131     fhPtHbpXEUeLeftCharged(0),      fhZTUeLeftCharged(0),
 132     fhPtHbpZTUeLeftCharged(0),
 133     fhPtTrigPout(0),                fhPtTrigCharged(0),
 134     fhDeltaPhiChargedMult(0x0),     fhDeltaEtaChargedMult(0x0),
 135     fhXEMult(0x0),                  fhXEUeMult(0x0),
 136     fhZTMult(0x0),                  fhZTUeMult(0x0),
 137     fhAssocPtBkg(0),                fhDeltaPhiDeltaEtaAssocPtBin(0),
 138     fhDeltaPhiAssocPtBin(0),
 139     fhDeltaPhiAssocPtBinDEta08(0),  fhDeltaPhiAssocPtBinDEta0(0),
 140     fhDeltaPhiAssocPtBinHMPID(0),   fhDeltaPhiAssocPtBinHMPIDAcc(0),
 141     fhDeltaPhiBradAssocPtBin(0),    fhDeltaPhiBrad(0),
 142     fhXEAssocPtBin(0),              fhZTAssocPtBin(0),
 143     fhXEVZ(0),                      fhZTVZ(0),
 144     fhDeltaPhiDeltaEtaNeutral(0),
 145     fhPhiNeutral(0),                fhEtaNeutral(0),
 146     fhDeltaPhiNeutral(0),           fhDeltaEtaNeutral(0),
 147     fhDeltaPhiNeutralPt(0),         fhDeltaPhiUeNeutralPt(0),
 148     fhXENeutral(0),                 fhXEUeNeutral(0),
 149     fhPtHbpXENeutral(0),            fhPtHbpXEUeNeutral(0),
 150     fhZTNeutral(0),                 fhZTUeNeutral(0),
 151     fhPtHbpZTNeutral(0),            fhPtHbpZTUeNeutral(0),
 152     fhDeltaPhiUeLeftNeutral(0),     fhXEUeLeftNeutral(0),
 153     fhPtHbpXEUeLeftNeutral(0),      fhZTUeLeftNeutral(0),
 154     fhPtHbpZTUeLeftNeutral(0),      fhPtPi0DecayRatio(0),
 155     fhDeltaPhiPi0DecayCharged(0),   fhXEPi0DecayCharged(0),        fhZTPi0DecayCharged(0),
 156     fhDeltaPhiPi0DecayNeutral(0),   fhXEPi0DecayNeutral(0),        fhZTPi0DecayNeutral(0),
 157     fhDeltaPhiDecayCharged(),       fhXEDecayCharged(),            fhZTDecayCharged(),
 158     fhDeltaPhiDecayChargedAssocPtBin(),
 159     fhMCPtTrigger(),                fhMCPhiTrigger(),              fhMCEtaTrigger(),
 160     fhMCPtTriggerNotLeading(),      fhMCPhiTriggerNotLeading(),    fhMCEtaTriggerNotLeading(),
 161     fhMCEtaCharged(),               fhMCPhiCharged(),
 162     fhMCDeltaEtaCharged(),          fhMCDeltaPhiCharged(),
 163     fhMCDeltaPhiDeltaEtaCharged(),  fhMCDeltaPhiChargedPt(),
 164     fhMCPtXECharged(),              fhMCPtXEUeCharged(),
 165     fhMCPtXEUeLeftCharged(),
 166     fhMCPtHbpXECharged(),           fhMCPtHbpXEUeCharged(),
 167     fhMCPtHbpXEUeLeftCharged(),
 168     fhMCUePart(),
 169     fhMCPtZTCharged(),              fhMCPtZTUeCharged(),
 170     fhMCPtZTUeLeftCharged(),
 171     fhMCPtHbpZTCharged(),           fhMCPtHbpZTUeCharged(),
 172     fhMCPtHbpZTUeLeftCharged(),
 173     fhMCPtTrigPout(),               fhMCPtAssocDeltaPhi(),
 174     //Mixing
 175     fhNEventsTrigger(0),            fhNtracksMB(0),                 fhNclustersMB(0),
 176     fhMixDeltaPhiCharged(0),        fhMixDeltaPhiDeltaEtaCharged(0),
 177     fhMixXECharged(0),              fhMixXEUeCharged(0),            fhMixHbpXECharged(0),
 178     fhMixDeltaPhiChargedAssocPtBin(),
 179     fhMixDeltaPhiChargedAssocPtBinDEta08(),
 180     fhMixDeltaPhiChargedAssocPtBinDEta0(),
 181     fhMixDeltaPhiDeltaEtaChargedAssocPtBin(),
 182     fhEventBin(0),                  fhEventMixBin(0),               fhEventMBBin(0)
 183 {
 184   //Default Ctor
 185
 186   //Initialize parameters
 187   InitParameters();
 188
 189   for(Int_t i = 0; i < fgkNmcTypes; i++)
 190   {
 191     fhPtTriggerMC[i] = 0;
 192     fhXEChargedMC[i] = 0;
 193     fhDeltaPhiChargedMC[i] = 0;
 194     for(Int_t ib = 0; ib < 4; ib++)  fhPtDecayTriggerMC[ib][i] = 0;
 195   }
 196
 197   for(Int_t ib = 0; ib < 4; ib++)  fhPtDecayTrigger[ib] = 0;
 198
 199   for(Int_t i = 0; i < 7; i++)
 200   {
 201     fhPtTriggerPileUp             [i] = 0 ;
 202     fhDeltaPhiChargedPileUp       [i] = 0 ; fhDeltaEtaChargedPileUp       [i] = 0 ;
 203     fhXEChargedPileUp             [i] = 0 ; fhXEUeChargedPileUp           [i] = 0 ;
 204     fhZTChargedPileUp             [i] = 0 ; fhZTUeChargedPileUp           [i] = 0 ;
 205     fhPtTrigChargedPileUp         [i] = 0 ;
 206     fhDeltaPhiChargedPtA3GeVPileUp[i] = 0 ; fhDeltaEtaChargedPtA3GeVPileUp[i] = 0 ;
 207   }
 208
 209 }
 210
 211 //_________________________________________________________________
 212 AliAnaParticleHadronCorrelation::~AliAnaParticleHadronCorrelation()
 213 {
 214   // Remove event containers
 215
 216   if(DoOwnMix())
 217   {
 218     if(fListMixTrackEvents)
 219     {
 220       for(Int_t iz=0; iz < GetNZvertBin(); iz++)
 221       {
 222         for(Int_t ic=0; ic < GetNCentrBin(); ic++)
 223         {
 224           for(Int_t irp=0; irp<GetNRPBin(); irp++)
 225           {
 226             Int_t bin = GetEventMixBin(ic, iz, irp);
 227             fListMixTrackEvents[bin]->Delete() ;
 228             delete fListMixTrackEvents[bin] ;
 229           }
 230         }
 231       }
 232     }
 233
 234     delete[] fListMixTrackEvents;
 235
 236     if(fListMixCaloEvents)
 237     {
 238       for(Int_t iz=0; iz < GetNZvertBin(); iz++)
 239       {
 240         for(Int_t ic=0; ic < GetNCentrBin(); ic++)
 241         {
 242           for(Int_t irp=0; irp<GetNRPBin(); irp++)
 243           {
 244             Int_t bin = GetEventMixBin(ic, iz, irp);
 245             fListMixCaloEvents[bin]->Delete() ;
 246             delete fListMixCaloEvents[bin] ;
 247           }
 248         }
 249       }
 250     }
 251
 252     delete[] fListMixCaloEvents;
 253
 254   }
 255 }
 256
 257 //__________________________________________________________________________________________________________________________________________
 258 void AliAnaParticleHadronCorrelation::FillChargedAngularCorrelationHistograms(Float_t ptAssoc,  Float_t ptTrig,      Int_t   bin,
 259                                                                               Float_t phiAssoc, Float_t phiTrig,     Float_t deltaPhi,
 260                                                                               Float_t etaAssoc, Float_t etaTrig,
 261                                                                               Int_t   decayTag, Float_t hmpidSignal, Int_t  outTOF,
 262                                                                               Int_t   cen,      Int_t   mcTag)
 263 {
 264   // Fill angular correlation related histograms
 265
 266   Float_t deltaEta    = etaTrig-etaAssoc;
 267   Float_t deltaPhiOrg = phiTrig-phiAssoc;
 268
 269   fhEtaCharged       ->Fill(ptAssoc,etaAssoc);
 270   fhPhiCharged       ->Fill(ptAssoc,phiAssoc);
 271   fhDeltaEtaCharged  ->Fill(ptTrig ,deltaEta);
 272   fhDeltaPhiCharged  ->Fill(ptTrig ,deltaPhi);
 273   fhDeltaPhiChargedPt->Fill(ptAssoc, deltaPhi);
 274   fhDeltaPhiDeltaEtaCharged->Fill(deltaPhi, deltaEta);
 275
 276   if(ptAssoc > 3 )
 277   {
 278     fhDeltaEtaChargedPtA3GeV        ->Fill(ptTrig  ,deltaEta);
 279     fhDeltaPhiChargedPtA3GeV        ->Fill(ptTrig  ,deltaPhi);
 280     fhDeltaPhiDeltaEtaChargedPtA3GeV->Fill(deltaPhi, deltaEta);
 281   }
 282
 283   // Pile up studies
 284
 285   if(IsPileUpAnalysisOn())
 286   {
 287     if     (outTOF==1)
 288     {
 289       fhDeltaPhiChargedOtherBC->Fill(ptTrig ,deltaPhi) ;
 290       if(ptAssoc > 3 ) fhDeltaPhiChargedPtA3GeVOtherBC->Fill(ptTrig ,deltaPhi) ;
 291     }
 292     else if(outTOF==0)
 293     {
 294         fhDeltaPhiChargedBC0->Fill(ptTrig ,deltaPhi) ;
 295         if(ptAssoc > 3 ) fhDeltaPhiChargedPtA3GeVBC0->Fill(ptTrig ,deltaPhi) ;
 296     }
 297
 298     Int_t vtxBC = GetReader()->GetVertexBC();
 299     if(vtxBC == 0 || vtxBC==AliVTrack::kTOFBCNA)
 300     {
 301       fhDeltaPhiChargedVtxBC0->Fill(ptTrig ,deltaPhi) ;
 302       if(ptAssoc > 3 ) fhDeltaPhiChargedPtA3GeVVtxBC0->Fill(ptTrig ,deltaPhi) ;
 303     }
 304
 305     if(GetReader()->IsPileUpFromSPD())               { fhDeltaEtaChargedPileUp[0]->Fill(ptTrig ,deltaEta) ; fhDeltaPhiChargedPileUp[0]->Fill(ptTrig ,deltaPhi) ; }
 306     if(GetReader()->IsPileUpFromEMCal())             { fhDeltaEtaChargedPileUp[1]->Fill(ptTrig ,deltaEta) ; fhDeltaPhiChargedPileUp[1]->Fill(ptTrig ,deltaPhi) ; }
 307     if(GetReader()->IsPileUpFromSPDOrEMCal())        { fhDeltaEtaChargedPileUp[2]->Fill(ptTrig ,deltaEta) ; fhDeltaPhiChargedPileUp[2]->Fill(ptTrig ,deltaPhi) ; }
 308     if(GetReader()->IsPileUpFromSPDAndEMCal())       { fhDeltaEtaChargedPileUp[3]->Fill(ptTrig ,deltaEta) ; fhDeltaPhiChargedPileUp[3]->Fill(ptTrig ,deltaPhi) ; }
 309     if(GetReader()->IsPileUpFromSPDAndNotEMCal())    { fhDeltaEtaChargedPileUp[4]->Fill(ptTrig ,deltaEta) ; fhDeltaPhiChargedPileUp[4]->Fill(ptTrig ,deltaPhi) ; }
 310     if(GetReader()->IsPileUpFromEMCalAndNotSPD())    { fhDeltaEtaChargedPileUp[5]->Fill(ptTrig ,deltaEta) ; fhDeltaPhiChargedPileUp[5]->Fill(ptTrig ,deltaPhi) ; }
 311     if(GetReader()->IsPileUpFromNotSPDAndNotEMCal()) { fhDeltaEtaChargedPileUp[6]->Fill(ptTrig ,deltaEta) ; fhDeltaPhiChargedPileUp[6]->Fill(ptTrig ,deltaPhi) ; }
 312
 313     if(ptAssoc > 3 )
 314     {
 315       if(GetReader()->IsPileUpFromSPD())               { fhDeltaEtaChargedPtA3GeVPileUp[0]->Fill(ptTrig ,deltaEta) ; fhDeltaPhiChargedPtA3GeVPileUp[0]->Fill(ptTrig ,deltaPhi) ; }
 316       if(GetReader()->IsPileUpFromEMCal())             { fhDeltaEtaChargedPtA3GeVPileUp[1]->Fill(ptTrig ,deltaEta) ; fhDeltaPhiChargedPtA3GeVPileUp[1]->Fill(ptTrig ,deltaPhi) ; }
 317       if(GetReader()->IsPileUpFromSPDOrEMCal())        { fhDeltaEtaChargedPtA3GeVPileUp[2]->Fill(ptTrig ,deltaEta) ; fhDeltaPhiChargedPtA3GeVPileUp[2]->Fill(ptTrig ,deltaPhi) ; }
 318       if(GetReader()->IsPileUpFromSPDAndEMCal())       { fhDeltaEtaChargedPtA3GeVPileUp[3]->Fill(ptTrig ,deltaEta) ; fhDeltaPhiChargedPtA3GeVPileUp[3]->Fill(ptTrig ,deltaPhi) ; }
 319       if(GetReader()->IsPileUpFromSPDAndNotEMCal())    { fhDeltaEtaChargedPtA3GeVPileUp[4]->Fill(ptTrig ,deltaEta) ; fhDeltaPhiChargedPtA3GeVPileUp[4]->Fill(ptTrig ,deltaPhi) ; }
 320       if(GetReader()->IsPileUpFromEMCalAndNotSPD())    { fhDeltaEtaChargedPtA3GeVPileUp[5]->Fill(ptTrig ,deltaEta) ; fhDeltaPhiChargedPtA3GeVPileUp[5]->Fill(ptTrig ,deltaPhi) ; }
 321       if(GetReader()->IsPileUpFromNotSPDAndNotEMCal()) { fhDeltaEtaChargedPtA3GeVPileUp[6]->Fill(ptTrig ,deltaEta) ; fhDeltaPhiChargedPtA3GeVPileUp[6]->Fill(ptTrig ,deltaPhi) ; }
 322     }
 323   }
 324
 325   if(IsDataMC())
 326   {
 327     Int_t mcIndex = GetMCTagHistogramIndex(mcTag);
 328     fhDeltaPhiChargedMC[mcIndex]->Fill(ptTrig , deltaPhi);
 329     if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCDecayPairLost) )
 330     {
 331       // check index in GetMCTagIndexHistogram
 332       if     ( mcIndex == 2 ) fhDeltaPhiChargedMC[8]->Fill(ptTrig , deltaPhi); // pi0 decay
 333       else if( mcIndex == 4 ) fhDeltaPhiChargedMC[9]->Fill(ptTrig , deltaPhi); // eta decay
 334     }
 335   }
 336
 337   if(fDecayTrigger && decayTag > 0)
 338   {
 339     for(Int_t ibit = 0; ibit<fNDecayBits; ibit++)
 340     {
 341       if(GetNeutralMesonSelection()->CheckDecayBit(decayTag,fDecayBits[ibit])) fhDeltaPhiDecayCharged[ibit]->Fill(ptTrig,deltaPhi);
 342     }
 343   }
 344
 345   Double_t  dphiBrad = -100;
 346   if(fFillBradHisto)
 347   {
 348     dphiBrad = atan2(sin(deltaPhiOrg), cos(deltaPhiOrg))/TMath::Pi();//-1 to 1
 349     if( TMath::Abs(dphiBrad) > 0.325 && TMath::Abs(dphiBrad) < 0.475 )  //Hardcoded values, BAD, FIXME
 350     {
 351       fhAssocPtBkg->Fill(ptTrig, ptAssoc);
 352     }
 353
 354     if( dphiBrad < -1./3 ) dphiBrad += 2;
 355     fhDeltaPhiBrad->Fill(ptTrig, dphiBrad);
 356   }
 357
 358   // Fill histograms in bins of associated particle pT
 359   if(bin>=0)
 360   {
 361       fhDeltaPhiDeltaEtaAssocPtBin    [bin]->Fill(deltaPhi,deltaEta);
 362
 363       fhDeltaPhiAssocPtBin            [bin]->Fill(ptTrig, deltaPhi);
 364
 365     if(fFillEtaGapsHisto)
 366     {
 367       if(TMath::Abs(deltaEta)> 0.8)
 368         fhDeltaPhiAssocPtBinDEta08      [bin]->Fill(ptTrig, deltaPhi);
 369
 370       if(TMath::Abs(deltaEta)< 0.01)
 371         fhDeltaPhiAssocPtBinDEta0       [bin]->Fill(ptTrig, deltaPhi);
 372     }
 373
 374     if (fFillBradHisto)
 375       fhDeltaPhiBradAssocPtBin        [bin]->Fill(ptTrig, dphiBrad);
 376
 377     if(fDecayTrigger && decayTag > 0 && fNDecayBits > 0 &&
 378        GetNeutralMesonSelection()->CheckDecayBit(decayTag,fDecayBits[0]))
 379       fhDeltaPhiDecayChargedAssocPtBin[bin]->Fill(ptTrig, deltaPhi);
 380
 381     if(fHMPIDCorrelation)
 382     {
 383       if( hmpidSignal > 0 )
 384       {
 385         //printf("Track pt %f with HMPID signal %f \n",pt,hmpidSignal);
 386         fhDeltaPhiAssocPtBinHMPID[bin]->Fill(ptTrig, deltaPhi);
 387       }
 388
 389       if(phiAssoc > 5*TMath::DegToRad() && phiAssoc < 20*TMath::DegToRad())
 390       {
 391         //printf("Track pt %f in HMPID acceptance phi %f \n ",pt,phi*TMath::RadToDeg() );
 392         fhDeltaPhiAssocPtBinHMPIDAcc[bin]->Fill(ptTrig, deltaPhi);
 393       }
 394     }
 395   }
 396
 397   //fill different multiplicity/centrality histogram
 398   if(IsHighMultiplicityAnalysisOn() && cen >= 0 && cen < GetNCentrBin())
 399   {
 400     fhDeltaPhiChargedMult[cen]->Fill(ptTrig,deltaPhi);
 401     fhDeltaEtaChargedMult[cen]->Fill(ptTrig,deltaEta);
 402   }
 403 }
 404
 405 //___________________________________________________________________________________________________________________________________
 406 Bool_t AliAnaParticleHadronCorrelation::FillChargedMCCorrelationHistograms(Float_t mcAssocPt, Float_t mcAssocPhi, Float_t mcAssocEta,
 407                                                                            Float_t mcTrigPt,  Float_t mcTrigPhi,  Float_t mcTrigEta,
 408                                                                            Int_t histoIndex, Bool_t lostDecayPair)
 409 {
 410   // Fill MC histograms independently of AOD or ESD
 411
 412   Bool_t lead = kTRUE;
 413
 414   // In case we requested the trigger to be a leading particle,
 415   // check if this is true at the MC level.
 416   // Not sure if it is correct to skip or not skip this.
 417   // Absolute leading?
 418   if( fMakeAbsoluteLeading && mcAssocPt > mcTrigPt )     lead = kFALSE; // skip event
 419
 420   // Skip this event if near side associated particle pt larger than trigger
 421   if( mcAssocPhi < 0 ) mcAssocPhi+=TMath::TwoPi();
 422
 423   Float_t mcdeltaPhi= mcTrigPhi-mcAssocPhi;
 424   if(mcdeltaPhi <= -TMath::PiOver2()) mcdeltaPhi+=TMath::TwoPi();
 425   if(mcdeltaPhi > 3*TMath::PiOver2()) mcdeltaPhi-=TMath::TwoPi();
 426
 427   if( fMakeNearSideLeading)
 428   {
 429     if( mcAssocPt > mcTrigPt && mcdeltaPhi < TMath::PiOver2() ) lead = kFALSE; // skip event
 430   }
 431
 432   //
 433   // Select only hadrons in pt range
 434   if ( mcAssocPt < fMinAssocPt || mcAssocPt > fMaxAssocPt ) return lead ; // exclude but continue
 435   if ( mcAssocPt < GetReader()->GetCTSPtMin())              return lead ;
 436
 437
 438   //
 439   // Remove trigger itself for correlation when use charged triggers
 440   if(TMath::Abs(mcAssocPt -mcTrigPt ) < 1e-6 &&
 441      mcdeltaPhi                       < 1e-6 &&
 442      TMath::Abs(mcAssocEta-mcTrigEta) < 1e-6)            return lead ; // exclude but continue
 443
 444   Float_t mcxE    =-mcAssocPt/mcTrigPt*TMath::Cos(mcdeltaPhi);// -(mcAssocPx*pxprim+mcAssocPy*pyprim)/(mcTrigPt*mcTrigPt);
 445   Float_t mchbpXE =-100 ;
 446   if(mcxE > 0 ) mchbpXE = TMath::Log(1./mcxE);
 447
 448   Float_t mczT    = mcAssocPt/mcTrigPt ;
 449   Float_t mchbpZT =-100 ;
 450   if(mczT > 0 ) mchbpZT = TMath::Log(1./mczT);
 451
 452   Double_t mcpout = mcAssocPt*TMath::Sin(mcdeltaPhi) ;
 453
 454   AliDebug(1,Form("Charged hadron: track Pt %f, track Phi %f, phi trigger %f. Cuts:  delta phi  %2.2f < %2.2f < %2.2f",
 455                   mcAssocPt,mcAssocPhi, mcTrigPhi,fDeltaPhiMinCut, mcdeltaPhi, fDeltaPhiMaxCut));
 456
 457
 458   // Fill Histograms
 459   fhMCEtaCharged     [histoIndex]->Fill(mcAssocPt, mcAssocEta);
 460   fhMCPhiCharged     [histoIndex]->Fill(mcAssocPt, mcAssocPhi);
 461   fhMCDeltaEtaCharged[histoIndex]->Fill(mcTrigPt , mcTrigEta-mcAssocEta);
 462   fhMCDeltaPhiCharged[histoIndex]->Fill(mcTrigPt , mcdeltaPhi);
 463   fhMCPtAssocDeltaPhi[histoIndex]->Fill(mcAssocPt, mcdeltaPhi);
 464
 465   fhMCDeltaPhiDeltaEtaCharged[histoIndex]->Fill(mcdeltaPhi,mcTrigEta-mcAssocEta);
 466
 467   //delta phi cut for correlation
 468   if( (mcdeltaPhi > fDeltaPhiMinCut) && (mcdeltaPhi < fDeltaPhiMaxCut) )
 469   {
 470     fhMCDeltaPhiChargedPt[histoIndex]->Fill(mcAssocPt,mcdeltaPhi);
 471     fhMCPtXECharged      [histoIndex]->Fill(mcTrigPt, mcxE);
 472     fhMCPtHbpXECharged   [histoIndex]->Fill(mcTrigPt, mchbpXE);
 473     fhMCPtZTCharged      [histoIndex]->Fill(mcTrigPt, mczT);
 474     fhMCPtHbpZTCharged   [histoIndex]->Fill(mcTrigPt, mchbpZT);
 475     fhMCPtTrigPout       [histoIndex]->Fill(mcTrigPt, mcpout) ;
 476   }
 477
 478   if(lostDecayPair)
 479   {
 480     // check index in GetMCTagIndexHistogram
 481     if(histoIndex == 2 && 8 >= fMCGenTypeMin && 8 <= fMCGenTypeMax )
 482     {
 483       // pi0 decay
 484       fhMCEtaCharged     [8]->Fill(mcAssocPt, mcAssocEta);
 485       fhMCPhiCharged     [8]->Fill(mcAssocPt, mcAssocPhi);
 486       fhMCDeltaEtaCharged[8]->Fill(mcTrigPt , mcTrigEta-mcAssocEta);
 487       fhMCDeltaPhiCharged[8]->Fill(mcTrigPt , mcdeltaPhi);
 488       fhMCPtAssocDeltaPhi[8]->Fill(mcAssocPt, mcdeltaPhi);
 489
 490       fhMCDeltaPhiDeltaEtaCharged[8]->Fill(mcdeltaPhi,mcTrigEta-mcAssocEta);
 491
 492       //delta phi cut for correlation
 493       if( (mcdeltaPhi > fDeltaPhiMinCut) && (mcdeltaPhi < fDeltaPhiMaxCut) )
 494       {
 495         fhMCDeltaPhiChargedPt[8]->Fill(mcAssocPt,mcdeltaPhi);
 496         fhMCPtXECharged      [8]->Fill(mcTrigPt, mcxE);
 497         fhMCPtHbpXECharged   [8]->Fill(mcTrigPt, mchbpXE);
 498         fhMCPtZTCharged      [8]->Fill(mcTrigPt, mczT);
 499         fhMCPtHbpZTCharged   [8]->Fill(mcTrigPt, mchbpZT);
 500         fhMCPtTrigPout       [8]->Fill(mcTrigPt, mcpout) ;
 501       }
 502     } // pi0 decay lost pair
 503     if(histoIndex == 4 && 9 >= fMCGenTypeMin && 9 <= fMCGenTypeMax )
 504     {
 505       // eta decay
 506       fhMCEtaCharged     [9]->Fill(mcAssocPt, mcAssocEta);
 507       fhMCPhiCharged     [9]->Fill(mcAssocPt, mcAssocPhi);
 508       fhMCDeltaEtaCharged[9]->Fill(mcTrigPt , mcTrigEta-mcAssocEta);
 509       fhMCDeltaPhiCharged[9]->Fill(mcTrigPt , mcdeltaPhi);
 510       fhMCPtAssocDeltaPhi[9]->Fill(mcAssocPt, mcdeltaPhi);
 511
 512       fhMCDeltaPhiDeltaEtaCharged[9]->Fill(mcdeltaPhi,mcTrigEta-mcAssocEta);
 513
 514       //delta phi cut for correlation
 515       if( (mcdeltaPhi > fDeltaPhiMinCut) && (mcdeltaPhi < fDeltaPhiMaxCut) )
 516       {
 517         fhMCDeltaPhiChargedPt[9]->Fill(mcAssocPt,mcdeltaPhi);
 518         fhMCPtXECharged      [9]->Fill(mcTrigPt, mcxE);
 519         fhMCPtHbpXECharged   [9]->Fill(mcTrigPt, mchbpXE);
 520         fhMCPtZTCharged      [9]->Fill(mcTrigPt, mczT);
 521         fhMCPtHbpZTCharged   [9]->Fill(mcTrigPt, mchbpZT);
 522         fhMCPtTrigPout       [9]->Fill(mcTrigPt, mcpout) ;
 523       }
 524     } // eta decay lost pair
 525   }
 526   // Underlying event
 527
 528   // Right
 529   if ( (mcdeltaPhi > fUeDeltaPhiMinCut) && (mcdeltaPhi < fUeDeltaPhiMaxCut) )
 530   {
 531     //Double_t randomphi = gRandom->Uniform(TMath::Pi()/2,3*TMath::Pi()/2);
 532     Double_t randomphi = gRandom->Uniform(fDeltaPhiMinCut,fDeltaPhiMaxCut);
 533     Double_t mcUexE = -(mcAssocPt/mcTrigPt)*TMath::Cos(randomphi);
 534     Double_t mcUezT =   mcAssocPt/mcTrigPt;
 535
 536     if(mcUexE < 0.)
 537       AliWarning(Form("Careful!!, negative xE %2.2f for right UE cos(dPhi %2.2f) = %2.2f, check correlation dPhi limits %f to %f",
 538                       mcUexE,randomphi*TMath::RadToDeg(),TMath::Cos(randomphi),fDeltaPhiMinCut*TMath::RadToDeg(),fDeltaPhiMaxCut*TMath::RadToDeg()));
 539
 540     fhMCPtXEUeCharged[histoIndex]->Fill(mcTrigPt,mcUexE);
 541     if(mcUexE > 0) fhMCPtHbpXEUeCharged[histoIndex]->Fill(mcTrigPt,TMath::Log(1/mcUexE));
 542
 543     fhMCPtZTUeCharged[histoIndex]->Fill(mcTrigPt,mcUezT);
 544     if(mcUezT > 0) fhMCPtHbpZTUeCharged[histoIndex]->Fill(mcTrigPt,TMath::Log(1/mcUezT));
 545
 546     fhMCUePart[histoIndex]->Fill(mcTrigPt);
 547
 548     if(lostDecayPair)
 549     {
 550       // check index in GetMCTagIndexHistogram
 551       if(histoIndex == 2 && 8 >= fMCGenTypeMin && 8 <= fMCGenTypeMax )
 552       {
 553         // pi0 decay
 554         fhMCPtXEUeCharged[8]->Fill(mcTrigPt,mcUexE);
 555         if(mcUexE > 0) fhMCPtHbpXEUeCharged[8]->Fill(mcTrigPt,TMath::Log(1/mcUexE));
 556
 557         fhMCPtZTUeCharged[8]->Fill(mcTrigPt,mcUezT);
 558         if(mcUezT > 0) fhMCPtHbpZTUeCharged[8]->Fill(mcTrigPt,TMath::Log(1/mcUezT));
 559
 560         fhMCUePart[8]->Fill(mcTrigPt);
 561       }
 562       else if(histoIndex == 4 && 9 >= fMCGenTypeMin && 9 <= fMCGenTypeMax )
 563       {
 564         // eta decay
 565         fhMCPtXEUeCharged[9]->Fill(mcTrigPt,mcUexE);
 566         if(mcUexE > 0) fhMCPtHbpXEUeCharged[9]->Fill(mcTrigPt,TMath::Log(1/mcUexE));
 567
 568         fhMCPtZTUeCharged[9]->Fill(mcTrigPt,mcUezT);
 569         if(mcUezT > 0) fhMCPtHbpZTUeCharged[9]->Fill(mcTrigPt,TMath::Log(1/mcUezT));
 570
 571         fhMCUePart[9]->Fill(mcTrigPt);
 572       }
 573     }
 574   }
 575
 576   if(fMakeSeveralUE)
 577   {
 578     // Left
 579     if((mcdeltaPhi<-fUeDeltaPhiMinCut) || (mcdeltaPhi >2*fUeDeltaPhiMaxCut))
 580     {
 581       Double_t randomphi = gRandom->Uniform(fDeltaPhiMinCut,fDeltaPhiMaxCut);
 582       Double_t mcUexE = -(mcAssocPt/mcTrigPt)*TMath::Cos(randomphi);
 583       Double_t mcUezT =   mcAssocPt/mcTrigPt;
 584
 585       if(mcUexE < 0.)
 586         AliWarning(Form("Careful!!, negative xE %2.2f for left UE cos(dPhi %2.2f) = %2.2f, check correlation dPhi limits %f to %f",
 587                         mcUexE,randomphi*TMath::RadToDeg(),TMath::Cos(randomphi),fDeltaPhiMinCut*TMath::RadToDeg(),fDeltaPhiMaxCut*TMath::RadToDeg()));
 588
 589       fhMCPtXEUeLeftCharged[histoIndex]->Fill(mcTrigPt,mcUexE);
 590       if(mcUexE > 0) fhMCPtHbpXEUeLeftCharged[histoIndex]->Fill(mcTrigPt,TMath::Log(1/mcUexE));
 591
 592       fhMCPtZTUeLeftCharged[histoIndex]->Fill(mcTrigPt,mcUezT);
 593       if(mcUezT > 0) fhMCPtHbpZTUeLeftCharged[histoIndex]->Fill(mcTrigPt,TMath::Log(1/mcUezT));
 594
 595       if(lostDecayPair)
 596       {
 597         // check index in GetMCTagIndexHistogram
 598         if(histoIndex == 2  && 8 >= fMCGenTypeMin && 8 <= fMCGenTypeMax )
 599         {
 600           // pi0 decay
 601           fhMCPtXEUeLeftCharged[8]->Fill(mcTrigPt,mcUexE);
 602           if(mcUexE > 0) fhMCPtHbpXEUeLeftCharged[8]->Fill(mcTrigPt,TMath::Log(1/mcUexE));
 603
 604           fhMCPtZTUeLeftCharged[8]->Fill(mcTrigPt,mcUezT);
 605           if(mcUezT > 0) fhMCPtHbpZTUeLeftCharged[8]->Fill(mcTrigPt,TMath::Log(1/mcUezT));
 606         }
 607         else if(histoIndex == 4  && 9 >= fMCGenTypeMin && 9 <= fMCGenTypeMax )
 608         {
 609           // eta decay
 610           fhMCPtXEUeLeftCharged[9]->Fill(mcTrigPt,mcUexE);
 611           if(mcUexE > 0) fhMCPtHbpXEUeLeftCharged[9]->Fill(mcTrigPt,TMath::Log(1/mcUexE));
 612
 613           fhMCPtZTUeLeftCharged[9]->Fill(mcTrigPt,mcUezT);
 614           if(mcUezT > 0) fhMCPtHbpZTUeLeftCharged[9]->Fill(mcTrigPt,TMath::Log(1/mcUezT));
 615         }
 616       }
 617     }
 618   }
 619
 620   return lead;
 621 }
 622
 623 //___________________________________________________________________________________________________________________
 624 void AliAnaParticleHadronCorrelation::FillChargedMomentumImbalanceHistograms(Float_t ptTrig,   Float_t ptAssoc,
 625                                                                              Float_t deltaPhi,
 626                                                                              Int_t   cen,      Int_t   charge,
 627                                                                              Int_t   bin,      Int_t   decayTag,
 628                                                                              Int_t   outTOF,   Int_t   mcTag)
 629
 630 {
 631   // Fill mostly momentum imbalance related histograms
 632
 633   Float_t zT =   ptAssoc/ptTrig ;
 634   Float_t xE   =-ptAssoc/ptTrig*TMath::Cos(deltaPhi); // -(px*pxTrig+py*pyTrig)/(ptTrig*ptTrig);
 635   Float_t pout = ptAssoc*TMath::Sin(deltaPhi) ;
 636
 637   if(xE < 0.)
 638     AliWarning(Form("Careful!!, negative xE %2.2f for right UE cos(dPhi %2.2f) = %2.2f, check correlation dPhi limits %f to %f",
 639                xE,deltaPhi*TMath::RadToDeg(),TMath::Cos(deltaPhi),fDeltaPhiMinCut*TMath::RadToDeg(),fDeltaPhiMaxCut*TMath::RadToDeg()));
 640
 641   Float_t hbpXE = -100;
 642   Float_t hbpZT = -100;
 643
 644   if(xE > 0 ) hbpXE = TMath::Log(1./xE);
 645   if(zT > 0 ) hbpZT = TMath::Log(1./zT);
 646
 647   fhXECharged         ->Fill(ptTrig , xE);
 648   fhPtHbpXECharged    ->Fill(ptTrig , hbpXE);
 649   fhZTCharged         ->Fill(ptTrig , zT);
 650   fhPtHbpZTCharged    ->Fill(ptTrig , hbpZT);
 651   fhPtTrigPout        ->Fill(ptTrig , pout) ;
 652   fhPtTrigCharged     ->Fill(ptTrig , ptAssoc) ;
 653   if((deltaPhi > 5*TMath::Pi()/6.)   && (deltaPhi < 7*TMath::Pi()/6.))
 654   {
 655     fhXECharged_Cone2         ->Fill(ptTrig , xE);
 656     fhPtHbpXECharged_Cone2    ->Fill(ptTrig , hbpXE);
 657   }
 658
 659   // MC
 660   if(IsDataMC())
 661   {
 662     Int_t mcIndex = GetMCTagHistogramIndex(mcTag);
 663     fhXEChargedMC[mcIndex]->Fill(ptTrig , xE);
 664     if( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCDecayPairLost) )
 665     {
 666       // check index in GetMCTagIndexHistogram
 667       if      ( mcIndex == 2 ) fhXEChargedMC[8]->Fill(ptTrig , xE); // pi0 decay
 668       else if ( mcIndex == 4 ) fhXEChargedMC[9]->Fill(ptTrig , xE); // eta decay
 669     }
 670   }
 671
 672   // Pile up studies
 673   if(IsPileUpAnalysisOn())
 674   {
 675     if     (outTOF==1)
 676     {
 677       fhXEChargedOtherBC    ->Fill(ptTrig,xE);
 678       fhZTChargedOtherBC    ->Fill(ptTrig,zT);
 679       fhPtTrigChargedOtherBC->Fill(ptTrig,ptAssoc);
 680     }
 681     else if(outTOF==0)
 682     {
 683       fhXEChargedBC0    ->Fill(ptTrig,xE);
 684       fhZTChargedBC0    ->Fill(ptTrig,zT);
 685       fhPtTrigChargedBC0->Fill(ptTrig,ptAssoc);
 686     }
 687
 688     Int_t vtxBC = GetReader()->GetVertexBC();
 689     if(vtxBC == 0 || vtxBC==AliVTrack::kTOFBCNA)
 690     {
 691       fhXEChargedVtxBC0    ->Fill(ptTrig,xE);
 692       fhZTChargedVtxBC0    ->Fill(ptTrig,zT);
 693       fhPtTrigChargedVtxBC0->Fill(ptTrig,ptAssoc);
 694     }
 695
 696     if(GetReader()->IsPileUpFromSPD())                { fhXEChargedPileUp[0]->Fill(ptTrig,xE); fhZTChargedPileUp[0]->Fill(ptTrig,zT); fhPtTrigChargedPileUp[0]->Fill(ptTrig,ptAssoc); }
 697     if(GetReader()->IsPileUpFromEMCal())              { fhXEChargedPileUp[1]->Fill(ptTrig,xE); fhZTChargedPileUp[1]->Fill(ptTrig,zT); fhPtTrigChargedPileUp[1]->Fill(ptTrig,ptAssoc); }
 698     if(GetReader()->IsPileUpFromSPDOrEMCal())         { fhXEChargedPileUp[2]->Fill(ptTrig,xE); fhZTChargedPileUp[2]->Fill(ptTrig,zT); fhPtTrigChargedPileUp[2]->Fill(ptTrig,ptAssoc); }
 699     if(GetReader()->IsPileUpFromSPDAndEMCal())        { fhXEChargedPileUp[3]->Fill(ptTrig,xE); fhZTChargedPileUp[3]->Fill(ptTrig,zT); fhPtTrigChargedPileUp[3]->Fill(ptTrig,ptAssoc); }
 700     if(GetReader()->IsPileUpFromSPDAndNotEMCal())     { fhXEChargedPileUp[4]->Fill(ptTrig,xE); fhZTChargedPileUp[4]->Fill(ptTrig,zT); fhPtTrigChargedPileUp[4]->Fill(ptTrig,ptAssoc); }
 701     if(GetReader()->IsPileUpFromEMCalAndNotSPD())     { fhXEChargedPileUp[5]->Fill(ptTrig,xE); fhZTChargedPileUp[5]->Fill(ptTrig,zT); fhPtTrigChargedPileUp[5]->Fill(ptTrig,ptAssoc); }
 702     if(GetReader()->IsPileUpFromNotSPDAndNotEMCal())  { fhXEChargedPileUp[6]->Fill(ptTrig,xE); fhZTChargedPileUp[6]->Fill(ptTrig,zT); fhPtTrigChargedPileUp[6]->Fill(ptTrig,ptAssoc); }
 703   }
 704
 705   if(fDecayTrigger && decayTag > 0)
 706   {
 707     for(Int_t ibit = 0; ibit<fNDecayBits; ibit++)
 708     {
 709       if(GetNeutralMesonSelection()->CheckDecayBit(decayTag,fDecayBits[ibit]))
 710       {
 711         fhXEDecayCharged[ibit]->Fill(ptTrig,xE);
 712         fhZTDecayCharged[ibit]->Fill(ptTrig,zT);
 713       }
 714     }
 715   } // photon decay pi0/eta trigger
 716
 717   if(bin >= 0 && fFillMomImbalancePtAssocBinsHisto)//away side
 718   {
 719     fhXEAssocPtBin[bin]->Fill(ptTrig, xE) ;
 720     fhZTAssocPtBin[bin]->Fill(ptTrig, zT) ;
 721   }
 722
 723   if(fCorrelVzBin)
 724   {
 725     Int_t vz = GetEventVzBin();
 726     fhXEVZ[vz]->Fill(ptTrig, xE) ;
 727     fhZTVZ[vz]->Fill(ptTrig, zT) ;
 728   }
 729
 730   if(charge > 0)
 731   {
 732     fhXEPosCharged->Fill(ptTrig,xE) ;
 733     fhZTPosCharged->Fill(ptTrig,zT) ;
 734   }
 735   else
 736   {
 737     fhXENegCharged->Fill(ptTrig,xE) ;
 738     fhZTNegCharged->Fill(ptTrig,zT) ;
 739   }
 740
 741   //fill different multiplicity/centrality histogram
 742   if(IsHighMultiplicityAnalysisOn() && cen >= 0 && cen < GetNCentrBin())
 743   {
 744     fhXEMult[cen]->Fill(ptTrig,xE);
 745     fhZTMult[cen]->Fill(ptTrig,zT);
 746   } //multiplicity/centrality events selection
 747 }
 748
 749 //_______________________________________________________________________________________________________________________
 750 void AliAnaParticleHadronCorrelation::FillChargedUnderlyingEventHistograms(Float_t ptTrig,   Float_t ptAssoc,
 751                                                                            Float_t deltaPhi, Int_t cen, Int_t outTOF)
 752 {
 753   // Fill underlying event histograms
 754
 755   fhUePart->Fill(ptTrig);
 756
 757   fhDeltaPhiUeChargedPt->Fill(ptAssoc,deltaPhi);
 758
 759   Double_t randomphi = gRandom->Uniform(fDeltaPhiMinCut,fDeltaPhiMaxCut);
 760   Double_t uexE = -(ptAssoc/ptTrig)*TMath::Cos(randomphi);
 761   Double_t uezT =   ptAssoc/ptTrig;
 762
 763   if(uexE < 0.)
 764     AliWarning(Form("Careful!!, negative xE %2.2f for right UE cos(dPhi %2.2f) = %2.2f, check correlation dPhi limits %f to %f",
 765                uexE,randomphi*TMath::RadToDeg(),TMath::Cos(randomphi),fDeltaPhiMinCut*TMath::RadToDeg(),fDeltaPhiMaxCut*TMath::RadToDeg()));
 766
 767   fhXEUeCharged->Fill(ptTrig,uexE);
 768   if(uexE > 0) fhPtHbpXEUeCharged->Fill(ptTrig,TMath::Log(1/uexE));
 769
 770   fhZTUeCharged->Fill(ptTrig,uezT);
 771   if(uezT > 0) fhPtHbpZTUeCharged->Fill(ptTrig,TMath::Log(1/uezT));
 772
 773   // Pile up studies
 774
 775   if(IsPileUpAnalysisOn())
 776   {
 777     if     (outTOF==1)
 778     {
 779       fhXEUeChargedOtherBC->Fill(ptTrig,uexE);
 780       fhZTUeChargedOtherBC->Fill(ptTrig,uezT);
 781     }
 782     else if(outTOF==0)
 783     {
 784       fhXEUeChargedBC0->Fill(ptTrig,uexE);
 785       fhZTUeChargedBC0->Fill(ptTrig,uezT);
 786     }
 787
 788     Int_t vtxBC = GetReader()->GetVertexBC();
 789     if(vtxBC == 0 || vtxBC==AliVTrack::kTOFBCNA)
 790     {
 791       fhXEUeChargedVtxBC0->Fill(ptTrig,uexE);
 792       fhZTUeChargedVtxBC0->Fill(ptTrig,uezT);
 793     }
 794
 795     if(GetReader()->IsPileUpFromSPD())               { fhXEUeChargedPileUp[0]->Fill(ptTrig,uexE); fhZTUeChargedPileUp[0]->Fill(ptTrig,uezT);}
 796     if(GetReader()->IsPileUpFromEMCal())             { fhXEUeChargedPileUp[1]->Fill(ptTrig,uexE); fhZTUeChargedPileUp[1]->Fill(ptTrig,uezT);}
 797     if(GetReader()->IsPileUpFromSPDOrEMCal())        { fhXEUeChargedPileUp[2]->Fill(ptTrig,uexE); fhZTUeChargedPileUp[2]->Fill(ptTrig,uezT);}
 798     if(GetReader()->IsPileUpFromSPDAndEMCal())       { fhXEUeChargedPileUp[3]->Fill(ptTrig,uexE); fhZTUeChargedPileUp[3]->Fill(ptTrig,uezT);}
 799     if(GetReader()->IsPileUpFromSPDAndNotEMCal())    { fhXEUeChargedPileUp[4]->Fill(ptTrig,uexE); fhZTUeChargedPileUp[4]->Fill(ptTrig,uezT);}
 800     if(GetReader()->IsPileUpFromEMCalAndNotSPD())    { fhXEUeChargedPileUp[5]->Fill(ptTrig,uexE); fhZTUeChargedPileUp[5]->Fill(ptTrig,uezT);}
 801     if(GetReader()->IsPileUpFromNotSPDAndNotEMCal()) { fhXEUeChargedPileUp[6]->Fill(ptTrig,uexE); fhZTUeChargedPileUp[6]->Fill(ptTrig,uezT);}
 802   }
 803
 804   //fill different multiplicity/centrality histogram
 805   if(IsHighMultiplicityAnalysisOn() && cen >= 0 && cen < GetNCentrBin())
 806   {
 807     fhXEUeMult[cen]->Fill(ptTrig,uexE);
 808     fhZTUeMult[cen]->Fill(ptTrig,uezT);
 809   }
 810 }
 811
 812 //_____________________________________________________________________________________________________
 813 void AliAnaParticleHadronCorrelation::FillChargedUnderlyingEventSidesHistograms(Float_t ptTrig,
 814                                                                                 Float_t ptAssoc,
 815                                                                                 Float_t deltaPhi)
 816 {
 817   // Fill underlying event histograms to the left and right of trigger
 818   // Right cone is the default UE.
 819
 820   if((deltaPhi<-fUeDeltaPhiMinCut) || (deltaPhi >2*fUeDeltaPhiMaxCut))
 821   {
 822     fhDeltaPhiUeLeftCharged->Fill(ptAssoc,deltaPhi);
 823     Double_t randomphi = gRandom->Uniform(fDeltaPhiMinCut,fDeltaPhiMaxCut);
 824     Double_t uexE = -(ptAssoc/ptTrig)*TMath::Cos(randomphi);
 825     Double_t uezT =   ptAssoc/ptTrig;
 826
 827     if(uexE < 0.)
 828       AliWarning(Form("Careful!!, negative xE %2.2f for left UE cos(dPhi %2.2f) = %2.2f, check correlation dPhi limits %f to %f",
 829                       uexE,randomphi*TMath::RadToDeg(),TMath::Cos(randomphi),fDeltaPhiMinCut*TMath::RadToDeg(),fDeltaPhiMaxCut*TMath::RadToDeg()));
 830
 831     fhXEUeLeftCharged->Fill(ptTrig,uexE);
 832     if(uexE > 0) fhPtHbpXEUeLeftCharged->Fill(ptTrig,TMath::Log(1/uexE));
 833
 834     fhZTUeLeftCharged->Fill(ptTrig,uezT);
 835     if(uezT > 0) fhPtHbpZTUeLeftCharged->Fill(ptTrig,TMath::Log(1/uezT));
 836
 837     fhDeltaPhiUeLeftCharged->Fill(ptAssoc, deltaPhi);
 838   }
 839
 840   if((deltaPhi<-fUeDeltaPhiMinCut) && (deltaPhi >-TMath::Pi()/2))
 841   {
 842     fhDeltaPhiUeLeftDownCharged->Fill(ptAssoc,deltaPhi);
 843     Double_t randomphi = gRandom->Uniform(fDeltaPhiMinCut,fDeltaPhiMaxCut);
 844     Double_t uexE = -(ptAssoc/ptTrig)*TMath::Cos(randomphi);
 845
 846     if(uexE < 0.)
 847       AliWarning(Form("Careful!!, negative xE %2.2f for left-down UE cos(dPhi %2.2f) = %2.2f, check correlation dPhi limits %f to %f",
 848                       uexE,randomphi*TMath::RadToDeg(),TMath::Cos(randomphi),fDeltaPhiMinCut*TMath::RadToDeg(),fDeltaPhiMaxCut*TMath::RadToDeg()));
 849
 850     fhXEUeLeftDownCharged->Fill(ptTrig,uexE);
 851   }
 852
 853   if((deltaPhi>2*fUeDeltaPhiMaxCut) && (deltaPhi <3*TMath::Pi()/2))
 854   {
 855     fhDeltaPhiUeLeftUpCharged->Fill(ptAssoc,deltaPhi);
 856     Double_t randomphi = gRandom->Uniform(fDeltaPhiMinCut,fDeltaPhiMaxCut);
 857     Double_t uexE = -(ptAssoc/ptTrig)*TMath::Cos(randomphi);
 858
 859     if(uexE < 0.)
 860       AliWarning(Form("Careful!!, negative xE %2.2f for left-up UE cos(dPhi %2.2f) = %2.2f, check correlation dPhi limits %f to %f",
 861                       uexE,randomphi*TMath::RadToDeg(),TMath::Cos(randomphi),fDeltaPhiMinCut*TMath::RadToDeg(),fDeltaPhiMaxCut*TMath::RadToDeg()));
 862
 863     fhXEUeLeftUpCharged->Fill(ptTrig,uexE);
 864   }
 865
 866   if((deltaPhi > TMath::Pi()/2) && (deltaPhi < fUeDeltaPhiMaxCut))
 867   {
 868     fhDeltaPhiUeRightUpCharged->Fill(ptAssoc,deltaPhi);
 869     Double_t randomphi = gRandom->Uniform(fDeltaPhiMinCut,fDeltaPhiMaxCut);
 870     Double_t uexE = -(ptAssoc/ptTrig)*TMath::Cos(randomphi);
 871
 872     if(uexE < 0.)
 873       AliWarning(Form("Careful!!, negative xE %2.2f for right-up UE cos(dPhi %2.2f) = %2.2f, check correlation dPhi limits %f to %f",
 874                       uexE,randomphi*TMath::RadToDeg(),TMath::Cos(randomphi),fDeltaPhiMinCut*TMath::RadToDeg(),fDeltaPhiMaxCut*TMath::RadToDeg()));
 875
 876     fhXEUeRightUpCharged->Fill(ptTrig,uexE);
 877   }
 878
 879   if((deltaPhi > fUeDeltaPhiMinCut) && (deltaPhi < TMath::Pi()/2))
 880   {
 881     fhDeltaPhiUeRightDownCharged->Fill(ptAssoc,deltaPhi);
 882     Double_t randomphi = gRandom->Uniform(fDeltaPhiMinCut,fDeltaPhiMaxCut);
 883     Double_t uexE = -(ptAssoc/ptTrig)*TMath::Cos(randomphi);
 884
 885     if(uexE < 0.)
 886       AliWarning(Form("Careful!!, negative xE %2.2f for right-down UE cos(dPhi %2.2f) = %2.2f, check correlation dPhi limits %f to %f",
 887              uexE,randomphi*TMath::RadToDeg(),TMath::Cos(randomphi),fDeltaPhiMinCut*TMath::RadToDeg(),fDeltaPhiMaxCut*TMath::RadToDeg()));
 888
 889     fhXEUeRightDownCharged->Fill(ptTrig,uexE);
 890   }
 891 }
 892
 893 //_____________________________________________________________________________________________________________________________________
 894 void AliAnaParticleHadronCorrelation::FillDecayPhotonCorrelationHistograms(Float_t ptAssoc, Float_t phiAssoc, Bool_t bChargedOrNeutral)
 895 {
 896   // Do correlation with decay photons of triggered pi0 or eta
 897
 898   // Calculate the correlation parameters
 899   Float_t ptDecay1 = fDecayMom1.Pt();
 900   Float_t ptDecay2 = fDecayMom2.Pt();
 901
 902   Float_t zTDecay1 = -100, zTDecay2 = -100;
 903   if(ptDecay1 > 0) zTDecay1 = ptAssoc/ptDecay1 ;
 904   if(ptDecay2 > 0) zTDecay2 = ptAssoc/ptDecay2 ;
 905
 906   Float_t deltaPhiDecay1 = fDecayMom1.Phi()-phiAssoc;
 907   if(deltaPhiDecay1< -TMath::PiOver2()) deltaPhiDecay1+=TMath::TwoPi();
 908   if(deltaPhiDecay1>3*TMath::PiOver2()) deltaPhiDecay1-=TMath::TwoPi();
 909
 910   Float_t deltaPhiDecay2 = fDecayMom2.Phi()-phiAssoc;
 911   if(deltaPhiDecay2< -TMath::PiOver2()) deltaPhiDecay2+=TMath::TwoPi();
 912   if(deltaPhiDecay2>3*TMath::PiOver2()) deltaPhiDecay2-=TMath::TwoPi();
 913
 914   Float_t xEDecay1   =-zTDecay1*TMath::Cos(deltaPhiDecay1); // -(px*pxTrig+py*pyTrig)/(ptTrig*ptTrig);
 915   Float_t xEDecay2   =-zTDecay2*TMath::Cos(deltaPhiDecay2); // -(px*pxTrig+py*pyTrig)/(ptTrig*ptTrig);
 916
 917   if(bChargedOrNeutral) // correlate with charges
 918   {
 919     fhDeltaPhiPi0DecayCharged->Fill(ptDecay1, deltaPhiDecay1);
 920     fhDeltaPhiPi0DecayCharged->Fill(ptDecay2, deltaPhiDecay2);
 921
 922     AliDebug(2,Form("deltaPhoton1 = %f, deltaPhoton2 = %f", deltaPhiDecay1, deltaPhiDecay2));
 923
 924     if( (deltaPhiDecay1 > fDeltaPhiMinCut) && ( deltaPhiDecay1 < fDeltaPhiMaxCut) )
 925     {
 926       fhZTPi0DecayCharged->Fill(ptDecay1,zTDecay1);
 927       fhXEPi0DecayCharged->Fill(ptDecay1,xEDecay1);
 928     }
 929     if( (deltaPhiDecay2 > fDeltaPhiMinCut) && ( deltaPhiDecay2 < fDeltaPhiMaxCut) )
 930     {
 931       fhZTPi0DecayCharged->Fill(ptDecay2,zTDecay2);
 932       fhXEPi0DecayCharged->Fill(ptDecay2,xEDecay2);
 933     }
 934   }
 935   else // correlate with neutrals
 936   {
 937     fhDeltaPhiPi0DecayNeutral->Fill(ptDecay1, deltaPhiDecay1);
 938     fhDeltaPhiPi0DecayNeutral->Fill(ptDecay2, deltaPhiDecay2);
 939
 940     AliDebug(2,Form("deltaPhoton1 = %f, deltaPhoton2 = %f", deltaPhiDecay1, deltaPhiDecay2));
 941
 942     if( (deltaPhiDecay1 > fDeltaPhiMinCut) && ( deltaPhiDecay1 < fDeltaPhiMaxCut) )
 943     {
 944       fhZTPi0DecayNeutral->Fill(ptDecay1,zTDecay1);
 945       fhXEPi0DecayNeutral->Fill(ptDecay1,xEDecay1);
 946     }
 947     if( (deltaPhiDecay2 > fDeltaPhiMinCut) && ( deltaPhiDecay2 < fDeltaPhiMaxCut) )
 948     {
 949       fhZTPi0DecayNeutral->Fill(ptDecay2,zTDecay2);
 950       fhXEPi0DecayNeutral->Fill(ptDecay2,xEDecay2);
 951     }
 952   }
 953 }
 954
 955 //_____________________________________________________________________________________________________________________________
 956 void AliAnaParticleHadronCorrelation::FillNeutralUnderlyingEventSidesHistograms(Float_t ptTrig,   Float_t ptAssoc,
 957                                                                                 Float_t zT,       Float_t hbpZT,
 958                                                                                 Float_t deltaPhi)
 959 {
 960   // Fill underlying event histograms to the left of trigger
 961   // Right is the default case
 962
 963   Double_t randomphi = gRandom->Uniform(fDeltaPhiMinCut,fDeltaPhiMaxCut);
 964
 965   Float_t xE  =-ptAssoc/ptTrig*TMath::Cos(randomphi); // -(px*pxTrig+py*pyTrig)/(ptTrig*ptTrig);
 966   Float_t hbpXE = -100;
 967   if(xE > 0 ) hbpXE = TMath::Log(1./xE);
 968
 969   if((deltaPhi<-fUeDeltaPhiMinCut) && (deltaPhi >-fUeDeltaPhiMaxCut))
 970   {
 971     fhDeltaPhiUeLeftNeutral->Fill(ptAssoc, deltaPhi);
 972     fhXEUeLeftNeutral      ->Fill(ptTrig , xE);
 973     fhPtHbpXEUeLeftNeutral ->Fill(ptTrig , hbpXE);
 974     fhZTUeLeftNeutral      ->Fill(ptTrig , zT);
 975     fhPtHbpZTUeLeftNeutral ->Fill(ptTrig , hbpZT);
 976   }
 977 }
 978
 979 //______________________________________________________
 980 void AliAnaParticleHadronCorrelation::FillEventMixPool()
 981 {
 982   // Fill the pool with tracks or clusters if requested
 983
 984   if ( !DoOwnMix() ) return;
 985
 986   FillChargedEventMixPool();
 987
 988   // Do the cluster pool filling only if requested
 989   // or in case of isolation cut using clusters in the cone.
 990   Bool_t isoCase = OnlyIsolated() && (GetIsolationCut()->GetParticleTypeInCone() != AliIsolationCut::kOnlyCharged);
 991
 992   if( !fFillNeutralEventMixPool && !isoCase) return;
 993
 994   FillNeutralEventMixPool();
 995 }
 996
 997 //_____________________________________________________________
 998 void AliAnaParticleHadronCorrelation::FillChargedEventMixPool()
 999 {
1000   // Mixed event pool filling for tracks
1001
1002   if(fUseMixStoredInReader && GetReader()->GetLastTracksMixedEvent() == GetEventNumber())
1003   {
1004     //printf("%s : Pool already filled for this event !!!\n",GetInputAODName().Data());
1005     return ; // pool filled previously for another trigger
1006   }
1007
1008   AliAnalysisManager   * manager      = AliAnalysisManager::GetAnalysisManager();
1009   AliInputEventHandler * inputHandler = dynamic_cast<AliInputEventHandler*>(manager->GetInputEventHandler());
1010
1011   if(!inputHandler) return ;
1012
1013   // Do mixing only with MB event (or the chosen mask), if not skip
1014   if( !(inputHandler->IsEventSelected( ) & GetReader()->GetMixEventTriggerMask()) ) return ;
1015
1016   Int_t eventBin = GetEventMixBin();
1017
1018   //Check that the bin exists, if not (bad determination of RP, centrality or vz bin) do nothing
1019   if(eventBin < 0) return;
1020
1021   fhEventMBBin->Fill(eventBin);
1022
1023   TObjArray * mixEventTracks = new TObjArray;
1024
1025   if(fUseMixStoredInReader)
1026   {
1027     fListMixTrackEvents[eventBin] = GetReader()->GetListWithMixedEventsForTracks(eventBin);
1028   }
1029
1030   if(!fListMixTrackEvents[eventBin]) fListMixTrackEvents[eventBin] = new TList();
1031
1032   //printf("%s ***** Pool Event bin : %d - nTracks %d\n",GetInputAODName().Data(),eventBin, GetCTSTracks()->GetEntriesFast());
1033
1034   TList * pool = fListMixTrackEvents[eventBin];
1035
1036   for(Int_t ipr = 0;ipr < GetCTSTracks()->GetEntriesFast() ; ipr ++ )
1037   {
1038     AliVTrack * track = (AliVTrack *) (GetCTSTracks()->At(ipr)) ;
1039
1040     fTrackVector.SetXYZ(track->Px(),track->Py(),track->Pz());
1041     Float_t pt   = fTrackVector.Pt();
1042
1043     //Select only hadrons in pt range
1044     if(pt < fMinAssocPt || pt > fMaxAssocPt) continue ;
1045
1046     AliAODPWG4Particle * mixedTrack = new AliAODPWG4Particle(track->Px(),track->Py(),track->Pz(),0);
1047     mixedTrack->SetDetectorTag(kCTS);
1048     mixedTrack->SetChargedBit(track->Charge()>0);
1049     mixEventTracks->Add(mixedTrack);
1050   }
1051
1052   fhNtracksMB->Fill(mixEventTracks->GetEntriesFast(),eventBin);
1053
1054   //Set the event number where the last event was added, to avoid double pool filling
1055   GetReader()->SetLastTracksMixedEvent(GetEventNumber());
1056
1057   //printf("Add event to pool with %d tracks \n ",mixEventTracks->GetEntries());
1058   pool->AddFirst(mixEventTracks);
1059   mixEventTracks = 0;
1060
1061   //printf("Pool size %d, max %d\n",pool->GetSize(), GetNMaxEvMix());
1062
1063   if(pool->GetSize() > GetNMaxEvMix())
1064   {//Remove last event
1065     TClonesArray * tmp = static_cast<TClonesArray*>(pool->Last()) ;
1066     pool->RemoveLast() ;
1067     delete tmp ;
1068   }
1069 }
1070
1071 //_____________________________________________________________
1072 void AliAnaParticleHadronCorrelation::FillNeutralEventMixPool()
1073 {
1074   // Mixed event pool filling for neutral clusters
1075   // Right now only for EMCAL and in isolation case
1076
1077   //printf("FillNeutralEventMixPool for %s\n",GetInputAODName().Data());
1078
1079   if(fUseMixStoredInReader && GetReader()->GetLastCaloMixedEvent() == GetEventNumber())
1080   {
1081     //printf("%s : Pool already filled for this event !!!\n",GetInputAODName().Data());
1082     return ; // pool filled previously for another trigger
1083   }
1084
1085   TObjArray * pl = GetEMCALClusters();
1086   //if (GetAODObjArrayName.Contains("PHOS") )pl    = GetPHOSClusters();
1087   //else                                     pl    = GetEMCALClusters();
1088
1089   AliAnalysisManager   * manager      = AliAnalysisManager::GetAnalysisManager();
1090   AliInputEventHandler * inputHandler = dynamic_cast<AliInputEventHandler*>(manager->GetInputEventHandler());
1091
1092   if(!inputHandler) return ;
1093
1094   // Do mixing only with MB event (or the chosen mask), if not skip
1095   if( !(inputHandler->IsEventSelected( ) & GetReader()->GetMixEventTriggerMask()) ) return ;
1096
1097   Int_t eventBin = GetEventMixBin();
1098
1099   //Check that the bin exists, if not (bad determination of RP, centrality or vz bin) do nothing
1100   if(eventBin < 0) return;
1101
1102   TObjArray * mixEventCalo = new TObjArray;
1103
1104   if(fUseMixStoredInReader)
1105   {
1106     fListMixCaloEvents[eventBin] = GetReader()->GetListWithMixedEventsForCalo(eventBin);
1107   }
1108
1109   if(!fListMixCaloEvents[eventBin]) fListMixCaloEvents[eventBin] = new TList();
1110
1111   TList * poolCalo = fListMixCaloEvents[eventBin];
1112
1113   for(Int_t ipr = 0;ipr <  pl->GetEntriesFast() ; ipr ++ )
1114   {
1115     AliVCluster * calo = (AliVCluster *) (pl->At(ipr)) ;
1116
1117     // remove matched clusters
1118     if( IsTrackMatched( calo, GetReader()->GetInputEvent() ) ) continue ;
1119
1120     //Cluster momentum calculation
1121     if(GetReader()->GetDataType() != AliCaloTrackReader::kMC)
1122     {
1123       calo->GetMomentum(fMomentum,GetVertex(0)) ;
1124     }//Assume that come from vertex in straight line
1125     else
1126     {
1127       Double_t vertex[]={0,0,0};
1128       calo->GetMomentum(fMomentum,vertex) ;
1129     }
1130
1131     Float_t pt = fMomentum.Pt();
1132     //Select only clusters in pt range
1133     if(pt < fMinAssocPt || pt > fMaxAssocPt) continue ;
1134
1135     AliAODPWG4Particle * mixedCalo = new AliAODPWG4Particle(fMomentum);
1136     mixedCalo->SetDetectorTag(kEMCAL);
1137     mixEventCalo->Add(mixedCalo);
1138   }
1139
1140   fhNclustersMB->Fill(mixEventCalo->GetEntriesFast(),eventBin);
1141
1142   //Set the event number where the last event was added, to avoid double pool filling
1143   GetReader()->SetLastCaloMixedEvent(GetEventNumber());
1144
1145   //printf("Add event to pool with %d clusters \n ",mixEventCalo->GetEntries());
1146   poolCalo->AddFirst(mixEventCalo);
1147   mixEventCalo = 0;
1148
1149   //printf("Pool size %d, max %d\n",poolCalo->GetSize(), GetNMaxEvMix());
1150
1151   if(poolCalo->GetSize() > GetNMaxEvMix())
1152   {//Remove last event
1153     TClonesArray * tmp = static_cast<TClonesArray*>(poolCalo->Last()) ;
1154     poolCalo->RemoveLast() ;
1155     delete tmp ;
1156   }
1157 }
1158
1159 //_________________________________________________________________________________________________________________
1160 Bool_t AliAnaParticleHadronCorrelation::FindLeadingOppositeHadronInWindow(AliAODPWG4ParticleCorrelation * particle)
1161 {
1162   // Select events where the leading charged particle in the opposite hemisphere
1163   // to the trigger particle is in a window centered at 180 from the trigger
1164
1165   Float_t etaTrig    = particle->Eta();
1166   Float_t ptTrig     = particle->Pt();
1167   Float_t phiTrig    = particle->Phi();
1168   if(phiTrig < 0 ) phiTrig+= TMath::TwoPi();
1169
1170   Float_t ptLeadHad  = 0 ;
1171   Float_t dphiLeadHad= -100 ;
1172   Float_t phiLeadHad = -100 ;
1173   Float_t etaLeadHad = -100 ;
1174   Int_t   nTrack     = 0;
1175
1176   for(Int_t ipr = 0;ipr < GetCTSTracks()->GetEntriesFast() ; ipr ++ )
1177   {
1178     AliVTrack * track = (AliVTrack *) (GetCTSTracks()->At(ipr)) ;
1179
1180     fTrackVector.SetXYZ(track->Px(),track->Py(),track->Pz());
1181
1182     Float_t pt   = fTrackVector.Pt();
1183     Float_t phi  = fTrackVector.Phi() ;
1184     if(phi < 0 ) phi+= TMath::TwoPi();
1185
1186     Float_t deltaPhi = phiTrig-phi;
1187     //
1188     // Calculate deltaPhi shift so that for the particles on the opposite side
1189     // it is defined between 90 and 270 degrees
1190     // Shift [-360,-90]  to [0, 270]
1191     // and [270,360] to [-90,0]
1192     if(deltaPhi <= -TMath::PiOver2()) deltaPhi+=TMath::TwoPi();
1193     if(deltaPhi > 3*TMath::PiOver2()) deltaPhi-=TMath::TwoPi();
1194
1195     if(pt > ptLeadHad && deltaPhi > TMath::PiOver2()) // in opposite hemisphere
1196     {
1197       ptLeadHad  = pt ;
1198       phiLeadHad = phi;
1199       dphiLeadHad= deltaPhi;
1200       etaLeadHad = fTrackVector.Eta();
1201       nTrack++;
1202     }
1203   }// track loop
1204
1205   if(fFillLeadHadOppositeHisto)
1206   {
1207     if(nTrack == 0)
1208     {
1209       fhPtNoLeadingOppositeHadron    ->Fill(ptTrig);
1210       fhEtaPhiNoLeadingOppositeHadron->Fill(etaTrig,phiTrig);
1211     }
1212     else
1213     {
1214       fhPtLeadingOppositeHadron       ->Fill(ptTrig,  ptLeadHad);
1215       fhPtDiffPhiLeadingOppositeHadron->Fill(ptTrig,dphiLeadHad);
1216       fhPtDiffEtaLeadingOppositeHadron->Fill(ptTrig, etaLeadHad-etaTrig);
1217     }
1218   }
1219
1220
1221     AliDebug(1,Form("pT %2.2f, phi %2.2f, eta %2.2f, nTracks away %d, total tracks %d",
1222                     ptLeadHad,phiLeadHad*TMath::RadToDeg(),etaLeadHad,nTrack, GetTrackMultiplicity()));
1223     AliDebug(1,Form("\t pT trig %2.2f, Dphi (trigger-hadron) %2.2f, Deta (trigger-hadron) %2.2f",
1224            ptTrig, dphiLeadHad*TMath::RadToDeg(), etaLeadHad-etaTrig));
1225     AliDebug(1,Form("\t cuts pT: min %2.2f, max %2.2f; DPhi: min %2.2f, max %2.2f",
1226                     fMinLeadHadPt,fMaxLeadHadPt,fMinLeadHadPhi*TMath::RadToDeg(),fMaxLeadHadPhi*TMath::RadToDeg()));
1227
1228
1229   // reject the trigger if the leading hadron is not in the requested pt or phi window and
1230
1231   if( nTrack == 0 ) return kFALSE; // No track found in opposite hemisphere
1232
1233   if( ptLeadHad < fMinLeadHadPt || ptLeadHad > fMaxLeadHadPt ) return kFALSE;
1234
1235   //printf("Accept leading hadron pT \n");
1236
1237   if( dphiLeadHad < fMinLeadHadPhi || dphiLeadHad > fMaxLeadHadPhi ) return kFALSE;
1238
1239   //printf("Accept leading hadron phi \n");
1240
1241
1242   return kTRUE ;
1243 }
1244
1245 //____________________________________________________________
1246 TObjString* AliAnaParticleHadronCorrelation::GetAnalysisCuts()
1247 {
1248   //Save parameters used for analysis
1249   TString parList ; //this will be list of parameters used for this analysis.
1250   const Int_t buffersize = 560;
1251   char onePar[buffersize] ;
1252
1253   snprintf(onePar,buffersize,"--- AliAnaPaticleHadronCorrelation ---:") ;
1254   parList+=onePar ;
1255   snprintf(onePar,buffersize," %3.2f < Pt associated < %3.2f; ", fMinAssocPt,   fMaxAssocPt) ;
1256   parList+=onePar ;
1257   snprintf(onePar,buffersize," %3.2f < Phi trigger particle-Hadron < %3.2f; ",    fDeltaPhiMinCut,   fDeltaPhiMaxCut) ;
1258   parList+=onePar ;
1259   snprintf(onePar,buffersize," %3.2f < Phi trigger particle-UeHadron <  %3.2f; ", fUeDeltaPhiMinCut, fUeDeltaPhiMaxCut) ;
1260   parList+=onePar ;
1261   snprintf(onePar,buffersize,"Isolated Trigger?  %d;",    fSelectIsolated) ;
1262   parList+=onePar ;
1263   snprintf(onePar,buffersize,"Several UE?  %d;",          fMakeSeveralUE) ;
1264   parList+=onePar ;
1265   snprintf(onePar,buffersize,"Name of AOD Pi0 Branch %s;", fPi0AODBranchName.Data());
1266   parList+=onePar ;
1267   snprintf(onePar,buffersize,"Do Decay-hadron correlation ?  pi0 %d, decay %d;", fPi0Trigger, fDecayTrigger) ;
1268   parList+=onePar ;
1269   snprintf(onePar,buffersize,"Select absolute leading for cluster triggers ? %d or Near Side Leading %d;",
1270            fMakeAbsoluteLeading, fMakeNearSideLeading) ;
1271   parList+=onePar ;
1272   snprintf(onePar,buffersize,"Associated particle pt bins  %d: ", fNAssocPtBins) ;
1273   parList+=onePar ;
1274   for (Int_t ibin = 0; ibin<fNAssocPtBins; ibin++) {
1275     snprintf(onePar,buffersize,"bin %d = [%2.1f,%2.1f];", ibin, fAssocPtBinLimit[ibin], fAssocPtBinLimit[ibin+1]) ;
1276   }
1277   parList+=onePar ;
1278
1279   //Get parameters set in base class.
1280   parList += GetBaseParametersList() ;
1281
1282   //Get parameters set in FiducialCut class (not available yet)
1283   //parlist += GetFidCut()->GetFidCutParametersList()
1284
1285   return new TObjString(parList) ;
1286
1287 }
1288
1289 //________________________________________________________________
1290 TList *  AliAnaParticleHadronCorrelation::GetCreateOutputObjects()
1291 {
1292   // Create histograms to be saved in output file and
1293   // store them in fOutputContainer
1294
1295   TList * outputContainer = new TList() ;
1296   outputContainer->SetName("CorrelationHistos") ;
1297
1298   Int_t   nptbins = GetHistogramRanges()->GetHistoPtBins(); Int_t  nphibins = GetHistogramRanges()->GetHistoPhiBins(); Int_t   netabins = GetHistogramRanges()->GetHistoEtaBins();
1299   Float_t ptmax   = GetHistogramRanges()->GetHistoPtMax();  Float_t phimax  = GetHistogramRanges()->GetHistoPhiMax();  Float_t etamax   = GetHistogramRanges()->GetHistoEtaMax();
1300   Float_t ptmin   = GetHistogramRanges()->GetHistoPtMin();  Float_t phimin  = GetHistogramRanges()->GetHistoPhiMin();  Float_t etamin   = GetHistogramRanges()->GetHistoEtaMin();
1301
1302   Int_t  ndeltaphibins = GetHistogramRanges()->GetHistoDeltaPhiBins(); Int_t   ndeltaetabins = GetHistogramRanges()->GetHistoDeltaEtaBins();
1303   Float_t deltaphimax  = GetHistogramRanges()->GetHistoDeltaPhiMax();  Float_t deltaetamax   = GetHistogramRanges()->GetHistoDeltaEtaMax();
1304   Float_t deltaphimin  = GetHistogramRanges()->GetHistoDeltaPhiMin();  Float_t deltaetamin   = GetHistogramRanges()->GetHistoDeltaEtaMin();
1305
1306   Int_t ntrbins = GetHistogramRanges()->GetHistoTrackMultiplicityBins(); Int_t nclbins = GetHistogramRanges()->GetHistoNClustersBins();
1307   Int_t trmax   = GetHistogramRanges()->GetHistoTrackMultiplicityMax();  Int_t clmax   = GetHistogramRanges()->GetHistoNClustersMax();
1308   Int_t trmin   = GetHistogramRanges()->GetHistoTrackMultiplicityMin();  Int_t clmin   = GetHistogramRanges()->GetHistoNClustersMin();
1309
1310   Int_t  nxeztbins = GetHistogramRanges()->GetHistoRatioBins();  Int_t  nhbpbins = GetHistogramRanges()->GetHistoHBPBins();
1311   Float_t xeztmax  = GetHistogramRanges()->GetHistoRatioMax();   Float_t hbpmax  = GetHistogramRanges()->GetHistoHBPMax();
1312   Float_t xeztmin  = GetHistogramRanges()->GetHistoRatioMin();   Float_t hbpmin  = GetHistogramRanges()->GetHistoHBPMin();
1313
1314   Int_t nMixBins = GetNCentrBin()*GetNZvertBin()*GetNRPBin();
1315
1316   TString nameMC[]     = {"Photon","Pi0","Pi0Decay","Eta","EtaDecay","OtherDecay","Electron","Hadron","Pi0DecayLostPair","EtaDecayLostPair"};
1317   TString pileUpName[] = {"SPD","EMCAL","SPDOrEMCAL","SPDAndEMCAL","SPDAndNotEMCAL","EMCALAndNotSPD","NotSPDAndNotEMCAL"} ;
1318
1319   // For vz dependent histograms, if option ON
1320   Int_t   nz  = 1  ;
1321   if(fCorrelVzBin) nz = GetNZvertBin();
1322   TString sz  = "" ;
1323   TString tz  = "" ;
1324
1325   // Fill histograms for neutral clusters in mixing?
1326   Bool_t isoCase = OnlyIsolated() && (GetIsolationCut()->GetParticleTypeInCone() != AliIsolationCut::kOnlyCharged);
1327   Bool_t neutralMix = fFillNeutralEventMixPool || isoCase ;
1328
1329   fhPtTriggerInput  = new TH1F("hPtTriggerInput","Input trigger #it{p}_{T}", nptbins,ptmin,ptmax);
1330   fhPtTriggerInput->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1331   outputContainer->Add(fhPtTriggerInput);
1332
1333   if( fM02MaxCut > 0 && fM02MinCut > 0 )
1334   {
1335     fhPtTriggerSSCut  = new TH1F("hPtTriggerSSCut","Trigger #it{p}_{T} after #lambda^{2}_{0} cut", nptbins,ptmin,ptmax);
1336     fhPtTriggerSSCut->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1337     outputContainer->Add(fhPtTriggerSSCut);
1338   }
1339
1340   if( OnlyIsolated() )
1341   {
1342     fhPtTriggerIsoCut  = new TH1F("hPtTriggerIsoCut","Trigger #it{p}_{T} after isolation (and #lambda^{2}_{0} cut)", nptbins,ptmin,ptmax);
1343     fhPtTriggerIsoCut->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1344     outputContainer->Add(fhPtTriggerIsoCut);
1345   }
1346
1347   fhPtTriggerFidCut  = new TH1F("hPtTriggerFidCut","Trigger #it{p}_{T} after fiducial (isolation and #lambda^{2}_{0}) cut", nptbins,ptmin,ptmax);
1348   fhPtTriggerFidCut->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1349   outputContainer->Add(fhPtTriggerFidCut);
1350
1351   fhPtTrigger  = new TH1F("hPtTrigger","#it{p}_{T} distribution of trigger particles (after opposite hadron leading cut and rest)", nptbins,ptmin,ptmax);
1352   fhPtTrigger->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1353   outputContainer->Add(fhPtTrigger);
1354
1355   if(IsDataMC())
1356   {
1357     for(Int_t i=0; i < fgkNmcTypes; i++)
1358     {
1359       fhPtTriggerMC[i]  = new TH1F(Form("hPtTrigger_MC%s",nameMC[i].Data()),
1360                                    Form("#it{p}_{T} distribution of trigger particles, trigger origin is %s",nameMC[i].Data()),
1361                                    nptbins,ptmin,ptmax);
1362       fhPtTriggerMC[i]->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1363       outputContainer->Add(fhPtTriggerMC[i]);
1364     }
1365   }
1366
1367   if(fDecayTrigger)
1368   {
1369     for(Int_t ibit = 0; ibit < fNDecayBits; ibit++)
1370     {
1371       fhPtDecayTrigger[ibit]  = new TH1F(Form("hPtDecayTrigger_bit%d",fDecayBits[ibit]),
1372                                          Form("#it{p}_{T} distribution of trigger particles, decay Bit %d",fDecayBits[ibit]),
1373                                          nptbins,ptmin,ptmax);
1374       fhPtDecayTrigger[ibit]->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1375       outputContainer->Add(fhPtDecayTrigger[ibit]);
1376
1377       if(IsDataMC())
1378       {
1379         for(Int_t i=0; i < fgkNmcTypes; i++)
1380         {
1381           fhPtDecayTriggerMC[ibit][i]  = new TH1F(Form("hPtDecayTrigger_bit%d_MC%s",fDecayBits[ibit], nameMC[i].Data()),
1382                                             Form("#it{p}_{T} distribution of trigger particles, decay Bit %d, trigger origin is %s",fDecayBits[ibit], nameMC[i].Data()),
1383                                             nptbins,ptmin,ptmax);
1384           fhPtDecayTriggerMC[ibit][i]->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1385           outputContainer->Add(fhPtDecayTriggerMC[ibit][i]);
1386         }
1387       }
1388     }
1389   }
1390
1391   if(fCorrelVzBin)
1392   {
1393     fhPtTriggerVzBin  = new TH2F("hPtTriggerVzBin","#it{p}_{T} distribution of trigger particles vs vz bin", nptbins,ptmin,ptmax,GetNZvertBin(),0,GetNZvertBin());
1394     fhPtTriggerVzBin->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1395     fhPtTriggerVzBin->SetYTitle("#it{v}_{#it{z}} bin");
1396     outputContainer->Add(fhPtTriggerVzBin);
1397   }
1398
1399   fhPtTriggerBin  = new TH2F ("hPtTriggerBin","#it{p}_{T} distribution of trigger particles", nptbins,ptmin,ptmax,nMixBins,0,nMixBins);
1400   fhPtTriggerBin->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1401   fhPtTriggerBin->SetYTitle("Bin");
1402   outputContainer->Add(fhPtTriggerBin);
1403
1404   fhPhiTrigger  = new TH2F ("hPhiTrigger","#phi distribution of trigger Particles",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
1405   fhPhiTrigger->SetYTitle("#phi (rad)");
1406   outputContainer->Add(fhPhiTrigger);
1407
1408   fhEtaTrigger  = new TH2F ("hEtaTrigger","#eta distribution of trigger",nptbins,ptmin,ptmax, netabins,etamin,etamax);
1409   fhEtaTrigger->SetYTitle("#eta ");
1410   outputContainer->Add(fhEtaTrigger);
1411
1412   if(IsHighMultiplicityAnalysisOn())
1413   {
1414     fhPtTriggerCentrality   = new TH2F("hPtTriggerCentrality","Trigger particle #it{p}_{T} vs centrality",nptbins,ptmin,ptmax,100,0.,100) ;
1415     fhPtTriggerCentrality->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1416     fhPtTriggerCentrality->SetYTitle("Centrality (%)");
1417     outputContainer->Add(fhPtTriggerCentrality) ;
1418
1419     fhPtTriggerEventPlane  = new TH2F("hPtTriggerEventPlane","Trigger particle #it{p}_{T} vs event plane angle",nptbins,ptmin,ptmax, 100,0.,TMath::Pi()) ;
1420     fhPtTriggerEventPlane->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1421     fhPtTriggerEventPlane->SetXTitle("EP angle (rad)");
1422     outputContainer->Add(fhPtTriggerEventPlane) ;
1423
1424     fhTriggerEventPlaneCentrality  = new TH2F("hTriggerEventPlaneCentrality","Trigger particle centrality vs event plane angle",100,0.,100,100,0.,TMath::Pi()) ;
1425     fhTriggerEventPlaneCentrality->SetXTitle("Centrality (%)");
1426     fhTriggerEventPlaneCentrality->SetYTitle("EP angle (rad)");
1427     outputContainer->Add(fhTriggerEventPlaneCentrality) ;
1428   }
1429
1430   // Leading hadron in oposite side
1431   if(fFillLeadHadOppositeHisto)
1432   {
1433     fhPtLeadingOppositeHadron  = new TH2F("hPtTriggerPtLeadingOppositeHadron","Leading hadron opposite to trigger vs trigger #it{p}_{T}",
1434                                           nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
1435     fhPtLeadingOppositeHadron->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1436     fhPtLeadingOppositeHadron->SetYTitle("#it{p}_{T}^{lead hadron} (GeV/#it{c})");
1437     outputContainer->Add(fhPtLeadingOppositeHadron);
1438
1439     fhPtNoLeadingOppositeHadron  = new TH1F("hPtTriggerNoLeadingOppositeHadron","No Leading hadron opposite to trigger #it{p}_{T}",
1440                                           nptbins,ptmin,ptmax);
1441     fhPtNoLeadingOppositeHadron->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1442     outputContainer->Add(fhPtNoLeadingOppositeHadron);
1443
1444     fhEtaPhiNoLeadingOppositeHadron  = new TH2F("hEtaPhiTriggerNoLeadingOppositeHadron","No Leading hadron opposite to trigger #eta:#phi",
1445                                             netabins,etamin,etamax,nphibins,phimin,phimax);
1446     fhEtaPhiNoLeadingOppositeHadron->SetXTitle("#eta");
1447     fhEtaPhiNoLeadingOppositeHadron->SetYTitle("#phi");
1448     outputContainer->Add(fhEtaPhiNoLeadingOppositeHadron);
1449
1450
1451     fhPtDiffPhiLeadingOppositeHadron  = new TH2F("hPtTriggerDiffPhiTriggerLeadingOppositeHadron","#phi_{trigger}-#phi_{leading opposite hadron} vs #it{p}_{T}^{trig}",
1452                                                  nptbins,ptmin,ptmax,ndeltaphibins,deltaphimin,deltaphimax);
1453     fhPtDiffPhiLeadingOppositeHadron->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1454     fhPtDiffPhiLeadingOppositeHadron->SetYTitle("#phi_{trigger}-#phi_{leading opposite hadron} (rad)");
1455     outputContainer->Add(fhPtDiffPhiLeadingOppositeHadron);
1456
1457     fhPtDiffEtaLeadingOppositeHadron  = new TH2F("hPtTriggerDiffEtaTriggerPhiLeadingOppositeHadron","#eta_{trigger}-#eta_{leading opposite hadron} vs #it{p}_{T}^{trig}",
1458                                                  nptbins,ptmin,ptmax,ndeltaetabins,deltaetamin,deltaetamax);
1459     fhPtDiffEtaLeadingOppositeHadron->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1460     fhPtDiffEtaLeadingOppositeHadron->SetYTitle("#eta_{trigger}-#eta_{leading opposite hadron}");
1461     outputContainer->Add(fhPtDiffEtaLeadingOppositeHadron);
1462   }
1463
1464   //Correlation with charged hadrons
1465
1466   fhDeltaPhiDeltaEtaCharged  = new TH2F
1467   ("hDeltaPhiDeltaEtaCharged","#eta_{trigger} - #eta_{h^{#pm}} vs #phi_{trigger} - #phi_{h^{#pm}}",
1468    ndeltaphibins ,deltaphimin,deltaphimax,ndeltaetabins,deltaetamin,deltaetamax);
1469   fhDeltaPhiDeltaEtaCharged->SetXTitle("#Delta #phi (rad)");
1470   fhDeltaPhiDeltaEtaCharged->SetYTitle("#Delta #eta");
1471
1472   fhDeltaPhiDeltaEtaChargedPtA3GeV  = new TH2F
1473   ("hDeltaPhiDeltaEtaChargedPtA3GeV","#eta_{trigger} - #eta_{h^{#pm}} vs #phi_{trigger} - #phi_{h^{#pm}, #it{p}_{TA}>3 GeV/#it{c}}",
1474    ndeltaphibins ,deltaphimin,deltaphimax,ndeltaetabins,deltaetamin,deltaetamax);
1475   fhDeltaPhiDeltaEtaChargedPtA3GeV->SetXTitle("#Delta #phi (rad)");
1476   fhDeltaPhiDeltaEtaChargedPtA3GeV->SetYTitle("#Delta #eta");
1477
1478   fhPhiCharged  = new TH2F
1479   ("hPhiCharged","#phi_{h^{#pm}}  vs #it{p}_{T #pm}",
1480    nptbins,ptmin,ptmax,180,0,TMath::TwoPi());
1481   fhPhiCharged->SetYTitle("#phi_{h^{#pm}} (rad)");
1482   fhPhiCharged->SetXTitle("#it{p}_{T #pm} (GeV/#it{c})");
1483
1484   fhEtaCharged  = new TH2F
1485   ("hEtaCharged","#eta_{h^{#pm}}  vs #it{p}_{T #pm}",
1486    nptbins,ptmin,ptmax,100,-1.,1.);
1487   fhEtaCharged->SetYTitle("#eta_{h^{#pm}} (rad)");
1488   fhEtaCharged->SetXTitle("#it{p}_{T #pm} (GeV/#it{c})");
1489
1490   fhDeltaPhiCharged  = new TH2F
1491   ("hDeltaPhiCharged","#phi_{trigger} - #phi_{h^{#pm}} vs #it{p}_{T trigger}",
1492    nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1493   fhDeltaPhiCharged->SetYTitle("#Delta #phi (rad)");
1494   fhDeltaPhiCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1495
1496   fhDeltaPhiChargedPtA3GeV  = new TH2F
1497   ("hDeltaPhiChargedPtA3GeV","#phi_{trigger} - #phi_{h^{#pm}} vs #it{p}_{T trigger}, #it{p}_{TA}>3 GeV/#it{c}",
1498    nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1499   fhDeltaPhiChargedPtA3GeV->SetYTitle("#Delta #phi (rad)");
1500   fhDeltaPhiChargedPtA3GeV->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1501
1502
1503   fhDeltaPhiChargedPt  = new TH2F
1504   ("hDeltaPhiChargedPt","#phi_{trigger} - #phi_{#h^{#pm}} vs #it{p}_{T h^{#pm}}",
1505    nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1506   fhDeltaPhiChargedPt->SetYTitle("#Delta #phi (rad)");
1507   fhDeltaPhiChargedPt->SetXTitle("#it{p}_{T h^{#pm}} (GeV/#it{c})");
1508
1509   fhDeltaEtaCharged  = new TH2F
1510   ("hDeltaEtaCharged","#eta_{trigger} - #eta_{h^{#pm}} vs #it{p}_{T trigger}",
1511    nptbins,ptmin,ptmax,ndeltaetabins,deltaetamin,deltaetamax);
1512   fhDeltaEtaCharged->SetYTitle("#Delta #eta");
1513   fhDeltaEtaCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1514
1515   fhDeltaEtaChargedPtA3GeV  = new TH2F
1516   ("hDeltaEtaChargedPtA3GeV","#eta_{trigger} - #eta_{h^{#pm}} vs #it{p}_{T trigger}, #it{p}_{TA}>3 GeV/#it{c}",
1517    nptbins,ptmin,ptmax,ndeltaetabins,deltaetamin,deltaetamax);
1518   fhDeltaEtaChargedPtA3GeV->SetYTitle("#Delta #eta");
1519   fhDeltaEtaChargedPtA3GeV->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1520
1521   fhXECharged  =
1522   new TH2F("hXECharged","#it{x}_{#it{E}} for charged tracks",
1523            nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1524   fhXECharged->SetYTitle("#it{x}_{#it{E}}");
1525   fhXECharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1526
1527   fhXECharged_Cone2  =
1528   new TH2F("hXECharged_Cone2","#it{x}_{#it{E}} for charged tracks in cone 2 (5#pi/6-7#pi/6)",
1529            nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1530   fhXECharged_Cone2->SetYTitle("#it{x}_{#it{E}}");
1531   fhXECharged_Cone2->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1532
1533   fhXEPosCharged  =
1534   new TH2F("hXEPositiveCharged","#it{x}_{#it{E}} for positive charged tracks",
1535            nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1536   fhXEPosCharged->SetYTitle("#it{x}_{#it{E}}");
1537   fhXEPosCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1538
1539   fhXENegCharged  =
1540   new TH2F("hXENegativeCharged","#it{x}_{#it{E}} for negative charged tracks",
1541            nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1542   fhXENegCharged->SetYTitle("#it{x}_{#it{E}}");
1543   fhXENegCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1544
1545   fhPtHbpXECharged  =
1546   new TH2F("hHbpXECharged","#xi = ln(1/#it{x}_{#it{E}}) with charged hadrons",
1547            nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
1548   fhPtHbpXECharged->SetYTitle("ln(1/#it{x}_{#it{E}})");
1549   fhPtHbpXECharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1550
1551   fhPtHbpXECharged_Cone2  =
1552   new TH2F("hHbpXECharged_Cone2","#xi = ln(1/#it{x}_{#it{E}}) with charged hadrons in cone 2 (5#pi/6-7#pi/6)",
1553            nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
1554   fhPtHbpXECharged_Cone2->SetYTitle("ln(1/#it{x}_{#it{E}})");
1555   fhPtHbpXECharged_Cone2->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1556
1557   fhZTCharged  =
1558   new TH2F("hZTCharged","#it{z}_{T} for charged tracks",
1559            nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1560   fhZTCharged->SetYTitle("#it{z}_{T}");
1561   fhZTCharged->SetXTitle("#it{p}_{T trigger}");
1562
1563   fhZTPosCharged  =
1564   new TH2F("hZTPositiveCharged","#it{z}_{T} for positive charged tracks",
1565            nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1566   fhZTPosCharged->SetYTitle("#it{z}_{T}");
1567   fhZTPosCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1568
1569   fhZTNegCharged  =
1570   new TH2F("hZTNegativeCharged","#it{z}_{T} for negative charged tracks",
1571            nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1572   fhZTNegCharged->SetYTitle("#it{z}_{T}");
1573   fhZTNegCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1574
1575   fhPtHbpZTCharged  =
1576   new TH2F("hHbpZTCharged","#xi = ln(1/#it{z}_{T}) with charged hadrons",
1577            nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
1578   fhPtHbpZTCharged->SetYTitle("ln(1/#it{z}_{T})");
1579   fhPtHbpZTCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1580
1581   fhPtTrigPout  =
1582   new TH2F("hPtTrigPout","Pout with triggers",
1583            nptbins,ptmin,ptmax,nptbins,-1.*ptmax/2.,ptmax/2.);
1584   fhPtTrigPout->SetYTitle("#it{p}_{out} (GeV/#it{c})");
1585   fhPtTrigPout->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1586
1587   fhPtTrigCharged  =
1588   new TH2F("hPtTrigCharged","trigger and charged tracks pt distribution",
1589            nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
1590   fhPtTrigCharged->SetYTitle("#it{p}_{T h^{#pm}} (GeV/#it{c})");
1591   fhPtTrigCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1592
1593   outputContainer->Add(fhDeltaPhiDeltaEtaCharged);
1594   outputContainer->Add(fhDeltaPhiDeltaEtaChargedPtA3GeV);
1595   outputContainer->Add(fhPhiCharged) ;
1596   outputContainer->Add(fhEtaCharged) ;
1597   outputContainer->Add(fhDeltaPhiCharged) ;
1598   outputContainer->Add(fhDeltaPhiChargedPtA3GeV) ;
1599   outputContainer->Add(fhDeltaEtaCharged) ;
1600   outputContainer->Add(fhDeltaEtaChargedPtA3GeV) ;
1601   outputContainer->Add(fhDeltaPhiChargedPt) ;
1602
1603   outputContainer->Add(fhXECharged) ;
1604   outputContainer->Add(fhXECharged_Cone2) ;
1605
1606   if(IsDataMC())
1607   {
1608     for(Int_t i=0; i < fgkNmcTypes; i++)
1609     {
1610
1611       fhDeltaPhiChargedMC[i]  = new TH2F(Form("hDeltaPhiCharged_MC%s",nameMC[i].Data()),
1612                                          Form("#Delta #phi for charged tracks, trigger origin is %s",nameMC[i].Data()),
1613                                          nptbins,ptmin,ptmax,ndeltaphibins ,deltaphimin,deltaphimax);
1614       fhDeltaPhiChargedMC[i]->SetYTitle("#it{x}_{#it{E}}");
1615       fhDeltaPhiChargedMC[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1616       outputContainer->Add(fhDeltaPhiChargedMC[i]) ;
1617
1618       fhXEChargedMC[i]  = new TH2F(Form("hXECharged_MC%s",nameMC[i].Data()),
1619                                    Form("#it{x}_{#it{E}} for charged tracks, trigger origin is %s",nameMC[i].Data()),
1620                                    nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1621       fhXEChargedMC[i]->SetYTitle("#it{x}_{#it{E}}");
1622       fhXEChargedMC[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1623       outputContainer->Add(fhXEChargedMC[i]) ;
1624     }
1625   }
1626
1627   outputContainer->Add(fhXEPosCharged) ;
1628   outputContainer->Add(fhXENegCharged) ;
1629   outputContainer->Add(fhPtHbpXECharged) ;
1630   outputContainer->Add(fhPtHbpXECharged_Cone2) ;
1631
1632   outputContainer->Add(fhZTCharged) ;
1633   outputContainer->Add(fhZTPosCharged) ;
1634   outputContainer->Add(fhZTNegCharged) ;
1635   outputContainer->Add(fhPtHbpZTCharged) ;
1636
1637   outputContainer->Add(fhPtTrigPout) ;
1638   outputContainer->Add(fhPtTrigCharged) ;
1639
1640   TString right = "";
1641   if(fMakeSeveralUE) right = "Right";
1642
1643   fhUePart  =  new TH1F("hUePart","UE particles distribution vs pt trig",
1644                         nptbins,ptmin,ptmax);
1645   fhUePart->SetYTitle("dNch");
1646   fhUePart->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1647
1648   fhDeltaPhiUeChargedPt  = new TH2F
1649   (Form("hDeltaPhiUe%sChargedPt",right.Data()),"#phi_{trigger} - #phi_{#Ueh^{#pm}} vs #it{p}_{T Ueh^{#pm}}",
1650    nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1651   fhDeltaPhiUeChargedPt->SetYTitle("#Delta #phi (rad)");
1652   fhDeltaPhiUeChargedPt->SetXTitle("#it{p}_{T h^{#pm}} (GeV/#it{c})");
1653
1654   fhXEUeCharged  =
1655   new TH2F(Form("hXEUeCharged%s",right.Data()),"#it{x}_{#it{E}} for Underlying Event",
1656            nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1657   fhXEUeCharged->SetYTitle("#it{x}_{#it{E}}");
1658   fhXEUeCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1659
1660   fhPtHbpXEUeCharged  =
1661   new TH2F(Form("hHbpXEUeCharged%s",right.Data()),"#xi = ln(1/#it{x}_{#it{E}}) for Underlying Event",
1662            nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
1663   fhPtHbpXEUeCharged->SetYTitle("ln(1/#it{x}_{#it{E}})");
1664   fhPtHbpXEUeCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1665
1666   fhZTUeCharged  =
1667   new TH2F(Form("hZTUeCharged%s",right.Data()),"#it{z}_{T} for Underlying Event",
1668            nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1669   fhZTUeCharged->SetYTitle("#it{z}_{T}");
1670   fhZTUeCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1671
1672   fhPtHbpZTUeCharged  =
1673   new TH2F(Form("hHbpZTUeCharged%s",right.Data()),"#xi = ln(1/#it{z}_{T}) for Underlying Event",
1674            nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
1675   fhPtHbpZTUeCharged->SetYTitle("ln(1/#it{x}_{#it{E}})");
1676   fhPtHbpZTUeCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1677
1678   outputContainer->Add(fhUePart);
1679   outputContainer->Add(fhDeltaPhiUeChargedPt) ;
1680   outputContainer->Add(fhXEUeCharged) ;
1681   outputContainer->Add(fhPtHbpXEUeCharged) ;
1682   outputContainer->Add(fhZTUeCharged) ;
1683   outputContainer->Add(fhPtHbpZTUeCharged) ;
1684
1685   if(fMakeSeveralUE)
1686   {
1687     fhDeltaPhiUeLeftCharged  = new TH2F
1688     ("hDeltaPhiUeLeftChargedPt","#phi_{trigger} - #phi_{#Ueh^{#pm}} vs #it{p}_{T Ueh^{#pm}} with UE left side range of trigger particles",
1689      nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1690     fhDeltaPhiUeLeftCharged->SetYTitle("#Delta #phi (rad)");
1691     fhDeltaPhiUeLeftCharged->SetXTitle("#it{p}_{T h^{#pm}} (GeV/#it{c})");
1692     outputContainer->Add(fhDeltaPhiUeLeftCharged) ;
1693
1694     fhDeltaPhiUeLeftUpCharged  = new TH2F
1695     ("hDeltaPhiUeLeftUpChargedPt","#phi_{trigger} - #phi_{#Ueh^{#pm}} vs #it{p}_{T Ueh^{#pm}} with UE left Up side range of trigger particles",
1696      nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1697     fhDeltaPhiUeLeftUpCharged->SetYTitle("#Delta #phi (rad)");
1698     fhDeltaPhiUeLeftUpCharged->SetXTitle("#it{p}_{T h^{#pm}} (GeV/#it{c})");
1699     outputContainer->Add(fhDeltaPhiUeLeftUpCharged) ;
1700
1701     fhDeltaPhiUeRightUpCharged  = new TH2F
1702     ("hDeltaPhiUeRightUpChargedPt","#phi_{trigger} - #phi_{#Ueh^{#pm}} vs #it{p}_{T Ueh^{#pm}} with UE right Up side range of trigger particles",
1703      nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1704     fhDeltaPhiUeRightUpCharged->SetYTitle("#Delta #phi (rad)");
1705     fhDeltaPhiUeRightUpCharged->SetXTitle("#it{p}_{T h^{#pm}} (GeV/#it{c})");
1706     outputContainer->Add(fhDeltaPhiUeRightUpCharged) ;
1707
1708     fhDeltaPhiUeLeftDownCharged  = new TH2F
1709     ("hDeltaPhiUeLeftDownChargedPt","#phi_{trigger} - #phi_{#Ueh^{#pm}} vs #it{p}_{T Ueh^{#pm}} with UE left Down side range of trigger particles",
1710      nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1711     fhDeltaPhiUeLeftDownCharged->SetYTitle("#Delta #phi (rad)");
1712     fhDeltaPhiUeLeftDownCharged->SetXTitle("#it{p}_{T h^{#pm}} (GeV/#it{c})");
1713     outputContainer->Add(fhDeltaPhiUeLeftDownCharged) ;
1714
1715     fhDeltaPhiUeRightDownCharged  = new TH2F
1716     ("hDeltaPhiUeRightDownChargedPt","#phi_{trigger} - #phi_{#Ueh^{#pm}} vs #it{p}_{T Ueh^{#pm}} with UE right Down side range of trigger particles",
1717      nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1718     fhDeltaPhiUeRightDownCharged->SetYTitle("#Delta #phi (rad)");
1719     fhDeltaPhiUeRightDownCharged->SetXTitle("#it{p}_{T h^{#pm}} (GeV/#it{c})");
1720     outputContainer->Add(fhDeltaPhiUeRightDownCharged) ;
1721
1722     fhXEUeLeftCharged  =
1723     new TH2F("hXEUeChargedLeft","#it{x}_{#it{E}} with UE left side of trigger",
1724              nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1725     fhXEUeLeftCharged->SetYTitle("#it{x}_{#it{E} Ueh^{#pm}}");
1726     fhXEUeLeftCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1727     outputContainer->Add(fhXEUeLeftCharged) ;
1728
1729     fhXEUeLeftUpCharged  =
1730     new TH2F("hXEUeChargedLeftUp","#it{x}_{#it{E}} with UE left Up side of trigger",
1731              nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1732     fhXEUeLeftUpCharged->SetYTitle("#it{x}_{#it{E} Ueh^{#pm}}");
1733     fhXEUeLeftUpCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1734     outputContainer->Add(fhXEUeLeftUpCharged) ;
1735
1736     fhXEUeRightUpCharged  =
1737     new TH2F("hXEUeChargedRightUp","#it{x}_{#it{E} h^{#pm}} with UE right Up side of trigger",
1738              nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1739     fhXEUeRightUpCharged->SetYTitle("#it{z}_{trigger Ueh^{#pm}}");
1740     fhXEUeRightUpCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1741     outputContainer->Add(fhXEUeRightUpCharged) ;
1742
1743     fhXEUeLeftDownCharged  =
1744     new TH2F("hXEUeChargedLeftDown","#it{x}_{#it{E}} with UE left Down side of trigger",
1745              nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1746     fhXEUeLeftDownCharged->SetYTitle("#it{x}_{#it{E} Ueh^{#pm}}");
1747     fhXEUeLeftDownCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1748     outputContainer->Add(fhXEUeLeftDownCharged) ;
1749
1750     fhXEUeRightDownCharged  =
1751     new TH2F("hXEUeChargedRightDown","#it{x}_{#it{E} h^{#pm}} with UE right Down side of trigger",
1752              nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1753     fhXEUeRightDownCharged->SetYTitle("#it{z}_{trigger Ueh^{#pm}}");
1754     fhXEUeRightDownCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1755     outputContainer->Add(fhXEUeRightDownCharged) ;
1756
1757     fhPtHbpXEUeLeftCharged  =
1758     new TH2F("hHbpXEUeChargedLeft","#xi = ln(1/#it{x}_{#it{E}}) with charged UE left side of trigger",
1759              nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
1760     fhPtHbpXEUeLeftCharged->SetYTitle("ln(1/#it{x}_{#it{E}})");
1761     fhPtHbpXEUeLeftCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1762     outputContainer->Add(fhPtHbpXEUeLeftCharged) ;
1763
1764     fhZTUeLeftCharged  =
1765     new TH2F("hZTUeChargedLeft","#it{z}_{trigger h^{#pm}} = #it{p}_{T Ueh^{#pm}} / #it{p}_{T trigger} with UE left side of trigger",
1766              nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1767     fhZTUeLeftCharged->SetYTitle("#it{z}_{trigger Ueh^{#pm}}");
1768     fhZTUeLeftCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1769     outputContainer->Add(fhZTUeLeftCharged) ;
1770
1771     fhPtHbpZTUeLeftCharged  =
1772     new TH2F("hHbpZTUeChargedLeft","#xi = ln(1/#it{z}_{T}) with charged UE left side of trigger",
1773              nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
1774     fhPtHbpZTUeLeftCharged->SetYTitle("ln(1/#it{z}_{T})");
1775     fhPtHbpZTUeLeftCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1776     outputContainer->Add(fhPtHbpZTUeLeftCharged) ;
1777   }
1778
1779   if(IsPileUpAnalysisOn())
1780   {
1781     fhDeltaPhiChargedOtherBC  = new TH2F
1782     ("hDeltaPhiChargedOtherBC","#phi_{trigger} - #phi_{h^{#pm}} vs #it{p}_{T trigger}, track BC!=0",
1783      nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1784     fhDeltaPhiChargedOtherBC->SetYTitle("#Delta #phi (rad)");
1785     fhDeltaPhiChargedOtherBC->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1786
1787     fhDeltaPhiChargedPtA3GeVOtherBC  = new TH2F
1788     ("hDeltaPhiChargedPtA3GeVOtherBC","#phi_{trigger} - #phi_{h^{#pm}} vs #it{p}_{T trigger}, #it{p}_{TA}>3 GeV/#it{c}, track BC!=0",
1789      nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1790     fhDeltaPhiChargedPtA3GeVOtherBC->SetYTitle("#Delta #phi (rad)");
1791     fhDeltaPhiChargedPtA3GeVOtherBC->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1792
1793     fhPtTrigChargedOtherBC  =
1794     new TH2F("hPtTrigChargedOtherBC","trigger and charged tracks pt distribution, track BC!=0",
1795              nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
1796     fhPtTrigChargedOtherBC->SetYTitle("#it{p}_{T h^{#pm}} (GeV/#it{c})");
1797     fhPtTrigChargedOtherBC->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1798
1799     fhXEChargedOtherBC  =
1800     new TH2F("hXEChargedOtherBC","#it{x}_{#it{E}} for charged tracks, track BC!=0",
1801              nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1802     fhXEChargedOtherBC->SetYTitle("#it{x}_{#it{E}}");
1803     fhXEChargedOtherBC->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1804
1805     fhXEUeChargedOtherBC  =
1806     new TH2F("hXEUeChargedOtherBC","#it{x}_{#it{E}} for Underlying Event, track BC!=0",
1807              nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1808     fhXEUeChargedOtherBC->SetYTitle("#it{x}_{#it{E}}");
1809     fhXEUeChargedOtherBC->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1810
1811     fhZTChargedOtherBC  =
1812     new TH2F("hZTChargedOtherBC","#it{z}_{T} for charged tracks, track BC!=0",
1813              nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1814     fhZTChargedOtherBC->SetYTitle("#it{z}_{T}");
1815     fhZTChargedOtherBC->SetXTitle("#it{p}_{T trigger}");
1816
1817     fhZTUeChargedOtherBC  =
1818     new TH2F("hZTUeChargedOtherBC","#it{z}_{T} for Underlying Event, track BC!=0",
1819              nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1820     fhZTUeChargedOtherBC->SetYTitle("#it{z}_{T}");
1821     fhZTUeChargedOtherBC->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1822
1823     outputContainer->Add(fhDeltaPhiChargedOtherBC) ;
1824     outputContainer->Add(fhDeltaPhiChargedPtA3GeVOtherBC) ;
1825     outputContainer->Add(fhXEChargedOtherBC) ;
1826     outputContainer->Add(fhXEUeChargedOtherBC) ;
1827     outputContainer->Add(fhZTChargedOtherBC) ;
1828     outputContainer->Add(fhZTUeChargedOtherBC) ;
1829     outputContainer->Add(fhPtTrigChargedOtherBC) ;
1830
1831     fhDeltaPhiChargedBC0  = new TH2F
1832     ("hDeltaPhiChargedBC0","#phi_{trigger} - #phi_{h^{#pm}} vs #it{p}_{T trigger}, track BC==0",
1833      nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1834     fhDeltaPhiChargedBC0->SetYTitle("#Delta #phi (rad)");
1835     fhDeltaPhiChargedBC0->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1836
1837     fhDeltaPhiChargedPtA3GeVBC0  = new TH2F
1838     ("hDeltaPhiChargedPtA3GeVBC0","#phi_{trigger} - #phi_{h^{#pm}} vs #it{p}_{T trigger}, #it{p}_{TA}>3 GeV/#it{c}, track BC==0",
1839      nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1840     fhDeltaPhiChargedPtA3GeVBC0->SetYTitle("#Delta #phi (rad)");
1841     fhDeltaPhiChargedPtA3GeVBC0->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1842
1843     fhPtTrigChargedBC0  =
1844     new TH2F("hPtTrigChargedBC0","trigger and charged tracks pt distribution, track BC==0",
1845              nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
1846     fhPtTrigChargedBC0->SetYTitle("#it{p}_{T h^{#pm}} (GeV/#it{c})");
1847     fhPtTrigChargedBC0->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1848
1849     fhXEChargedBC0  =
1850     new TH2F("hXEChargedBC0","#it{x}_{#it{E}} for charged tracks, track BC==0",
1851              nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1852     fhXEChargedBC0->SetYTitle("#it{x}_{#it{E}}");
1853     fhXEChargedBC0->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1854
1855     fhXEUeChargedBC0  =
1856     new TH2F("hXEUeChargedBC0","#it{x}_{#it{E}} for Underlying Event, track BC==0",
1857              nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1858     fhXEUeChargedBC0->SetYTitle("#it{x}_{#it{E}}");
1859     fhXEUeChargedBC0->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1860
1861     fhZTChargedBC0  =
1862     new TH2F("hZTChargedBC0","#it{z}_{T} for charged tracks, track BC==0",
1863              nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1864     fhZTChargedBC0->SetYTitle("#it{z}_{T}");
1865     fhZTChargedBC0->SetXTitle("#it{p}_{T trigger}");
1866
1867     fhZTUeChargedBC0  =
1868     new TH2F("hZTUeChargedBC0","#it{z}_{T} for Underlying Event, track BC==0",
1869              nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1870     fhZTUeChargedBC0->SetYTitle("#it{z}_{T}");
1871     fhZTUeChargedBC0->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1872
1873     outputContainer->Add(fhDeltaPhiChargedBC0) ;
1874     outputContainer->Add(fhDeltaPhiChargedPtA3GeVBC0) ;
1875     outputContainer->Add(fhXEChargedBC0) ;
1876     outputContainer->Add(fhXEUeChargedBC0) ;
1877     outputContainer->Add(fhZTChargedBC0) ;
1878     outputContainer->Add(fhZTUeChargedBC0) ;
1879     outputContainer->Add(fhPtTrigChargedBC0) ;
1880
1881     fhPtTriggerVtxBC0  = new TH1F("hPtTriggerVtxBC0","#it{p}_{T} distribution of trigger particles", nptbins,ptmin,ptmax);
1882     fhPtTriggerVtxBC0->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1883
1884     fhDeltaPhiChargedVtxBC0  = new TH2F
1885     ("hDeltaPhiChargedVtxBC0","#phi_{trigger} - #phi_{h^{#pm}} vs #it{p}_{T trigger}, track BC==0",
1886      nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1887     fhDeltaPhiChargedVtxBC0->SetYTitle("#Delta #phi (rad)");
1888     fhDeltaPhiChargedVtxBC0->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1889
1890     fhDeltaPhiChargedPtA3GeVVtxBC0  = new TH2F
1891     ("hDeltaPhiChargedPtA3GeVVtxBC0","#phi_{trigger} - #phi_{h^{#pm}} vs #it{p}_{T trigger}, #it{p}_{TA}>3 GeV/#it{c}, track BC==0",
1892      nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1893     fhDeltaPhiChargedPtA3GeVVtxBC0->SetYTitle("#Delta #phi (rad)");
1894     fhDeltaPhiChargedPtA3GeVVtxBC0->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1895
1896     fhPtTrigChargedVtxBC0  =
1897     new TH2F("hPtTrigChargedVtxBC0","trigger and charged tracks pt distribution, track BC==0",
1898              nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
1899     fhPtTrigChargedVtxBC0->SetYTitle("#it{p}_{T h^{#pm}} (GeV/#it{c})");
1900     fhPtTrigChargedVtxBC0->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1901
1902     fhXEChargedVtxBC0  =
1903     new TH2F("hXEChargedVtxBC0","#it{x}_{#it{E}} for charged tracks, track BC==0",
1904              nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1905     fhXEChargedVtxBC0->SetYTitle("#it{x}_{#it{E}}");
1906     fhXEChargedVtxBC0->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1907
1908     fhXEUeChargedVtxBC0  =
1909     new TH2F("hXEUeChargedVtxBC0","#it{x}_{#it{E}} for Underlying Event, track BC==0",
1910              nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1911     fhXEUeChargedVtxBC0->SetYTitle("#it{x}_{#it{E}}");
1912     fhXEUeChargedVtxBC0->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1913
1914     fhZTChargedVtxBC0  =
1915     new TH2F("hZTChargedVtxBC0","#it{z}_{T} for charged tracks, track BC==0",
1916              nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1917     fhZTChargedVtxBC0->SetYTitle("#it{z}_{T}");
1918     fhZTChargedVtxBC0->SetXTitle("#it{p}_{T trigger}");
1919
1920     fhZTUeChargedVtxBC0  =
1921     new TH2F("hZTUeChargedVtxBC0","#it{z}_{T} for Underlying Event, track BC==0",
1922              nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1923     fhZTUeChargedVtxBC0->SetYTitle("#it{z}_{T}");
1924     fhZTUeChargedVtxBC0->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1925
1926     outputContainer->Add(fhPtTriggerVtxBC0);
1927     outputContainer->Add(fhDeltaPhiChargedVtxBC0) ;
1928     outputContainer->Add(fhDeltaPhiChargedPtA3GeVVtxBC0) ;
1929     outputContainer->Add(fhXEChargedVtxBC0) ;
1930     outputContainer->Add(fhXEUeChargedVtxBC0) ;
1931     outputContainer->Add(fhZTChargedVtxBC0) ;
1932     outputContainer->Add(fhZTUeChargedVtxBC0) ;
1933     outputContainer->Add(fhPtTrigChargedVtxBC0) ;
1934
1935     for(Int_t i = 0 ; i < 7 ; i++)
1936     {
1937       fhPtTriggerPileUp[i]  = new TH1F(Form("hPtTriggerPileUp%s",pileUpName[i].Data()),
1938                                        Form("#it{p}_{T} distribution of trigger particles, %s Pile-Up event",pileUpName[i].Data()), nptbins,ptmin,ptmax);
1939       fhPtTriggerPileUp[i]->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1940       outputContainer->Add(fhPtTriggerPileUp[i]);
1941
1942       fhDeltaPhiChargedPileUp[i]  = new TH2F(Form("hDeltaPhiChargedPileUp%s",pileUpName[i].Data()),
1943                                              Form("#phi_{trigger} - #phi_{h^{#pm}} vs #it{p}_{T trigger}, %s Pile-Up event",pileUpName[i].Data()),
1944                                              nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1945       fhDeltaPhiChargedPileUp[i]->SetYTitle("#Delta #phi (rad)");
1946       fhDeltaPhiChargedPileUp[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1947       outputContainer->Add(fhDeltaPhiChargedPileUp[i]) ;
1948
1949       fhDeltaPhiChargedPtA3GeVPileUp[i]  = new TH2F(Form("hDeltaPhiChargedPtA3GeVPileUp%s",pileUpName[i].Data()),
1950                                                     Form("#phi_{trigger} - #phi_{h^{#pm}} vs #it{p}_{T trigger}, #it{p}_{TA}>3 GeV/#it{c}, %s Pile-Up event",pileUpName[i].Data()),
1951                                                     nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1952       fhDeltaPhiChargedPtA3GeVPileUp[i]->SetYTitle("#Delta #phi (rad)");
1953       fhDeltaPhiChargedPtA3GeVPileUp[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1954       outputContainer->Add(fhDeltaPhiChargedPtA3GeVPileUp[i]) ;
1955
1956       fhDeltaEtaChargedPileUp[i]  = new TH2F(Form("hDeltaEtaChargedPileUp%s",pileUpName[i].Data()),
1957                                              Form("#eta_{trigger} - #eta_{h^{#pm}} vs #it{p}_{T trigger}, %s Pile-Up event",pileUpName[i].Data()),
1958                                              nptbins,ptmin,ptmax,ndeltaetabins,deltaetamin,deltaetamax);
1959       fhDeltaEtaChargedPileUp[i]->SetYTitle("#Delta #eta");
1960       fhDeltaEtaChargedPileUp[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1961       outputContainer->Add(fhDeltaEtaChargedPileUp[i]) ;
1962
1963       fhDeltaEtaChargedPtA3GeVPileUp[i]  = new TH2F(Form("hDeltaEtaChargedPtA3GeVPileUp%s",pileUpName[i].Data()),
1964                                                     Form("#eta_{trigger} - #eta_{h^{#pm}} vs #it{p}_{T trigger}, #it{p}_{TA}>3 GeV/#it{c}, %s Pile-Up event",pileUpName[i].Data()),
1965                                                     nptbins,ptmin,ptmax,ndeltaetabins,deltaetamin,deltaetamax);
1966       fhDeltaEtaChargedPtA3GeVPileUp[i]->SetYTitle("#Delta #eta");
1967       fhDeltaEtaChargedPtA3GeVPileUp[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1968       outputContainer->Add(fhDeltaEtaChargedPtA3GeVPileUp[i]) ;
1969
1970       fhXEChargedPileUp[i]  = new TH2F(Form("hXEChargedPileUp%s",pileUpName[i].Data()),
1971                                        Form("#it{x}_{#it{E}} for charged tracks, %s Pile-Up event",pileUpName[i].Data()),
1972                                        nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1973       fhXEChargedPileUp[i]->SetYTitle("#it{x}_{#it{E}}");
1974       fhXEChargedPileUp[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1975       outputContainer->Add(fhXEChargedPileUp[i]) ;
1976
1977       fhXEUeChargedPileUp[i]  = new TH2F(Form("hXEUeChargedPileUp%s",pileUpName[i].Data()),
1978                                          Form("#it{x}_{#it{E}} for Underlying Event, %s Pile-Up event",pileUpName[i].Data()),
1979                                          nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1980       fhXEUeChargedPileUp[i]->SetYTitle("#it{x}_{#it{E}}");
1981       fhXEUeChargedPileUp[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1982       outputContainer->Add(fhXEUeChargedPileUp[i]) ;
1983
1984       fhZTChargedPileUp[i]  = new TH2F(Form("hZTChargedPileUp%s",pileUpName[i].Data()),
1985                                        Form("#it{z}_{T} for charged tracks, %s Pile-Up event",pileUpName[i].Data()),
1986                                        nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1987       fhZTChargedPileUp[i]->SetYTitle("#it{z}_{T}");
1988       fhZTChargedPileUp[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1989       outputContainer->Add(fhZTChargedPileUp[i]) ;
1990
1991       fhZTUeChargedPileUp[i]  = new TH2F(Form("hZTUeChargedPileUp%s",pileUpName[i].Data()),
1992                                          Form("#it{z}_{T} for Underlying Event, %s Pile-Up event",pileUpName[i].Data()),
1993                                          nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
1994       fhZTUeChargedPileUp[i]->SetYTitle("#it{z}_{T}");
1995       fhZTUeChargedPileUp[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1996       outputContainer->Add(fhZTUeChargedPileUp[i]) ;
1997
1998       fhPtTrigChargedPileUp[i]  = new TH2F(Form("hPtTrigChargedPileUp%s",pileUpName[i].Data()),
1999                                            Form("trigger and charged tracks pt distribution, %s Pile-Up event",pileUpName[i].Data()),
2000                                            nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2001       fhPtTrigChargedPileUp[i]->SetYTitle("#it{p}_{T h^{#pm}} (GeV/#it{c})");
2002       fhPtTrigChargedPileUp[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2003       outputContainer->Add(fhPtTrigChargedPileUp[i]) ;
2004
2005     }
2006   }
2007
2008   if(IsHighMultiplicityAnalysisOn())
2009   {
2010     Int_t nMultiBins = GetNCentrBin();
2011     fhDeltaPhiChargedMult = new TH2F*[nMultiBins] ;
2012     fhDeltaEtaChargedMult = new TH2F*[nMultiBins] ;
2013     fhXEMult              = new TH2F*[nMultiBins] ;
2014     fhXEUeMult            = new TH2F*[nMultiBins] ;
2015     fhZTMult              = new TH2F*[nMultiBins] ;
2016     fhZTUeMult            = new TH2F*[nMultiBins] ;
2017
2018     for(Int_t im=0; im<nMultiBins; im++)
2019     {
2020       fhDeltaPhiChargedMult[im]  = new TH2F
2021       (Form("hDeltaPhiCharged_Mult%d",im),Form("#Delta #phi charged Mult bin %d",im), nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2022       fhDeltaPhiChargedMult[im]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2023       fhDeltaPhiChargedMult[im]->SetYTitle("#Delta #phi (rad)");
2024
2025       fhDeltaEtaChargedMult[im]  = new TH2F
2026       (Form("hDeltaEtaCharged_Mult%d",im),Form("#Delta #eta charged Mult bin %d",im), nptbins,ptmin,ptmax, ndeltaetabins ,deltaetamin,deltaetamax);
2027       fhDeltaEtaChargedMult[im]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2028       fhDeltaEtaChargedMult[im]->SetYTitle("#Delta #eta");
2029
2030       fhXEMult[im]  = new TH2F
2031       (Form("hXECharged_Mult%d",im),Form("#it{x}_{E} charged Mult bin %d",im), nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2032       fhXEMult[im]->SetYTitle("#it{x}_{E}");
2033       fhXEMult[im]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2034
2035       fhXEUeMult[im]  = new TH2F
2036       (Form("hXEUeCharged_Mult%d",im),Form("#it{x}_{E} UE charged Mult bin %d",im), nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2037       fhXEUeMult[im]->SetYTitle("#it{x}_{E}");
2038       fhXEUeMult[im]->SetXTitle("#it{p}_{T trigger}(GeV/#it{c})");
2039
2040       fhZTMult[im]  = new TH2F
2041       (Form("hZTCharged_Mult%d",im),Form("#it{z}_{T} charged  Mult bin %d",im), nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2042       fhZTMult[im]->SetYTitle("#it{z}_{T}");
2043       fhZTMult[im]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2044
2045       fhZTUeMult[im]  = new TH2F
2046       (Form("hZTUeCharged_Mult%d",im),Form("#it{z}_{T} UE charged  Mult bin %d",im), nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2047       fhZTUeMult[im]->SetYTitle("#it{z}_{T}");
2048       fhZTUeMult[im]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2049
2050       outputContainer->Add(fhDeltaPhiChargedMult[im]) ;
2051       outputContainer->Add(fhDeltaEtaChargedMult[im]) ;
2052       outputContainer->Add(fhXEMult  [im]);
2053       outputContainer->Add(fhXEUeMult[im]);
2054       outputContainer->Add(fhZTMult  [im]);
2055       outputContainer->Add(fhZTUeMult[im]);
2056     }
2057   }
2058
2059   if(fFillBradHisto)
2060   {
2061     fhAssocPtBkg        = new TH2F("hAssocPtBkg", " Trigger #it{p}_{T} vs associated hadron #it{p}_{T} from background",
2062                                    nptbins, ptmin, ptmax,nptbins,ptmin,ptmax);
2063     fhAssocPtBkg->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2064     fhAssocPtBkg->SetYTitle("#it{p}_{T associated} (GeV/#it{c})");
2065     outputContainer->Add(fhAssocPtBkg) ;
2066
2067     fhDeltaPhiBrad = new TH2F("hDeltaPhiBrad","atan2(sin(#Delta #phi), cos(#Delta #phi))/#pi vs #it{p}_{T trigger} ",
2068                               nptbins, ptmin, ptmax,288, -1.0/3.0, 5.0/3.0);
2069     fhDeltaPhiBrad->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2070     fhDeltaPhiBrad->SetYTitle("atan2(sin(#Delta #phi), cos(#Delta #phi))/#pi");
2071     outputContainer->Add(fhDeltaPhiBrad) ;
2072   }
2073
2074   fhDeltaPhiDeltaEtaAssocPtBin = new TH2F*[fNAssocPtBins*nz];
2075   fhDeltaPhiAssocPtBin       = new TH2F*[fNAssocPtBins*nz];
2076   fhDeltaPhiAssocPtBinDEta08 = new TH2F*[fNAssocPtBins*nz];
2077   fhDeltaPhiAssocPtBinDEta0  = new TH2F*[fNAssocPtBins*nz];
2078   if(fFillMomImbalancePtAssocBinsHisto)
2079   {
2080     fhXEAssocPtBin           = new TH2F*[fNAssocPtBins*nz];
2081     fhZTAssocPtBin           = new TH2F*[fNAssocPtBins*nz];
2082   }
2083
2084   if(fCorrelVzBin)
2085   {
2086     fhXEVZ = new TH2F*[nz];
2087     fhZTVZ = new TH2F*[nz];
2088   }
2089
2090   if(fFillBradHisto)
2091     fhDeltaPhiBradAssocPtBin = new TH2F*[fNAssocPtBins*nz];
2092
2093
2094   fhDeltaPhiAssocPtBin       = new TH2F*[fNAssocPtBins*nz];
2095   if(fFillEtaGapsHisto)fhDeltaPhiAssocPtBinDEta08       = new TH2F*[fNAssocPtBins*nz];
2096   if(fDecayTrigger)    fhDeltaPhiDecayChargedAssocPtBin = new TH2F*[fNAssocPtBins*nz];
2097
2098
2099   if(fHMPIDCorrelation)
2100   {
2101     fhDeltaPhiAssocPtBinHMPID   = new TH2F*[fNAssocPtBins*nz];
2102     fhDeltaPhiAssocPtBinHMPIDAcc= new TH2F*[fNAssocPtBins*nz];
2103   }
2104
2105   for(Int_t i = 0 ; i < fNAssocPtBins ; i++)
2106   {
2107     for(Int_t z = 0 ; z < nz ; z++)
2108     {
2109       Int_t bin = i*nz+z;
2110
2111       if(fCorrelVzBin)
2112       {
2113         sz = Form("_vz%d",z);
2114         tz = Form(", #it{v}_{#it{z}} bin %d",z);
2115       }
2116
2117       //printf("iAssoc %d, Vz %d, bin %d - sz %s, tz %s \n",i,z,bin,sz.Data(),tz.Data());
2118
2119       fhDeltaPhiDeltaEtaAssocPtBin[bin]  = new TH2F(Form("hDeltaPhiDeltaEtaPtAssocPt%2.1f_%2.1f%s", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],sz.Data()),
2120                                                     Form("#Delta #phi vs #Delta #eta for associated #it{p}_{T} bin [%2.1f,%2.1f]%s", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],tz.Data()),
2121                                                     ndeltaphibins ,deltaphimin,deltaphimax,ndeltaetabins,deltaetamin,deltaetamax);
2122       fhDeltaPhiDeltaEtaAssocPtBin[bin]->SetXTitle("#Delta #phi (rad)");
2123       fhDeltaPhiDeltaEtaAssocPtBin[bin]->SetYTitle("#Delta #eta");
2124
2125       fhDeltaPhiAssocPtBin[bin] = new TH2F(Form("hDeltaPhiPtAssocPt%2.1f_%2.1f%s", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],sz.Data()),
2126                                            Form("#Delta #phi vs #it{p}_{T trigger} for associated #it{p}_{T} bin [%2.1f,%2.1f]%s", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],tz.Data()),
2127                                            nptbins, ptmin, ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2128       fhDeltaPhiAssocPtBin[bin]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2129       fhDeltaPhiAssocPtBin[bin]->SetYTitle("#Delta #phi (rad)");
2130
2131       outputContainer->Add(fhDeltaPhiDeltaEtaAssocPtBin[bin]) ;
2132       outputContainer->Add(fhDeltaPhiAssocPtBin[bin]) ;
2133
2134       if(fFillEtaGapsHisto)
2135       {
2136         fhDeltaPhiAssocPtBinDEta08[bin] = new TH2F(Form("hDeltaPhiDeltaEta0.8PtAssocPt%2.1f_%2.1f%s", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],sz.Data()),
2137                                                    Form("#Delta #phi vs #it{p}_{T trigger} for associated #it{p}_{T} bin [%2.1f,%2.1f]%s, for #Delta #eta > 0.8", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],tz.Data()),
2138                                                    nptbins, ptmin, ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2139         fhDeltaPhiAssocPtBinDEta08[bin]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2140         fhDeltaPhiAssocPtBinDEta08[bin]->SetYTitle("#Delta #phi (rad)");
2141
2142         fhDeltaPhiAssocPtBinDEta0[bin] = new TH2F(Form("hDeltaPhiDeltaEta0PtAssocPt%2.1f_%2.1f%s", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],tz.Data()),
2143                                                   Form("#Delta #phi vs #it{p}_{T trigger} for associated #it{p}_{T} bin [%2.1f,%2.1f]%s, for #Delta #eta = 0.", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],tz.Data()),
2144                                                   nptbins, ptmin, ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2145         fhDeltaPhiAssocPtBinDEta0[bin]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2146         fhDeltaPhiAssocPtBinDEta0[bin]->SetYTitle("#Delta #phi (rad)");
2147
2148         outputContainer->Add(fhDeltaPhiAssocPtBinDEta08[bin]) ;
2149         outputContainer->Add(fhDeltaPhiAssocPtBinDEta0[bin]) ;
2150       }
2151
2152       if(fDecayTrigger)
2153       {
2154         fhDeltaPhiDecayChargedAssocPtBin[bin] = new TH2F(Form("hDeltaPhiPtDecayChargedAssocPt%2.1f_%2.1f%s_bit%d", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],sz.Data(),fDecayBits[0]),
2155                                                          Form("#Delta #phi vs #it{p}_{T trigger} tagged as decay for associated #it{p}_{T} bin [%2.1f,%2.1f]%s, Bit %d", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],tz.Data(),fDecayBits[0]),
2156                                                          nptbins, ptmin, ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2157         fhDeltaPhiDecayChargedAssocPtBin[bin]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2158         fhDeltaPhiDecayChargedAssocPtBin[bin]->SetYTitle("#Delta #phi (rad)");
2159
2160         outputContainer->Add(fhDeltaPhiDecayChargedAssocPtBin[bin]) ;
2161       }
2162
2163       if(fFillBradHisto)
2164       {
2165         fhDeltaPhiBradAssocPtBin[bin] = new TH2F(Form("hDeltaPhiBradPtAssocPt%2.1f_%2.1f%s", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],sz.Data()),
2166                                                  Form("atan2(sin(#Delta #phi), cos(#Delta #phi))/#pi vs #it{p}_{T trigger} for associated #it{p}_{T} bin [%2.1f,%2.1f]%s", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],tz.Data()),
2167                                                  nptbins, ptmin, ptmax,288, -1.0/3.0, 5.0/3.0);
2168         fhDeltaPhiBradAssocPtBin[bin]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2169         fhDeltaPhiBradAssocPtBin[bin]->SetYTitle("atan2(sin(#Delta #phi), cos(#Delta #phi))/#pi");
2170         outputContainer->Add(fhDeltaPhiBradAssocPtBin[bin]) ;
2171       }
2172
2173       if(fHMPIDCorrelation)
2174       {
2175         fhDeltaPhiAssocPtBinHMPID[bin] = new TH2F(Form("hDeltaPhiPtAssocPt%2.1f_%2.1f%sHMPID", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],sz.Data()),
2176                                                   Form("#Delta #phi vs #it{p}_{T trigger} for associated #it{p}_{T} bin [%2.1f,%2.1f]%s, with track having HMPID signal", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],tz.Data()),
2177                                                   nptbins, ptmin, ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2178         fhDeltaPhiAssocPtBinHMPID[bin]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})" );
2179         fhDeltaPhiAssocPtBinHMPID[bin]->SetYTitle("#Delta #phi (rad)");
2180
2181         fhDeltaPhiAssocPtBinHMPIDAcc[bin] = new TH2F(Form("hDeltaPhiPtAssocPt%2.1f_%2.1f%sHMPIDAcc", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],sz.Data()),
2182                                                      Form("#Delta #phi vs #it{p}_{T trigger} for associated #it{p}_{T} bin [%2.1f,%2.1f]%s, with track within 5<phi<20 deg", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],tz.Data()),
2183                                                      nptbins, ptmin, ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2184         fhDeltaPhiAssocPtBinHMPIDAcc[bin]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2185         fhDeltaPhiAssocPtBinHMPIDAcc[bin]->SetYTitle("#Delta #phi (rad)");
2186
2187         outputContainer->Add(fhDeltaPhiAssocPtBinHMPID   [bin]) ;
2188         outputContainer->Add(fhDeltaPhiAssocPtBinHMPIDAcc[bin]) ;
2189       }
2190     }
2191   }
2192
2193   if(fFillMomImbalancePtAssocBinsHisto)
2194   {
2195     for(Int_t i = 0 ; i < fNAssocPtBins ; i++)
2196     {
2197       fhXEAssocPtBin[i]       = new TH2F(Form("hXEAssocPtBin%1.f_%1.f", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
2198                                          Form("#it{x}_{#it{E}} vs #it{p}_{T trigger} for associated #it{p}_{T} bin [%2.1f,%2.1f]", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
2199                                          nptbins, ptmin, ptmax,nxeztbins,xeztmin,xeztmax);
2200       fhXEAssocPtBin[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2201       fhXEAssocPtBin[i]->SetYTitle("#it{x}_{#it{E}}");
2202
2203       fhZTAssocPtBin[i]       = new TH2F(Form("hZTAssocPtBin%1.f_%1.f", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
2204                                          Form("#it{z}_{T} vs #it{p}_{T trigger} for associated #it{p}_{T} bin [%2.1f,%2.1f]", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
2205                                          nptbins, ptmin, ptmax,nxeztbins,xeztmin,xeztmax);
2206       fhZTAssocPtBin[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2207       fhZTAssocPtBin[i]->SetYTitle("#it{z}_{T}");
2208
2209
2210       outputContainer->Add(fhXEAssocPtBin[i]);
2211       outputContainer->Add(fhZTAssocPtBin[i]);
2212     }
2213   }
2214
2215   if(fCorrelVzBin)
2216   {
2217     for(Int_t z = 0 ; z < nz ; z++)
2218     {
2219       sz = Form("_vz%d",z);
2220       tz = Form(", #it{v}_{#it{z}} bin %d",z);
2221
2222       fhXEVZ[z]       = new TH2F(Form("hXE%s", sz.Data()),
2223                                  Form("#it{x}_{#it{E}} vs #it{p}_{T trigger}%s", tz.Data()),
2224                                  nptbins, ptmin, ptmax,nxeztbins,xeztmin,xeztmax);
2225       fhXEVZ[z]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2226       fhXEVZ[z]->SetYTitle("#it{x}_{#it{E}}");
2227
2228       fhZTVZ[z]       = new TH2F(Form("hZT%s",sz.Data()),
2229                                  Form("#it{z}_{T} vs #it{p}_{T trigger}%s", tz.Data()),
2230                                  nptbins, ptmin, ptmax,nxeztbins,xeztmin,xeztmax);
2231       fhZTVZ[z]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2232       fhZTVZ[z]->SetYTitle("#it{z}_{T}");
2233
2234       outputContainer->Add(fhXEVZ[z]);
2235       outputContainer->Add(fhZTVZ[z]);
2236     }
2237   }
2238
2239
2240   if(fPi0Trigger)
2241   {
2242     fhPtPi0DecayRatio  = new TH2F
2243     ("hPtPi0DecayRatio","#it{p}_{T} of #pi^{0} and the ratio of pt for two decay",
2244      nptbins,ptmin,ptmax, 100,0.,2.);
2245     fhPtPi0DecayRatio->SetXTitle("#it{p}_{T}^{#pi^{0}} (GeV/#it{c})");
2246     fhPtPi0DecayRatio->SetYTitle("#it{p}_{T}^{Decay}/#it{p}_{T}^{#pi^{0}}");
2247     outputContainer->Add(fhPtPi0DecayRatio) ;
2248
2249     fhDeltaPhiPi0DecayCharged  = new TH2F
2250     ("hDeltaPhiPi0DecayCharged","#phi_{Decay} - #phi_{h^{#pm}} vs #it{p}_{T Decay}",
2251      nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2252     fhDeltaPhiPi0DecayCharged->SetYTitle("#Delta #phi (rad)");
2253     fhDeltaPhiPi0DecayCharged->SetXTitle("#it{p}_{T Decay} (GeV/#it{c})");
2254
2255     fhXEPi0DecayCharged  =
2256     new TH2F("hXEPi0DecayCharged","#it{x}_{#it{E}}  Decay",
2257              nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2258     fhXEPi0DecayCharged->SetYTitle("#it{x}_{#it{E}}");
2259     fhXEPi0DecayCharged->SetXTitle("#it{p}_{T decay} (GeV/#it{c})");
2260
2261     fhZTPi0DecayCharged  =
2262     new TH2F("hZTPi0DecayCharged","#it{z}_{trigger h^{#pm}} = #it{p}_{T h^{#pm}} / #it{p}_{T Decay}",
2263              nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2264     fhZTPi0DecayCharged->SetYTitle("#it{z}_{decay h^{#pm}}");
2265     fhZTPi0DecayCharged->SetXTitle("#it{p}_{T decay} (GeV/#it{c})");
2266
2267     outputContainer->Add(fhDeltaPhiPi0DecayCharged) ;
2268     outputContainer->Add(fhXEPi0DecayCharged) ;
2269     outputContainer->Add(fhZTPi0DecayCharged) ;
2270   }
2271
2272   if(fDecayTrigger)
2273   {
2274     for(Int_t ibit = 0; ibit< fNDecayBits; ibit++)
2275     {
2276       fhDeltaPhiDecayCharged[ibit]  = new TH2F
2277       (Form("hDeltaPhiDecayCharged_bit%d",fDecayBits[ibit]),
2278        Form("#phi_{Decay} - #phi_{h^{#pm}} vs #it{p}_{T Decay}, Bit %d",fDecayBits[ibit]),
2279        nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2280       fhDeltaPhiDecayCharged[ibit]->SetYTitle("#Delta #phi (rad)");
2281       fhDeltaPhiDecayCharged[ibit]->SetXTitle("#it{p}_{T Decay} (GeV/#it{c})");
2282
2283       fhXEDecayCharged[ibit]  =
2284       new TH2F(Form("hXEDecayCharged_bit%d",fDecayBits[ibit]),
2285                Form("#it{x}_{#it{E}}  Decay, Bit %d",fDecayBits[ibit]),
2286                nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2287       fhXEDecayCharged[ibit]->SetYTitle("#it{x}_{#it{E}}");
2288       fhXEDecayCharged[ibit]->SetXTitle("#it{p}_{T decay} (GeV/#it{c})");
2289
2290       fhZTDecayCharged[ibit]  =
2291       new TH2F(Form("hZTDecayCharged_bit%d",fDecayBits[ibit]),
2292                Form("#it{z}_{trigger h^{#pm}} = #it{p}_{T h^{#pm}} / #it{p}_{T Decay}, Bit %d",fDecayBits[ibit]),
2293                nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2294       fhZTDecayCharged[ibit]->SetYTitle("#it{z}_{decay h^{#pm}}");
2295       fhZTDecayCharged[ibit]->SetXTitle("#it{p}_{T decay} (GeV/#it{c})");
2296
2297       outputContainer->Add(fhDeltaPhiDecayCharged[ibit]) ;
2298       outputContainer->Add(fhXEDecayCharged[ibit]) ;
2299       outputContainer->Add(fhZTDecayCharged[ibit]) ;
2300     }
2301   }
2302
2303   //Correlation with neutral hadrons
2304   if(fNeutralCorr)
2305   {
2306     fhDeltaPhiDeltaEtaNeutral  = new TH2F
2307     ("hDeltaPhiDeltaEtaNeutral","#phi_{trigger} - #phi_{h^{0}} vs #eta_{trigger} - #eta_{h^{0}}",
2308      ndeltaphibins ,deltaphimin,deltaphimax, ndeltaetabins ,deltaetamin,deltaetamax);
2309     fhDeltaPhiDeltaEtaNeutral->SetXTitle("#Delta #phi (rad)");
2310     fhDeltaPhiDeltaEtaNeutral->SetYTitle("#Delta #eta");
2311
2312     fhPhiNeutral  = new TH2F
2313     ("hPhiNeutral","#phi_{#pi^{0}}  vs #it{p}_{T #pi^{0}}",
2314      nptbins,ptmin,ptmax,180,0,TMath::TwoPi());
2315     fhPhiNeutral->SetYTitle("#phi_{#pi^{0}} (rad)");
2316     fhPhiNeutral->SetXTitle("#it{p}_{T #pi^{0}} (GeV/#it{c})");
2317
2318     fhEtaNeutral  = new TH2F
2319     ("hEtaNeutral","#eta_{#pi^{0}}  vs #it{p}_{T #pi^{0}}",
2320      nptbins,ptmin,ptmax,200,-1.,1.);
2321     fhEtaNeutral->SetYTitle("#eta_{#pi^{0}} (rad)");
2322     fhEtaNeutral->SetXTitle("#it{p}_{T #pi^{0}} (GeV/#it{c})");
2323
2324     fhDeltaPhiNeutral  = new TH2F
2325     ("hDeltaPhiNeutral","#phi_{trigger} - #phi_{#pi^{0}} vs #it{p}_{T trigger}",
2326      nptbins,ptmin,ptmax,ndeltaphibins ,deltaphimin,deltaphimax);
2327     fhDeltaPhiNeutral->SetYTitle("#Delta #phi (rad)");
2328     fhDeltaPhiNeutral->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2329
2330     fhDeltaPhiNeutralPt  = new TH2F
2331     ("hDeltaPhiNeutralPt","#phi_{trigger} - #phi_{#pi^{0}} vs #it{p}_{T #pi^{0}}}",
2332      nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2333     fhDeltaPhiNeutralPt->SetYTitle("#Delta #phi (rad)");
2334     fhDeltaPhiNeutralPt->SetXTitle("#it{p}_{T h^{0}} (GeV/#it{c})");
2335
2336     fhDeltaEtaNeutral  = new TH2F
2337     ("hDeltaEtaNeutral","#eta_{trigger} - #eta_{#pi^{0}} vs #it{p}_{T trigger}",
2338      nptbins,ptmin,ptmax, ndeltaetabins ,deltaetamin,deltaetamax);
2339     fhDeltaEtaNeutral->SetYTitle("#Delta #eta");
2340     fhDeltaEtaNeutral->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2341
2342     fhXENeutral  =
2343     new TH2F("hXENeutral","#it{x}_{#it{E}} for #pi^{0} associated",
2344              nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2345     fhXENeutral->SetYTitle("#it{x}_{#it{E}}");
2346     fhXENeutral->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2347
2348     fhPtHbpXENeutral  =
2349     new TH2F("hHbpXENeutral","#xi = ln(1/#it{x}_{#it{E}})for #pi^{0} associated",
2350              nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
2351     fhPtHbpXENeutral->SetYTitle("ln(1/#it{x}_{#it{E}})");
2352     fhPtHbpXENeutral->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2353
2354     fhZTNeutral  =
2355     new TH2F("hZTNeutral","#it{z}_{trigger #pi} = #it{p}_{T #pi^{0}} / #it{p}_{T trigger} for #pi^{0} associated",
2356              nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2357     fhZTNeutral->SetYTitle("#it{z}_{trigger #pi^{0}}");
2358     fhZTNeutral->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2359
2360     fhPtHbpZTNeutral  =
2361     new TH2F("hHbpZTNeutral","#xi = ln(1/#it{x}_{#it{E}}) for #pi^{0} associated",
2362              nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
2363     fhPtHbpZTNeutral->SetYTitle("ln(1/#it{z}_{T})");
2364     fhPtHbpZTNeutral->SetXTitle("#it{p}_{T trigger}");
2365
2366     outputContainer->Add(fhDeltaPhiDeltaEtaNeutral);
2367     outputContainer->Add(fhPhiNeutral) ;
2368     outputContainer->Add(fhEtaNeutral) ;
2369     outputContainer->Add(fhDeltaPhiNeutral) ;
2370     outputContainer->Add(fhDeltaPhiNeutralPt) ;
2371     outputContainer->Add(fhDeltaEtaNeutral) ;
2372     outputContainer->Add(fhXENeutral) ;
2373     outputContainer->Add(fhPtHbpXENeutral) ;
2374     outputContainer->Add(fhZTNeutral) ;
2375     outputContainer->Add(fhPtHbpZTNeutral) ;
2376
2377     fhDeltaPhiUeNeutralPt  = new TH2F
2378     (Form("hDeltaPhiUe%sNeutralPt",right.Data()),"#phi_{trigger} - #phi_{#pi^{0}} vs #it{p}_{T #pi^{0}}}",
2379      nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2380     fhDeltaPhiUeNeutralPt->SetYTitle("#Delta #phi (rad)");
2381     fhDeltaPhiUeNeutralPt->SetXTitle("#it{p}_{T h^{0}} (GeV/#it{c})");
2382
2383     fhXEUeNeutral  =
2384     new TH2F(Form("hXEUeNeutral%s",right.Data()),"#it{x}_{#it{E}} for #pi^{0} associated",
2385              nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2386     fhXEUeNeutral->SetYTitle("#it{x}_{#it{E}}");
2387     fhXEUeNeutral->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2388
2389     fhPtHbpXEUeNeutral  =
2390     new TH2F(Form("hHbpXEUeNeutral%s",right.Data()),"#xi = ln(1/#it{x}_{#it{E}}) for #pi^{0} associated",
2391              nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
2392     fhPtHbpXEUeNeutral->SetYTitle("ln(1/#it{x}_{#it{E}})");
2393     fhPtHbpXEUeNeutral->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2394
2395     fhZTUeNeutral  =
2396     new TH2F(Form("hZTUeNeutral%s",right.Data()),"#it{z}_{trigger #pi} = #it{p}_{T #pi^{0}} / #it{p}_{T trigger} for #pi^{0} associated",
2397              nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2398     fhZTUeNeutral->SetYTitle("#it{z}_{trigger #pi^{0}}");
2399     fhZTUeNeutral->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2400
2401     fhPtHbpZTUeNeutral  =
2402     new TH2F(Form("hHbpZTUeNeutral%s",right.Data()),"#xi = ln(1/#it{x}_{#it{E}}) for #pi^{0} associated",
2403              nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
2404     fhPtHbpXEUeNeutral->SetYTitle("ln(1/#it{z}_{T})");
2405     fhPtHbpXEUeNeutral->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2406
2407     outputContainer->Add(fhDeltaPhiUeNeutralPt) ;
2408     outputContainer->Add(fhXEUeNeutral) ;
2409     outputContainer->Add(fhPtHbpXEUeNeutral) ;
2410     outputContainer->Add(fhZTUeNeutral) ;
2411     outputContainer->Add(fhPtHbpZTUeNeutral) ;
2412
2413     if(fMakeSeveralUE)
2414     {
2415       fhDeltaPhiUeLeftNeutral  = new TH2F
2416       ("hDeltaPhiUeLeftNeutralPt","#phi_{trigger} - #phi_{#Ueh^{0}} vs #it{p}_{T h^{0}} with neutral UE left side range of trigger particles",
2417        nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2418       fhDeltaPhiUeLeftNeutral->SetYTitle("#Delta #phi (rad)");
2419       fhDeltaPhiUeLeftNeutral->SetXTitle("#it{p}_{T h^{0}} (GeV/#it{c})");
2420       outputContainer->Add(fhDeltaPhiUeLeftNeutral) ;
2421
2422       fhXEUeLeftNeutral  =
2423       new TH2F("hXEUeNeutralLeft","#it{x}_{#it{E}} = #it{p}_{T Ueh^{0}} / #it{p}_{T trigger} with neutral UE left side of trigger",
2424                nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2425       fhXEUeLeftNeutral->SetYTitle("#it{z}_{trigger Ueh^{0}}");
2426       fhXEUeLeftNeutral->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2427       outputContainer->Add(fhXEUeLeftNeutral) ;
2428
2429       fhPtHbpXEUeLeftNeutral  =
2430       new TH2F("hHbpXEUeNeutralLeft","#xi = ln(1/#it{x}_{#it{E}}) with neutral UE left side of trigger",
2431                nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
2432       fhPtHbpXEUeLeftNeutral->SetYTitle("ln(1/#it{x}_{#it{E}})");
2433       fhPtHbpXEUeLeftNeutral->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2434       outputContainer->Add(fhPtHbpXEUeLeftNeutral) ;
2435
2436       fhZTUeLeftNeutral  =
2437       new TH2F("hZTUeNeutralLeft","#it{z}_{trigger h^{0}} = #it{p}_{T Ueh^{0}} / #it{p}_{T trigger} with neutral UE left side of trigger",
2438                nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2439       fhZTUeLeftNeutral->SetYTitle("#it{z}_{trigger Ueh^{0}}");
2440       fhZTUeLeftNeutral->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2441       outputContainer->Add(fhZTUeLeftNeutral) ;
2442
2443       fhPtHbpZTUeLeftNeutral  =
2444       new TH2F("hHbpZTUeNeutralLeft","#xi = ln(1/#it{z}_{T}) with neutral UE left side of trigger",
2445                nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
2446       fhPtHbpZTUeLeftNeutral->SetYTitle("ln(1/#it{z}_{T})");
2447       fhPtHbpZTUeLeftNeutral->SetXTitle("#it{p}_{T trigger}");
2448       outputContainer->Add(fhPtHbpZTUeLeftNeutral) ;
2449     }
2450
2451     if(fPi0Trigger)
2452     {
2453       fhDeltaPhiPi0DecayNeutral  = new TH2F
2454       ("hDeltaPhiPi0DecayNeutral","#phi_{Decay} - #phi_{h^{0}} vs #it{p}_{T Decay}",
2455        nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2456       fhDeltaPhiPi0DecayNeutral->SetYTitle("#Delta #phi (rad)");
2457       fhDeltaPhiPi0DecayNeutral->SetXTitle("#it{p}_{T Decay} (GeV/#it{c})");
2458
2459       fhXEPi0DecayNeutral  =
2460       new TH2F("hXEPi0DecayNeutral","#it{x}_{#it{E}} for decay trigger",
2461                nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2462       fhXEPi0DecayNeutral->SetYTitle("#it{x}_{#it{E}}");
2463       fhXEPi0DecayNeutral->SetXTitle("#it{p}_{T decay}");
2464
2465       fhZTPi0DecayNeutral  =
2466       new TH2F("hZTPi0DecayNeutral","#it{z}_{trigger h^{0}} = #it{p}_{T h^{0}} / #it{p}_{T Decay}",
2467                nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2468       fhZTPi0DecayNeutral->SetYTitle("#it{z}_{h^{0}}");
2469       fhZTPi0DecayNeutral->SetXTitle("#it{p}_{T decay}");
2470
2471       outputContainer->Add(fhDeltaPhiPi0DecayNeutral) ;
2472       outputContainer->Add(fhXEPi0DecayNeutral) ;
2473       outputContainer->Add(fhZTPi0DecayNeutral) ;
2474     }
2475   }//Correlation with neutral hadrons
2476
2477   // If data is MC, fill more histograms, depending on origin
2478   if(IsDataMC())
2479   {
2480     for(Int_t i= fMCGenTypeMin; i <= fMCGenTypeMax; i++)
2481     {
2482       fhMCPtTrigger[i]  = new TH1F (Form("hMCPtTrigger_%s",nameMC[i].Data()),
2483                                  Form("MC %s: generated trigger #it{p}_{T}",nameMC[i].Data()),
2484                                  nptbins,ptmin,ptmax);
2485       fhMCPtTrigger[i]->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
2486
2487       fhMCPhiTrigger[i]  = new TH2F (Form("hMCPhiTrigger_%s",nameMC[i].Data()),
2488                                      Form("MC %s: generated trigger #phi",nameMC[i].Data()),
2489                                      nptbins,ptmin,ptmax, nphibins,phimin,phimax);
2490       fhMCPhiTrigger[i]->SetYTitle("#phi (rad)");
2491       fhMCPhiTrigger[i]->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
2492
2493       fhMCEtaTrigger[i]  = new TH2F (Form("hMCEtaTrigger_%s",nameMC[i].Data()),
2494                                      Form("MC %s: generated trigger #eta",nameMC[i].Data()),
2495                                      nptbins,ptmin,ptmax, netabins,etamin,etamax);
2496       fhMCEtaTrigger[i]->SetYTitle("#eta");
2497       fhMCEtaTrigger[i]->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
2498
2499       if(fMakeAbsoluteLeading || fMakeNearSideLeading)
2500       {
2501         fhMCPtTriggerNotLeading[i]  = new TH1F (Form("hMCPtTriggerNotLeading_%s",nameMC[i].Data()),
2502                                                 Form("MC %s: generated trigger #it{p}_{T}, when not leading of primaries",nameMC[i].Data()),
2503                                                 nptbins,ptmin,ptmax);
2504         fhMCPtTriggerNotLeading[i]->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
2505
2506         fhMCPhiTriggerNotLeading[i]  = new TH2F (Form("hMCPhiTriggerNotLeading_%s",nameMC[i].Data()),
2507                                                  Form("MC %s: generated trigger #phi, when not leading of primaries",nameMC[i].Data()),
2508                                                  nptbins,ptmin,ptmax, nphibins,phimin,phimax);
2509         fhMCPhiTriggerNotLeading[i]->SetYTitle("#phi (rad)");
2510         fhMCPhiTriggerNotLeading[i]->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
2511
2512
2513         fhMCEtaTriggerNotLeading[i]  = new TH2F (Form("hMCEtaTriggerNotLeading_%s",nameMC[i].Data()),
2514                                                  Form("MC %s: generated triogger #eta, when not leading of primaries",nameMC[i].Data()),
2515                                                  nptbins,ptmin,ptmax, netabins,etamin,etamax);
2516         fhMCEtaTriggerNotLeading[i]->SetYTitle("#eta ");
2517         fhMCEtaTriggerNotLeading[i]->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
2518       }
2519
2520       fhMCEtaCharged[i]  = new TH2F (Form("hMCEtaCharged_%s",nameMC[i].Data()),
2521                                      Form("MC %s: #eta_{h^{#pm}}  vs #it{p}_{T #pm}",nameMC[i].Data()),
2522                                      nptbins,ptmin,ptmax,100,-1.,1.);
2523       fhMCEtaCharged[i]->SetYTitle("#eta_{h^{#pm}} (rad)");
2524       fhMCEtaCharged[i]->SetXTitle("#it{p}_{T #pm} (GeV/#it{c})");
2525
2526       fhMCPhiCharged[i]  = new TH2F(Form("hMCPhiCharged_%s",nameMC[i].Data()),
2527                                     Form("MC %s: phi_{h^{#pm}}  vs #it{p}_{T #pm}",nameMC[i].Data()),
2528                                     nptbins,ptmin,ptmax,180,0,TMath::TwoPi());
2529       fhMCPhiCharged[i]->SetYTitle("MC #phi_{h^{#pm}} (rad)");
2530       fhMCPhiCharged[i]->SetXTitle("#it{p}_{T #pm} (GeV/#it{c})");
2531
2532       fhMCDeltaPhiDeltaEtaCharged[i]  = new TH2F (Form("hMCDeltaPhiDeltaEtaCharged_%s",nameMC[i].Data()),
2533                                                   Form("MC %s: phi_{trigger} - #phi_{h^{#pm}} vs #eta_{trigger} - #eta_{h^{#pm}}",nameMC[i].Data()),
2534                                                   ndeltaphibins ,deltaphimin,deltaphimax,ndeltaetabins ,deltaetamin,deltaetamax);
2535       fhMCDeltaPhiDeltaEtaCharged[i]->SetXTitle("#Delta #phi (rad)");
2536       fhMCDeltaPhiDeltaEtaCharged[i]->SetYTitle("#Delta #eta");
2537
2538       fhMCDeltaEtaCharged[i]  = new TH2F (Form("hMCDeltaEtaCharged_%s",nameMC[i].Data()),
2539                                           Form("MC %s: #eta_{trigger} - #eta_{h^{#pm}} vs #it{p}_{T trigger} and #it{p}_{T assoc}",nameMC[i].Data()),
2540                                           nptbins,ptmin,ptmax,ndeltaetabins ,deltaetamin,deltaetamax);
2541       fhMCDeltaEtaCharged[i]->SetYTitle("#Delta #eta");
2542       fhMCDeltaEtaCharged[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2543
2544       fhMCDeltaPhiCharged[i]  = new TH2F  (Form("hMCDeltaPhiCharged_%s",nameMC[i].Data()),
2545                                            Form("MC %s: #phi_{trigger} - #phi_{h^{#pm}} vs #it{p}_{T trigger}",nameMC[i].Data()),
2546                                            nptbins,ptmin,ptmax,ndeltaphibins ,deltaphimin,deltaphimax);
2547       fhMCDeltaPhiCharged[i]->SetYTitle("#Delta #phi (rad)");
2548       fhMCDeltaPhiCharged[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2549
2550       fhMCDeltaPhiChargedPt[i]  = new TH2F (Form("hMCDeltaPhiChargedPt_%s",nameMC[i].Data()),
2551                                             Form("MC %s: #phi_{trigger} - #phi_{#h^{#pm}} vs #it{p}_{T h^{#pm}}",nameMC[i].Data()),
2552                                             nptbins,ptmin,ptmax,ndeltaphibins ,deltaphimin,deltaphimax);
2553       fhMCDeltaPhiChargedPt[i]->SetYTitle("#Delta #phi (rad)");
2554       fhMCDeltaPhiChargedPt[i]->SetXTitle("#it{p}_{T h^{#pm}} (GeV/#it{c})");
2555
2556       fhMCPtXECharged[i]  = new TH2F (Form("hMCPtXECharged_%s",nameMC[i].Data()),
2557                                       Form("MC %s: #it{x}_{#it{E}} with charged hadrons",nameMC[i].Data()),
2558                                       nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2559       fhMCPtXECharged[i]->SetYTitle("#it{x}_{#it{E}}");
2560       fhMCPtXECharged[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2561
2562       fhMCPtHbpXECharged[i]  = new TH2F(Form("hMCHbpXECharged_%s",nameMC[i].Data()),
2563                                         Form("MC %s: #xi = ln(1/#it{x}_{#it{E}}) with charged hadrons",nameMC[i].Data()),
2564                                         nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
2565       fhMCPtHbpXECharged[i]->SetYTitle("ln(1/#it{x}_{#it{E}})");
2566       fhMCPtHbpXECharged[i]->SetXTitle("#it{p}_{T trigger}");
2567
2568       fhMCPtZTCharged[i]  = new TH2F(Form("hMCPtZTCharged_%s",nameMC[i].Data()),
2569                                      Form("MC %s: #it{z}_{T} with charged hadrons",nameMC[i].Data()),
2570                                      nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2571       fhMCPtZTCharged[i]->SetYTitle("#it{z}_{T}");
2572       fhMCPtZTCharged[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2573
2574       fhMCPtHbpZTCharged[i]  = new TH2F(Form("hMCHbpZTCharged_%s",nameMC[i].Data()),
2575                                         Form("MC %s: #xi = ln(1/#it{z}_{T}) with charged hadrons",nameMC[i].Data()),
2576                                              nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
2577       fhMCPtHbpZTCharged[i]->SetYTitle("ln(1/#it{z}_{T})");
2578       fhMCPtHbpZTCharged[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2579
2580       fhMCPtTrigPout[i]  = new TH2F(Form("hMCPtTrigPout_%s",nameMC[i].Data()),
2581                                     Form("MC %s: #it{p}_{out} with triggers",nameMC[i].Data()),
2582                                     nptbins,ptmin,ptmax,nptbins,-1.*ptmax/2.,ptmax/2.);
2583       fhMCPtTrigPout[i]->SetYTitle("#it{p}_{out} (GeV/#it{c})");
2584       fhMCPtTrigPout[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2585
2586       fhMCPtAssocDeltaPhi[i]  = new TH2F(Form("hMCPtAssocDeltaPhi_%s",nameMC[i].Data()),
2587                                          Form("MC %s: #Delta #phi with associated charged hadrons",nameMC[i].Data()),
2588                                          nptbins,ptmin,ptmax,ndeltaphibins ,deltaphimin,deltaphimax);
2589       fhMCPtAssocDeltaPhi[i]->SetYTitle("#Delta #phi (rad)");
2590       fhMCPtAssocDeltaPhi[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2591
2592       outputContainer->Add(fhMCPtTrigger[i]);
2593       outputContainer->Add(fhMCPhiTrigger[i]);
2594       outputContainer->Add(fhMCEtaTrigger[i]);
2595
2596       if(fMakeAbsoluteLeading || fMakeNearSideLeading)
2597       {
2598         outputContainer->Add(fhMCPtTriggerNotLeading[i]);
2599         outputContainer->Add(fhMCPhiTriggerNotLeading[i]);
2600         outputContainer->Add(fhMCEtaTriggerNotLeading[i]);
2601       }
2602
2603       outputContainer->Add(fhMCDeltaPhiDeltaEtaCharged[i]);
2604       outputContainer->Add(fhMCPhiCharged[i]) ;
2605       outputContainer->Add(fhMCEtaCharged[i]) ;
2606       outputContainer->Add(fhMCDeltaEtaCharged[i]) ;
2607       outputContainer->Add(fhMCDeltaPhiCharged[i]) ;
2608
2609       outputContainer->Add(fhMCDeltaPhiChargedPt[i]) ;
2610       outputContainer->Add(fhMCPtXECharged[i]) ;
2611       outputContainer->Add(fhMCPtZTCharged[i]) ;
2612       outputContainer->Add(fhMCPtHbpXECharged[i]) ;
2613       outputContainer->Add(fhMCPtHbpZTCharged[i]) ;
2614       outputContainer->Add(fhMCPtTrigPout[i]) ;
2615       outputContainer->Add(fhMCPtAssocDeltaPhi[i]) ;
2616
2617       // Underlying event
2618
2619       fhMCUePart[i]  =
2620       new TH1F(Form("hMCUePart_%s",nameMC[i].Data()),
2621                Form("MC %s: UE particles distribution vs #it{p}_{T trigger}",nameMC[i].Data()),
2622                nptbins,ptmin,ptmax);
2623       fhMCUePart[i]->SetYTitle("#it{dN}^{ch}");
2624       fhMCUePart[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2625
2626       fhMCPtXEUeCharged[i]  =
2627       new TH2F(Form("hMCPtXEUeCharged%s_%s",right.Data(),nameMC[i].Data()),
2628                Form("MC %s: #it{x}_{#it{E}} with charged hadrons, Underlying Event %s",nameMC[i].Data(),right.Data()),
2629                nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2630       fhMCPtXEUeCharged[i]->SetYTitle("#it{x}_{#it{E}}");
2631       fhMCPtXEUeCharged[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2632
2633       fhMCPtHbpXEUeCharged[i] =
2634       new TH2F(Form("hMCPtHbpXEUeCharged%s_%s",right.Data(),nameMC[i].Data()),
2635                Form("MC %s: #xi = ln(1/#it{x}_{#it{E}}) with charged hadrons, Underlying Event %s",nameMC[i].Data(),right.Data()),
2636                nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
2637       fhMCPtHbpXEUeCharged[i]->SetYTitle("ln(1/#it{x}_{#it{E}})");
2638       fhMCPtHbpXEUeCharged[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2639
2640       fhMCPtZTUeCharged[i] =
2641       new TH2F(Form("hMCPtZTUeCharged%s_%s",right.Data(),nameMC[i].Data()),
2642                Form("MC %s: #it{z}_{T} with charged hadrons, Underlying Event %s",nameMC[i].Data(),right.Data()),
2643                nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2644       fhMCPtZTUeCharged[i]->SetYTitle("#it{z}_{T}");
2645       fhMCPtZTUeCharged[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2646
2647       fhMCPtHbpZTUeCharged[i] =
2648       new TH2F(Form("hMCPtHbpZTUeCharged%s_%s",right.Data(),nameMC[i].Data()),
2649                Form("MC %s: #xi = ln(1/#it{z}_{T}) with charged hadrons, Underlying Event %s",nameMC[i].Data(),right.Data()),
2650                nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
2651       fhMCPtHbpZTUeCharged[i]->SetYTitle("ln(1/#it{z}_{T})");
2652       fhMCPtHbpZTUeCharged[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2653
2654       outputContainer->Add(fhMCUePart[i]);
2655       outputContainer->Add(fhMCPtXEUeCharged[i]) ;
2656       outputContainer->Add(fhMCPtZTUeCharged[i]) ;
2657       outputContainer->Add(fhMCPtHbpZTUeCharged[i]);
2658       outputContainer->Add(fhMCPtHbpXEUeCharged[i]);
2659
2660       if(fMakeSeveralUE)
2661       {
2662         fhMCPtXEUeLeftCharged[i]  = new TH2F(Form("hMCPtXEUeChargedLeft_%s",nameMC[i].Data()),
2663                                              Form("MC %s: #it{x}_{#it{E}} with charged hadrons, with UE left side range of trigger particles",nameMC[i].Data()),
2664                                              nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2665         fhMCPtXEUeLeftCharged[i]->SetYTitle("#it{x}_{#it{E}}");
2666         fhMCPtXEUeLeftCharged[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2667
2668         fhMCPtHbpXEUeLeftCharged[i] = new TH2F(Form("hMCPtHbpXEUeChargedLeft_%s",nameMC[i].Data()),
2669                                                Form("MC %s: #xi = ln(1/#it{x}_{#it{E}}) with charged hadrons, with UE left side range of trigger particles",nameMC[i].Data()),
2670                                                nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
2671         fhMCPtHbpXEUeLeftCharged[i]->SetYTitle("ln(1/#it{x}_{#it{E}})");
2672         fhMCPtHbpXEUeLeftCharged[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2673
2674         fhMCPtZTUeLeftCharged[i]  = new TH2F(Form("hMCPtZTUeChargedLeft_%s",nameMC[i].Data()),
2675                                              Form("MC %s: #it{z}_{T} with charged hadrons, with UE left side range of trigger particles",nameMC[i].Data()),
2676                                              nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2677         fhMCPtZTUeLeftCharged[i]->SetYTitle("#it{z}_{T}");
2678         fhMCPtZTUeLeftCharged[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2679
2680         fhMCPtHbpZTUeLeftCharged[i] = new TH2F(Form("hMCPtHbpZTUeChargedLeft_%s",nameMC[i].Data()),
2681                                                Form("MC %s: #xi = ln(1/#it{z}_{T}) with charged hadrons, with UE left side range of trigger particles",nameMC[i].Data()),
2682                                                nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
2683         fhMCPtHbpZTUeLeftCharged[i]->SetYTitle("ln(1/#it{z}_{T})");
2684         fhMCPtHbpZTUeLeftCharged[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2685
2686         outputContainer->Add(fhMCPtXEUeLeftCharged[i]) ;
2687         outputContainer->Add(fhMCPtZTUeLeftCharged[i]) ;
2688         outputContainer->Add(fhMCPtHbpXEUeLeftCharged[i]);
2689         outputContainer->Add(fhMCPtHbpZTUeLeftCharged[i]) ;
2690
2691       }
2692     }
2693   } //for MC histogram
2694
2695   if(DoOwnMix())
2696   {
2697     //create event containers
2698
2699     if(!fUseMixStoredInReader || (fUseMixStoredInReader && !GetReader()->ListWithMixedEventsForTracksExists()))
2700     {
2701       Int_t nvz = GetNZvertBin();
2702       Int_t nrp = GetNRPBin();
2703       Int_t nce = GetNCentrBin();
2704
2705       fListMixTrackEvents= new TList*[nvz*nrp*nce] ;
2706
2707       for( Int_t ice = 0 ; ice < nce ; ice++ )
2708       {
2709         for( Int_t ivz = 0 ; ivz < nvz ; ivz++ )
2710         {
2711           for( Int_t irp = 0 ; irp < nrp ; irp++ )
2712           {
2713             Int_t bin = GetEventMixBin(ice,ivz,irp); //ic*nvz*nrp+iz*nrp+irp;
2714
2715             //printf("GetCreateOutputObjects - Bins : cent %d, vz %d, RP %d, event %d/%d\n",
2716             //       ic,iz, irp, bin);
2717
2718             fListMixTrackEvents[bin] = new TList() ;
2719             fListMixTrackEvents[bin]->SetOwner(kFALSE);
2720           }
2721         }
2722       }
2723     }
2724
2725     fhPtTriggerMixed  = new TH1F ("hPtTriggerMixed","#it{p}_{T} distribution of trigger particles, used for mixing", nptbins,ptmin,ptmax);
2726     fhPtTriggerMixed->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
2727
2728     if(fCorrelVzBin)
2729     {
2730       fhPtTriggerMixedVzBin  = new TH2F ("hPtTriggerMixedVzBin","#it{p}_{T} distribution of trigger particles, used for mixing", nptbins,ptmin,ptmax,GetNZvertBin(),0,GetNZvertBin());
2731       fhPtTriggerMixedVzBin->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
2732       fhPtTriggerMixedVzBin->SetYTitle("#it{v}_{#it{z}} bin");
2733       outputContainer->Add(fhPtTriggerMixedVzBin);
2734     }
2735
2736     fhPtTriggerMixedBin  = new TH2F ("hPtTriggerMixedBin","#it{p}_{T} distribution of trigger particles vs mixing bin", nptbins,ptmin,ptmax,nMixBins,0,nMixBins);
2737     fhPtTriggerMixedBin->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
2738     fhPtTriggerMixedBin->SetYTitle("Bin");
2739
2740     fhPhiTriggerMixed  = new TH2F ("hPhiTriggerMixed","#phi distribution of trigger Particles, used for mixing",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
2741     fhPhiTriggerMixed->SetYTitle("#phi (rad)");
2742
2743     fhEtaTriggerMixed  = new TH2F ("hEtaTriggerMixed","#eta distribution of trigger, used for mixing",nptbins,ptmin,ptmax, netabins,etamin,etamax);
2744     fhEtaTriggerMixed->SetYTitle("#eta ");
2745
2746     outputContainer->Add(fhPtTriggerMixed);
2747     outputContainer->Add(fhPtTriggerMixedBin);
2748     outputContainer->Add(fhPhiTriggerMixed);
2749     outputContainer->Add(fhEtaTriggerMixed);
2750
2751     // Fill the cluster pool only in isolation analysis or if requested
2752     if( neutralMix && (!fUseMixStoredInReader || (fUseMixStoredInReader && !GetReader()->ListWithMixedEventsForCaloExists())))
2753     {
2754       Int_t nvz = GetNZvertBin();
2755       Int_t nrp = GetNRPBin();
2756       Int_t nce = GetNCentrBin();
2757
2758       fListMixCaloEvents= new TList*[nvz*nrp*nce] ;
2759
2760       for( Int_t ice = 0 ; ice < nce ; ice++ )
2761       {
2762         for( Int_t ivz = 0 ; ivz < nvz ; ivz++ )
2763         {
2764           for( Int_t irp = 0 ; irp < nrp ; irp++ )
2765           {
2766             Int_t bin = GetEventMixBin(ice,ivz,irp); //ic*nvz*nrp+iz*nrp+irp;
2767
2768             //printf("GetCreateOutputObjects - Bins : cent %d, vz %d, RP %d, event %d/%d\n",
2769             //       ic,iz, irp, bin);
2770
2771             fListMixCaloEvents[bin] = new TList() ;
2772             fListMixCaloEvents[bin]->SetOwner(kFALSE);
2773           }
2774         }
2775       }
2776     }
2777
2778     //Init the list in the reader if not done previously
2779     if(fUseMixStoredInReader)
2780     {
2781       if( !GetReader()->ListWithMixedEventsForTracksExists() )
2782         GetReader()->SetListWithMixedEventsForTracks(fListMixTrackEvents);
2783
2784       if( !GetReader()->ListWithMixedEventsForCaloExists()   )
2785         GetReader()->SetListWithMixedEventsForCalo  (fListMixCaloEvents );
2786     }
2787
2788     fhEventBin=new TH1I("hEventBin","Number of triggers per bin(cen,vz,rp)",
2789                         GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1,0,
2790                         GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1) ;
2791     fhEventBin->SetXTitle("event bin");
2792     outputContainer->Add(fhEventBin) ;
2793
2794     fhEventMixBin=new TH1I("hEventMixBin","Number of triggers mixed per event bin(cen,vz,rp)",
2795                            GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1,0,
2796                            GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1) ;
2797     fhEventMixBin->SetXTitle("event bin");
2798     outputContainer->Add(fhEventMixBin) ;
2799
2800     fhEventMBBin=new TH1I("hEventMBBin","Number of min bias events per bin(cen,vz,rp)",
2801                            GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1,0,
2802                            GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1) ;
2803     fhEventMBBin->SetXTitle("event bin");
2804     outputContainer->Add(fhEventMBBin) ;
2805
2806     fhNtracksMB=new TH2F("hNtracksMBEvent","Number of filtered tracks in MB event per event bin",ntrbins,trmin,trmax,
2807                          GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1,0,
2808                          GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1) ;
2809     fhNtracksMB->SetYTitle("event bin");
2810     fhNtracksMB->SetXTitle("#it{N}_{track}");
2811     outputContainer->Add(fhNtracksMB);
2812
2813     if( neutralMix )
2814     {
2815       fhNclustersMB=new TH2F("hNclustersMBEvent","Number of filtered clusters in MB events per event bin",nclbins,clmin,clmax,
2816                              GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1,0,
2817                              GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1) ;
2818       fhNclustersMB->SetYTitle("event bin");
2819       fhNclustersMB->SetXTitle("#it{N}_{cluster}");
2820       outputContainer->Add(fhNclustersMB);
2821     }
2822
2823     fhMixDeltaPhiCharged  = new TH2F
2824     ("hMixDeltaPhiCharged","Mixed event : #phi_{trigger} - #phi_{h^{#pm}} vs #it{p}_{T trigger}",
2825      nptbins,ptmin,ptmax,ndeltaphibins ,deltaphimin,deltaphimax);
2826     fhMixDeltaPhiCharged->SetYTitle("#Delta #phi (rad)");
2827     fhMixDeltaPhiCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2828     outputContainer->Add(fhMixDeltaPhiCharged);
2829
2830     fhMixDeltaPhiDeltaEtaCharged  = new TH2F
2831     ("hMixDeltaPhiDeltaEtaCharged","Mixed event : #phi_{trigger} - #phi_{h^{#pm}} vs #eta_{trigger} - #eta_{h^{#pm}}",
2832      ndeltaphibins ,deltaphimin,deltaphimax,ndeltaetabins ,deltaetamin,deltaetamax);
2833     fhMixDeltaPhiDeltaEtaCharged->SetXTitle("#Delta #phi (rad)");
2834     fhMixDeltaPhiDeltaEtaCharged->SetYTitle("#Delta #eta");
2835     outputContainer->Add(fhMixDeltaPhiDeltaEtaCharged);
2836
2837     fhMixXECharged  =
2838     new TH2F("hMixXECharged","Mixed event : #it{x}_{#it{E}} for charged tracks",
2839              nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2840     fhMixXECharged->SetYTitle("#it{x}_{#it{E}}");
2841     fhMixXECharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2842     outputContainer->Add(fhMixXECharged);
2843
2844     fhMixXEUeCharged  =
2845     new TH2F("hMixXEUeCharged","Mixed event : #it{x}_{#it{E}} for charged tracks in Ue region",
2846              nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
2847     fhMixXEUeCharged->SetYTitle("#it{x}_{#it{E}}");
2848     fhMixXEUeCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2849     outputContainer->Add(fhMixXEUeCharged);
2850
2851     fhMixHbpXECharged  =
2852     new TH2F("hMixHbpXECharged","mixed event : #xi = ln(1/#it{x}_{#it{E}}) with charged hadrons",
2853              nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
2854     fhMixHbpXECharged->SetYTitle("ln(1/#it{x}_{#it{E}})");
2855     fhMixHbpXECharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2856     outputContainer->Add(fhMixHbpXECharged);
2857
2858     fhMixDeltaPhiChargedAssocPtBin         = new TH2F*[fNAssocPtBins*nz];
2859     fhMixDeltaPhiChargedAssocPtBinDEta08   = new TH2F*[fNAssocPtBins*nz];
2860     fhMixDeltaPhiChargedAssocPtBinDEta0    = new TH2F*[fNAssocPtBins*nz];
2861     fhMixDeltaPhiDeltaEtaChargedAssocPtBin = new TH2F*[fNAssocPtBins*nz];
2862
2863     for(Int_t i = 0 ; i < fNAssocPtBins ; i++)
2864     {
2865       for(Int_t z = 0 ; z < nz ; z++)
2866       {
2867         Int_t bin = i*nz+z;
2868
2869         if(fCorrelVzBin)
2870         {
2871           sz = Form("_vz%d",z);
2872           tz = Form(", #it{v}_{#it{z}} bin %d",z);
2873         }
2874
2875         //printf("MIX : iAssoc %d, Vz %d, bin %d - sz %s, tz %s \n",i,z,bin,sz.Data(),tz.Data());
2876
2877         fhMixDeltaPhiDeltaEtaChargedAssocPtBin[bin] = new TH2F(Form("hMixDeltaPhiDeltaEtaChargedAssocPtBin%2.1f_%2.1f%s", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],sz.Data()),
2878                                                                Form("Mixed event #Delta #phi vs #it{p}_{T trigger} for associated #it{p}_{T} bin [%2.1f,%2.1f]%s", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],tz.Data()),
2879                                                                ndeltaphibins ,deltaphimin,deltaphimax,ndeltaetabins ,deltaetamin,deltaetamax);
2880         fhMixDeltaPhiDeltaEtaChargedAssocPtBin[bin]->SetXTitle("#Delta #phi (rad)");
2881         fhMixDeltaPhiDeltaEtaChargedAssocPtBin[bin]->SetYTitle("#Delta #eta");
2882
2883         outputContainer->Add(fhMixDeltaPhiDeltaEtaChargedAssocPtBin[bin]);
2884
2885         fhMixDeltaPhiChargedAssocPtBin[bin] = new TH2F(Form("hMixDeltaPhiChargedAssocPtBin%2.1f_%2.1f%s", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],sz.Data()),
2886                                                        Form("Mixed event #Delta #phi vs #it{p}_{T trigger} for associated #it{p}_{T} bin [%2.1f,%2.1f]%s", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],tz.Data()),
2887                                                        nptbins, ptmin, ptmax,  ndeltaphibins ,deltaphimin,deltaphimax);
2888         fhMixDeltaPhiChargedAssocPtBin[bin]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2889         fhMixDeltaPhiChargedAssocPtBin[bin]->SetYTitle("#Delta #phi (rad)");
2890
2891         outputContainer->Add(fhMixDeltaPhiChargedAssocPtBin[bin]);
2892
2893         if(fFillEtaGapsHisto)
2894         {
2895           fhMixDeltaPhiChargedAssocPtBinDEta08[bin] = new TH2F(Form("hMixDeltaPhiDeltaEta0.8ChargedAssocPtBin%2.1f_%2.1f%s", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],sz.Data()),
2896                                                                Form("Mixed event #Delta #phi vs #it{p}_{T trigger} for associated #it{p}_{T} bin [%2.1f,%2.1f]%s, for #Delta #eta > 0.8", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],tz.Data()),
2897                                                                nptbins, ptmin, ptmax,  ndeltaphibins ,deltaphimin,deltaphimax);
2898           fhMixDeltaPhiChargedAssocPtBinDEta08[bin]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2899           fhMixDeltaPhiChargedAssocPtBinDEta08[bin]->SetYTitle("#Delta #phi (rad)");
2900
2901           fhMixDeltaPhiChargedAssocPtBinDEta0[bin] = new TH2F(Form("hMixDeltaPhiDeltaEta0ChargedAssocPtBin%2.1f_%2.1f%s", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],sz.Data()),
2902                                                               Form("Mixed event #Delta #phi vs #it{p}_{T trigger} for associated #it{p}_{T} bin [%2.1f,%2.1f]%s, for #Delta #eta = 0", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],tz.Data()),
2903                                                               nptbins, ptmin, ptmax,  ndeltaphibins ,deltaphimin,deltaphimax);
2904           fhMixDeltaPhiChargedAssocPtBinDEta0[bin]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2905           fhMixDeltaPhiChargedAssocPtBinDEta0[bin]->SetYTitle("#Delta #phi (rad)");
2906
2907           outputContainer->Add(fhMixDeltaPhiChargedAssocPtBinDEta08[bin]);
2908           outputContainer->Add(fhMixDeltaPhiChargedAssocPtBinDEta0[bin]);
2909         }
2910       }
2911     }
2912   }
2913
2914   return outputContainer;
2915
2916 }
2917
2918 //_____________________________________________________________________________________________________________________
2919 Bool_t AliAnaParticleHadronCorrelation::GetDecayPhotonMomentum(Int_t indexPhoton1, Int_t indexPhoton2, Int_t idetector)
2920 {
2921   // Get the momentum of the pi0/eta assigned decay photons
2922   // In case of pi0/eta trigger, we may want to check their decay correlation,
2923   // get their decay children
2924
2925   if(indexPhoton1!=-1 || indexPhoton2!=-1) return kFALSE;
2926
2927   AliDebug(1,Form("indexPhoton1 = %d, indexPhoton2 = %d", indexPhoton1, indexPhoton2));
2928
2929   TObjArray * clusters  = 0x0 ;
2930   if(idetector==kEMCAL) clusters = GetEMCALClusters() ;
2931   else                  clusters = GetPHOSClusters()  ;
2932
2933   for(Int_t iclus = 0; iclus < clusters->GetEntriesFast(); iclus++)
2934   {
2935     AliVCluster * photon =  (AliVCluster*) (clusters->At(iclus));
2936
2937     if(photon->GetID()==indexPhoton1) photon->GetMomentum(fDecayMom1,GetVertex(0)) ;
2938     if(photon->GetID()==indexPhoton2) photon->GetMomentum(fDecayMom2,GetVertex(0)) ;
2939
2940     AliDebug(1,Form("Photon1 = %f, Photon2 = %f", fDecayMom1.Pt(), fDecayMom2.Pt()));
2941
2942   } //cluster loop
2943
2944   return kTRUE;
2945
2946 }
2947
2948 //_____________________________________________________________
2949 Int_t AliAnaParticleHadronCorrelation::GetMCTagHistogramIndex(Int_t mcTag)
2950 {
2951   // Index of MC histograms depending on MC origin
2952
2953   if     ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPrompt) ||
2954            GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCFragmentation)) return 0;
2955   else if( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0)          ) return 1;
2956   else if( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0Decay)     ) return 2;
2957   else if( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEta)          ) return 3;
2958   else if( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEtaDecay)     ) return 4;
2959   else if( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCOtherDecay)   ) return 5;
2960   else if(!GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCElectron)     ) return 6;
2961   else                                                                                    return 7;
2962
2963 }
2964
2965 //_________________________________________
2966 void AliAnaParticleHadronCorrelation::Init()
2967 {
2968   //Init
2969   //Do some checks
2970
2971   if(!GetReader()->IsCTSSwitchedOn())
2972     AliFatal("STOP!: You want to use CTS tracks in analysis but not read!! \n!!Check the configuration file!!");
2973 }
2974
2975 //____________________________________________________
2976 void AliAnaParticleHadronCorrelation::InitParameters()
2977 {
2978
2979   //Initialize the parameters of the analysis.
2980   SetInputAODName("Particle");
2981   SetAODObjArrayName("Hadrons");
2982   AddToHistogramsName("AnaHadronCorr_");
2983
2984   SetPtCutRange(0.,300);
2985   fDeltaPhiMinCut       = TMath::DegToRad()*120.;
2986   fDeltaPhiMaxCut       = TMath::DegToRad()*240. ;
2987   fSelectIsolated       = kFALSE;
2988   fMakeSeveralUE        = kFALSE;
2989   fUeDeltaPhiMinCut     = TMath::DegToRad()*60.;
2990   fUeDeltaPhiMaxCut     = TMath::DegToRad()*120 ;
2991
2992   fNeutralCorr          = kFALSE ;
2993   fPi0Trigger           = kFALSE ;
2994   fDecayTrigger         = kFALSE ;
2995   fHMPIDCorrelation     = kFALSE ;
2996
2997   fMakeAbsoluteLeading  = kTRUE;
2998   fMakeNearSideLeading  = kFALSE;
2999
3000   fNAssocPtBins         = 9   ;
3001   fAssocPtBinLimit[0]   = 0.2 ;
3002   fAssocPtBinLimit[1]   = 0.5 ;
3003   fAssocPtBinLimit[2]   = 1.0 ;
3004   fAssocPtBinLimit[3]   = 2.0 ;
3005   fAssocPtBinLimit[4]   = 3.0 ;
3006   fAssocPtBinLimit[5]   = 4.0 ;
3007   fAssocPtBinLimit[6]   = 5.0 ;
3008   fAssocPtBinLimit[7]   = 6.0 ;
3009   fAssocPtBinLimit[8]   = 7.0 ;
3010   fAssocPtBinLimit[9]   = 8.0 ;
3011   fAssocPtBinLimit[10]  = 9.0 ;
3012   fAssocPtBinLimit[11]  = 10.0 ;
3013   fAssocPtBinLimit[12]  = 12.0 ;
3014   fAssocPtBinLimit[13]  = 14.0 ;
3015   fAssocPtBinLimit[14]  = 16.0 ;
3016   fAssocPtBinLimit[15]  = 20.0 ;
3017   fAssocPtBinLimit[16]  = 30.0 ;
3018   fAssocPtBinLimit[17]  = 40.0 ;
3019   fAssocPtBinLimit[18]  = 50.0 ;
3020   fAssocPtBinLimit[19]  = 200.0 ;
3021
3022   fUseMixStoredInReader = kTRUE;
3023
3024   fM02MinCut   = -1 ;
3025   fM02MaxCut   = -1 ;
3026
3027   fSelectLeadingHadronAngle = kFALSE;
3028   fFillLeadHadOppositeHisto = kFALSE;
3029   fMinLeadHadPhi = 150*TMath::DegToRad();
3030   fMaxLeadHadPhi = 210*TMath::DegToRad();
3031
3032   fMinLeadHadPt  = 1;
3033   fMaxLeadHadPt  = 100;
3034
3035   fMCGenTypeMin =  0;
3036   fMCGenTypeMax = 10;
3037
3038   fNDecayBits = 1;
3039   fDecayBits[0] = AliNeutralMesonSelection::kPi0;
3040   fDecayBits[1] = AliNeutralMesonSelection::kEta;
3041   fDecayBits[2] = AliNeutralMesonSelection::kPi0Side;
3042   fDecayBits[3] = AliNeutralMesonSelection::kEtaSide;
3043 }
3044
3045 //_________________________________________________________________________
3046 Bool_t AliAnaParticleHadronCorrelation::IsTriggerTheEventLeadingParticle()
3047 {
3048   // Check if the what of the selected triggers is leading particle comparing
3049   // with all the triggers, all the tracks or all the clusters (if requested for the clusters).
3050
3051   Double_t ptTrig      = GetMinPt();
3052   Double_t phiTrig     = 0 ;
3053   fLeadingTriggerIndex =-1 ;
3054   Int_t index          =-1 ;
3055   AliAODPWG4ParticleCorrelation* pLeading = 0;
3056
3057   // Loop on stored AOD particles, find leading trigger on the selected list, with at least min pT.
3058
3059   for(Int_t iaod = 0; iaod < GetInputAODBranch()->GetEntriesFast() ; iaod++)
3060   {
3061     AliAODPWG4ParticleCorrelation* particle =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
3062     particle->SetLeadingParticle(kFALSE); // set it later
3063
3064     // Vertex cut in case of mixing
3065     Int_t check = CheckMixedEventVertex(particle->GetCaloLabel(0), particle->GetTrackLabel(0));
3066     if(check ==  0) continue;
3067     if(check == -1) return kFALSE; // not sure if it is correct.
3068
3069     // find the leading particles with highest momentum
3070     if (particle->Pt() > ptTrig)
3071     {
3072       ptTrig   = particle->Pt() ;
3073       phiTrig  = particle->Phi();
3074       index    = iaod     ;
3075       pLeading = particle ;
3076     }
3077   }// finish search of leading trigger particle on the AOD branch.
3078
3079   if(index < 0) return kFALSE;
3080
3081   //printf("AOD leading pT %2.2f, ID %d\n", pLeading->Pt(),pLeading->GetCaloLabel(0));
3082
3083   if(phiTrig < 0 ) phiTrig += TMath::TwoPi();
3084
3085   // Compare if it is the leading of all tracks
3086
3087   for(Int_t ipr = 0;ipr < GetCTSTracks()->GetEntriesFast() ; ipr ++ )
3088   {
3089     AliVTrack * track = (AliVTrack *) (GetCTSTracks()->At(ipr)) ;
3090
3091     if(track->GetID() == pLeading->GetTrackLabel(0) || track->GetID() == pLeading->GetTrackLabel(1) ||
3092        track->GetID() == pLeading->GetTrackLabel(2) || track->GetID() == pLeading->GetTrackLabel(3)   ) continue ;
3093
3094     fTrackVector.SetXYZ(track->Px(),track->Py(),track->Pz());
3095     Float_t pt   = fTrackVector.Pt();
3096     Float_t phi  = fTrackVector.Phi() ;
3097     if(phi < 0) phi+=TMath::TwoPi();
3098
3099     //jump out this event if near side associated particle pt larger than trigger
3100     if (fMakeNearSideLeading)
3101     {
3102       Float_t deltaPhi = phiTrig-phi;
3103       if(deltaPhi <= -TMath::PiOver2()) deltaPhi+=TMath::TwoPi();
3104       if(deltaPhi > 3*TMath::PiOver2()) deltaPhi-=TMath::TwoPi();
3105
3106       if(pt > ptTrig && deltaPhi < TMath::PiOver2())  return kFALSE;
3107     }
3108     //jump out this event if there is any other particle with pt larger than trigger
3109     else
3110     {
3111       if(pt > ptTrig)  return kFALSE ;
3112     }
3113    }// track loop
3114
3115   // Compare if it is leading of all calorimeter clusters
3116
3117   if(fCheckLeadingWithNeutralClusters)
3118   {
3119     // Select the calorimeter cluster list
3120     TObjArray * nePl = 0x0;
3121     if      (pLeading->GetDetectorTag() == kPHOS )
3122       nePl = GetPHOSClusters();
3123     else
3124       nePl = GetEMCALClusters();
3125
3126     if(!nePl) return kTRUE; // Do the selection just with the tracks if no calorimeter is available.
3127
3128     for(Int_t ipr = 0;ipr < nePl->GetEntriesFast() ; ipr ++ )
3129     {
3130       AliVCluster * cluster = (AliVCluster *) (nePl->At(ipr)) ;
3131
3132       if(cluster->GetID() == pLeading->GetCaloLabel(0) || cluster->GetID() == pLeading->GetCaloLabel(1) ) continue ;
3133
3134       cluster->GetMomentum(fMomentum,GetVertex(0));
3135
3136       Float_t pt   = fMomentum.Pt();
3137       Float_t phi  = fMomentum.Phi() ;
3138       if(phi < 0) phi+=TMath::TwoPi();
3139
3140       if(IsTrackMatched(cluster,GetReader()->GetInputEvent())) continue ; // avoid charged clusters, already covered by tracks, or cluster merging with track.
3141
3142       //jump out this event if near side associated particle pt larger than trigger
3143       // not really needed for calorimeter, unless DCal is included
3144       if (fMakeNearSideLeading)
3145       {
3146         Float_t deltaPhi = phiTrig-phi;
3147         if(deltaPhi <= -TMath::PiOver2()) deltaPhi+=TMath::TwoPi();
3148         if(deltaPhi > 3*TMath::PiOver2()) deltaPhi-=TMath::TwoPi();
3149
3150         if(pt > ptTrig && deltaPhi < TMath::PiOver2()) return kFALSE ;
3151       }
3152       //jump out this event if there is any other particle with pt larger than trigger
3153       else
3154       {
3155         if(pt > ptTrig)  return kFALSE ;
3156       }
3157     }// cluster loop
3158   } // check neutral clusters
3159
3160   fLeadingTriggerIndex = index ;
3161   pLeading->SetLeadingParticle(kTRUE);
3162
3163   AliDebug(1,Form("\t particle AOD with index %d is leading with pT %2.2f", fLeadingTriggerIndex, pLeading->Pt()));
3164
3165   return kTRUE;
3166
3167 }
3168
3169 //_________________________________________________________________
3170 void  AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms()
3171 {
3172   //Particle-Hadron Correlation Analysis, fill histograms
3173
3174   if(!GetInputAODBranch())
3175   {
3176     AliFatal(Form("No input particles in AOD with name branch < %s >, STOP",GetInputAODName().Data()));
3177     return ; // coverity
3178   }
3179
3180   Int_t naod = GetInputAODBranch()->GetEntriesFast();
3181   if( naod == 0 )
3182   {
3183     AliWarning("No particle AOD found!");
3184     return ; // no trigger particles found.
3185   }
3186
3187   AliDebug(1,Form("Begin hadron correlation analysis, fill histograms"));
3188   AliDebug(1,Form("n particle branch aod entries %d", naod));
3189   AliDebug(1,Form("In CTS aod entries %d",GetCTSTracks()->GetEntriesFast()));
3190
3191   //------------------------------------------------------
3192   // Find leading trigger if analysis request only leading,
3193   // if there is no leading trigger, then skip the event
3194
3195   Int_t iaod = 0 ;
3196   if( fMakeAbsoluteLeading || fMakeNearSideLeading )
3197   {
3198     Bool_t leading = IsTriggerTheEventLeadingParticle();
3199
3200     AliDebug(1,Form("AOD Leading trigger? %d, with index %d",leading,fLeadingTriggerIndex));
3201
3202     if(!leading)
3203     {
3204       AliDebug(1,"Leading was requested and not found");
3205       return ;
3206     }
3207     else
3208     {
3209       // Select only the leading in the trigger AOD loop
3210       naod = fLeadingTriggerIndex+1 ;
3211       iaod = fLeadingTriggerIndex   ;
3212     }
3213   }
3214
3215   //------------------------------------------------------
3216   // Get event multiplicity and bins
3217
3218   Float_t cen         = GetEventCentrality();
3219   Float_t ep          = GetEventPlaneAngle();
3220   if(IsHighMultiplicityAnalysisOn()) fhTriggerEventPlaneCentrality->Fill(cen,ep);
3221
3222   Int_t   mixEventBin = GetEventMixBin();
3223   Int_t   vzbin       = GetEventVzBin();
3224
3225   //------------------------------------------------------
3226   // Loop on trigger AOD
3227
3228   for( iaod = 0; iaod < naod; iaod++ )
3229   {
3230     AliAODPWG4ParticleCorrelation* particle =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
3231
3232     //
3233     // Trigger particle selection criteria:
3234     //
3235     Float_t pt = particle->Pt();
3236
3237     if(pt < GetMinPt() || pt > GetMaxPt() ) continue ;
3238
3239     fhPtTriggerInput->Fill(pt);
3240
3241     //
3242     // check if it was a calorimeter cluster
3243     // and if the shower shape cut was requested apply it.
3244     // Not needed if already done at the particle identification level,
3245     // but for isolation studies, it is preferred not to remove so we do it here
3246     //
3247
3248     AliDebug(1,Form("%s Trigger : min %f, max %f, det %d",
3249                     GetInputAODName().Data(),fM02MinCut,fM02MaxCut,particle->GetDetectorTag()));
3250
3251     if(fM02MaxCut > 0 && fM02MinCut > 0) //clID1 > 0 && clID2 < 0 &&
3252     {
3253 //      Int_t iclus = -1;
3254 //      TObjArray* clusters = 0x0;
3255 //      if     (particle->GetDetectorTag() == kEMCAL) clusters = GetEMCALClusters();
3256 //      else if(particle->GetDetectorTag() == kPHOS ) clusters = GetPHOSClusters();
3257 //
3258 //      if(clusters)
3259 //      {
3260 //        AliVCluster *cluster = FindCluster(clusters,clID1,iclus);
3261 //        Float_t m02 = cluster->GetM02();
3262 //        if(m02 > fM02MaxCut || m02 < fM02MinCut) continue ;
3263 //      }
3264
3265       Float_t m02 = particle->GetM02();
3266
3267       if(m02 > fM02MaxCut || m02 < fM02MinCut) continue ;
3268
3269       fhPtTriggerSSCut->Fill(pt);
3270
3271       AliDebug(1,"Pass the shower shape cut");
3272     }
3273
3274     //
3275     // Check if the particle is isolated or if we want to take the isolation into account
3276     // This bool is set in AliAnaParticleIsolation
3277     //
3278     if(OnlyIsolated())
3279     {
3280       if( !particle->IsIsolated() ) continue;
3281
3282       fhPtTriggerIsoCut->Fill(pt);
3283
3284       AliDebug(1,"Pass the isolation cut");
3285     }
3286
3287     //
3288     // Check if trigger is in fiducial region
3289     //
3290     if(IsFiducialCutOn())
3291     {
3292       Bool_t in = GetFiducialCut()->IsInFiducialCut(particle->Eta(),particle->Phi(),particle->GetDetectorTag()) ;
3293
3294       if(! in ) continue ;
3295
3296       AliDebug(1,"Pass the fiducial cut");
3297     }
3298
3299     fhPtTriggerFidCut->Fill(pt);
3300
3301     //---------------------------------------
3302     // Make correlation
3303
3304     // Find the leading hadron in the opposite hemisphere to the triggeer
3305     // and accept the trigger if leading is in defined window.
3306     Bool_t okLeadHad = kTRUE;
3307     if(fSelectLeadingHadronAngle || fFillLeadHadOppositeHisto)
3308     {
3309       okLeadHad = FindLeadingOppositeHadronInWindow(particle);
3310       if(!okLeadHad && fSelectLeadingHadronAngle) continue;
3311     }
3312
3313     //
3314     // Charged particles correlation
3315     //
3316     MakeChargedCorrelation(particle);
3317
3318     // MC
3319     Int_t mcIndex = -1;
3320     Int_t mcTag   = particle->GetTag();
3321     Bool_t lostDecayPair = kFALSE;
3322     if(IsDataMC())
3323     {
3324       mcIndex = GetMCTagHistogramIndex(mcTag);
3325       lostDecayPair = GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCDecayPairLost);
3326       MakeMCChargedCorrelation(particle->GetLabel(), mcIndex,lostDecayPair);
3327     }
3328
3329     // Do own mixed event with charged,
3330     // add event and remove previous or fill the mixed histograms
3331     if(DoOwnMix())
3332       MakeChargedMixCorrelation(particle);
3333
3334     //
3335     // Neutral particles correlation
3336     //
3337     if(fNeutralCorr)
3338       MakeNeutralCorrelation(particle);
3339
3340     //----------------------------------------------------------------
3341     // Fill trigger pT related histograms if not absolute leading
3342
3343     //
3344     // pT of the trigger, vs trigger origin if MC
3345     //
3346     fhPtTrigger->Fill(pt);
3347     if(IsDataMC() && mcIndex >=0 && mcIndex < fgkNmcTypes)
3348     {
3349       fhPtTriggerMC[mcIndex]->Fill(pt);
3350       if( lostDecayPair )
3351       {
3352         // check index in GetMCTagIndexHistogram
3353         if     ( mcIndex == 2 ) fhPtTriggerMC[8]->Fill(pt); // pi0 decay
3354         else if( mcIndex == 4 ) fhPtTriggerMC[9]->Fill(pt); // eta decay
3355       }
3356     }
3357
3358     if(fDecayTrigger)
3359     {
3360       Int_t decayTag = particle->DecayTag();
3361       if(decayTag < 0) decayTag = 0;
3362
3363       for(Int_t ibit = 0; ibit<fNDecayBits; ibit++)
3364       {
3365         if(GetNeutralMesonSelection()->CheckDecayBit(decayTag,fDecayBits[ibit]))
3366         {
3367           fhPtDecayTrigger[ibit]->Fill(pt);
3368
3369           if(IsDataMC() && mcIndex >=0 && mcIndex < fgkNmcTypes)
3370           {
3371             fhPtDecayTriggerMC[ibit][mcIndex]->Fill(pt);
3372             if( lostDecayPair )
3373             {
3374               // check index in GetMCTagIndexHistogram
3375               if( mcIndex == 2 )  fhPtDecayTriggerMC[ibit][8]->Fill(pt); // pi0 decay
3376               if( mcIndex == 4 )  fhPtDecayTriggerMC[ibit][9]->Fill(pt); // eta decay
3377             }
3378           }
3379         }// check bit
3380       }// bit loop
3381     }
3382
3383     //
3384     // Acceptance of the trigger
3385     //
3386     Float_t phi = particle->Phi();
3387     if( phi < 0 ) phi+=TMath::TwoPi();
3388     fhPhiTrigger->Fill(pt, phi);
3389
3390     fhEtaTrigger->Fill(pt, particle->Eta());
3391     //printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms() - Trigger particle : pt %f, eta %f, phi %f\n",particle->Pt(),particle->Eta(),phi);
3392
3393     //----------------------------------
3394     // Trigger particle pT vs event bins
3395
3396     fhPtTriggerBin->Fill(pt,mixEventBin);
3397     if(fCorrelVzBin)
3398       fhPtTriggerVzBin->Fill(pt,vzbin);
3399
3400     if(IsHighMultiplicityAnalysisOn())
3401     {
3402       fhPtTriggerCentrality->Fill(pt,cen);
3403       fhPtTriggerEventPlane->Fill(pt,ep);
3404     }
3405
3406     //----------------------------------
3407     // Trigger particle pT vs pile-up
3408
3409     if(IsPileUpAnalysisOn())
3410     {
3411       Int_t vtxBC = GetReader()->GetVertexBC();
3412       if(vtxBC == 0 || vtxBC==AliVTrack::kTOFBCNA)     fhPtTriggerVtxBC0->Fill(pt);
3413
3414       if(GetReader()->IsPileUpFromSPD())               fhPtTriggerPileUp[0]->Fill(pt);
3415       if(GetReader()->IsPileUpFromEMCal())             fhPtTriggerPileUp[1]->Fill(pt);
3416       if(GetReader()->IsPileUpFromSPDOrEMCal())        fhPtTriggerPileUp[2]->Fill(pt);
3417       if(GetReader()->IsPileUpFromSPDAndEMCal())       fhPtTriggerPileUp[3]->Fill(pt);
3418       if(GetReader()->IsPileUpFromSPDAndNotEMCal())    fhPtTriggerPileUp[4]->Fill(pt);
3419       if(GetReader()->IsPileUpFromEMCalAndNotSPD())    fhPtTriggerPileUp[5]->Fill(pt);
3420       if(GetReader()->IsPileUpFromNotSPDAndNotEMCal()) fhPtTriggerPileUp[6]->Fill(pt);
3421     }
3422   } // AOD trigger loop
3423
3424   //Reinit for next event
3425   fLeadingTriggerIndex = -1;
3426
3427   AliDebug(1,"End fill histograms");
3428 }
3429
3430 //_______________________________________________________________________________________________________
3431 void  AliAnaParticleHadronCorrelation::MakeChargedCorrelation(AliAODPWG4ParticleCorrelation *aodParticle)
3432 {
3433   // Charged Hadron Correlation Analysis
3434   AliDebug(1,"Make trigger particle - charged hadron correlation");
3435
3436   Float_t phiTrig = aodParticle->Phi();
3437   Float_t etaTrig = aodParticle->Eta();
3438   Float_t ptTrig  = aodParticle->Pt();
3439   Int_t    mcTag  = aodParticle->GetTag();
3440   Double_t bz     = GetReader()->GetInputEvent()->GetMagneticField();
3441
3442
3443   Int_t   decayTag = 0;
3444   if(fDecayTrigger)
3445   {
3446     //decay = aodParticle->IsTagged();
3447     decayTag = aodParticle->DecayTag();
3448     if(decayTag < 0) decayTag = 0;
3449 //    printf("Correlation: pT %2.2f, BTag %d, Tagged %d\n",ptTrig, decayTag, aodParticle->IsTagged());
3450 //    printf("\t check bit Pi0 %d, Eta %d,  Pi0Side %d, EtaSide %d\n",
3451 //           GetNeutralMesonSelection()->CheckDecayBit(decayTag,AliNeutralMesonSelection::kPi0),
3452 //           GetNeutralMesonSelection()->CheckDecayBit(decayTag,AliNeutralMesonSelection::kEta),
3453 //           GetNeutralMesonSelection()->CheckDecayBit(decayTag,AliNeutralMesonSelection::kPi0Side),
3454 //           GetNeutralMesonSelection()->CheckDecayBit(decayTag,AliNeutralMesonSelection::kEtaSide));
3455   }
3456
3457   Float_t pt       = -100. ;
3458   Float_t phi      = -100. ;
3459   Float_t eta      = -100. ;
3460   Float_t deltaPhi = -100. ;
3461
3462   TObjArray * reftracks = 0x0;
3463   Int_t nrefs           = 0;
3464
3465   // Mixed event settings
3466   Int_t evtIndex11   = -1 ; // cluster trigger or pi0 trigger
3467   Int_t evtIndex12   = -1 ; // pi0 trigger
3468   Int_t evtIndex13   = -1 ; // charged trigger
3469
3470   if (GetMixedEvent())
3471   {
3472     evtIndex11 = GetMixedEvent()->EventIndexForCaloCluster(aodParticle->GetCaloLabel(0)) ;
3473     evtIndex12 = GetMixedEvent()->EventIndexForCaloCluster(aodParticle->GetCaloLabel(1)) ;
3474     evtIndex13 = GetMixedEvent()->EventIndex(aodParticle->GetTrackLabel(0)) ;
3475   }
3476
3477   // Track multiplicity or cent bin
3478   Int_t cenbin = 0;
3479   if(IsHighMultiplicityAnalysisOn()) cenbin = GetEventCentralityBin();
3480
3481   //
3482   // In case of pi0/eta trigger, we may want to check their decay correlation,
3483   // get their decay children
3484   //
3485
3486   Bool_t decayFound = kFALSE;
3487   if( fPi0Trigger )
3488   {
3489     decayFound = GetDecayPhotonMomentum(aodParticle->GetCaloLabel(0),aodParticle->GetCaloLabel(1),aodParticle->GetDetectorTag());
3490     if(decayFound)
3491     {
3492       fhPtPi0DecayRatio->Fill(ptTrig, fDecayMom1.Pt()/ptTrig);
3493       fhPtPi0DecayRatio->Fill(ptTrig, fDecayMom2.Pt()/ptTrig);
3494     }
3495   }
3496
3497   //-----------------------------------------------------------------------
3498   // Track loop, select tracks with good pt, phi and fill AODs or histograms
3499   //-----------------------------------------------------------------------
3500
3501   for(Int_t ipr = 0;ipr < GetCTSTracks()->GetEntriesFast() ; ipr ++ )
3502   {
3503     AliVTrack * track = (AliVTrack *) (GetCTSTracks()->At(ipr)) ;
3504
3505     fTrackVector.SetXYZ(track->Px(),track->Py(),track->Pz());
3506     pt   = fTrackVector.Pt();
3507     eta  = fTrackVector.Eta();
3508     phi  = fTrackVector.Phi() ;
3509     if(phi < 0) phi+=TMath::TwoPi();
3510
3511     //Select only hadrons in pt range
3512     if(pt < fMinAssocPt || pt > fMaxAssocPt) continue ;
3513
3514     //remove trigger itself for correlation when use charged triggers
3515     if( track->GetID() == aodParticle->GetTrackLabel(0) || track->GetID() == aodParticle->GetTrackLabel(1) ||
3516         track->GetID() == aodParticle->GetTrackLabel(2) || track->GetID() == aodParticle->GetTrackLabel(3)   )
3517       continue ;
3518
3519     //Only for mixed event frame
3520     Int_t evtIndex2 = 0 ;
3521     if (GetMixedEvent())
3522     {
3523       evtIndex2 = GetMixedEvent()->EventIndex(track->GetID()) ;
3524       if (evtIndex11 == evtIndex2 || evtIndex12 == evtIndex2 || evtIndex13 == evtIndex2 ) // photon and track from different events
3525         continue ;
3526       //vertex cut
3527       if (TMath::Abs(GetVertex(evtIndex2)[2]) > GetZvertexCut())
3528         continue;
3529     }
3530
3531      AliDebug(2,Form("Selected charge for momentum imbalance: pt %2.2f, phi %2.2f, eta %2.2f",pt,phi,eta));
3532
3533     // ------------------------------
3534     // Track type bin or bits setting
3535     //
3536
3537     //
3538     // * Set the pt associated bin for the defined bins *
3539     //
3540     Int_t assocBin   = -1;
3541     for(Int_t i = 0 ; i < fNAssocPtBins ; i++)
3542     {
3543       if(pt > fAssocPtBinLimit[i] && pt < fAssocPtBinLimit[i+1]) assocBin= i;
3544     }
3545
3546     //
3547     // * Assign to the histogram array a bin corresponding
3548     // to a combination of pTa and vz bins *
3549     //
3550     Int_t nz = 1;
3551     Int_t vz = 0;
3552
3553     if(fCorrelVzBin)
3554     {
3555       nz = GetNZvertBin();
3556       vz = GetEventVzBin();
3557     }
3558
3559     Int_t bin = assocBin*nz+vz;
3560
3561     //printf("assoc Bin = %d, vZ bin  = %d, bin = %d \n", assocBin,GetEventVzBin(),bin);
3562
3563     //
3564     // * Get the status of the TOF bit *
3565     //
3566     ULong_t status = track->GetStatus();
3567     Bool_t okTOF = ( (status & AliVTrack::kTOFout) == AliVTrack::kTOFout ) ;
3568     //Double32_t tof = track->GetTOFsignal()*1e-3;
3569     Int_t trackBC = track->GetTOFBunchCrossing(bz);
3570
3571     Int_t outTOF = -1;
3572     if     (okTOF && trackBC!=0) outTOF = 1;
3573     else if(okTOF && trackBC==0) outTOF = 0;
3574
3575     //----------------
3576     // Fill Histograms
3577
3578     //
3579     // Azimuthal Angle histograms
3580     //
3581
3582     deltaPhi = phiTrig-phi;
3583
3584     //
3585     // Calculate deltaPhi shift so that for the particles on the opposite side
3586     // it is defined between 90 and 270 degrees
3587     // Shift [-360,-90]  to [0, 270]
3588     // and [270,360] to [-90,0]
3589     if(deltaPhi <= -TMath::PiOver2()) deltaPhi+=TMath::TwoPi();
3590     if(deltaPhi > 3*TMath::PiOver2()) deltaPhi-=TMath::TwoPi();
3591
3592     FillChargedAngularCorrelationHistograms(pt,  ptTrig,  bin, phi, phiTrig,  deltaPhi,
3593                                             eta, etaTrig, decayTag, track->GetHMPIDsignal(),
3594                                             outTOF, cenbin, mcTag);
3595
3596     //
3597     // Imbalance zT/xE/pOut histograms
3598     //
3599
3600     //
3601     // Delta phi cut for momentum imbalance correlation
3602     //
3603     if  ( (deltaPhi > fDeltaPhiMinCut)   && (deltaPhi < fDeltaPhiMaxCut)   )
3604       FillChargedMomentumImbalanceHistograms(ptTrig, pt, deltaPhi, cenbin, track->Charge(),
3605                                              assocBin, decayTag, outTOF, mcTag);
3606
3607     //
3608     // Underlying event, right side, default case
3609     //
3610     if ( (deltaPhi > fUeDeltaPhiMinCut) && (deltaPhi < fUeDeltaPhiMaxCut) )
3611       FillChargedUnderlyingEventHistograms(ptTrig, pt, deltaPhi, cenbin, outTOF);
3612
3613     //
3614     // Several UE calculation,  in different perpendicular regions, up to 6:
3615     // left, right, upper-left, lower left, upper-right, lower-right
3616     //
3617     if(fMakeSeveralUE)
3618       FillChargedUnderlyingEventSidesHistograms(ptTrig,pt,deltaPhi);
3619
3620     //
3621     if(fPi0Trigger && decayFound)
3622       FillDecayPhotonCorrelationHistograms(pt, phi, kTRUE) ;
3623
3624     //
3625     // Add track reference to array
3626     //
3627     if(fFillAODWithReferences)
3628     {
3629       nrefs++;
3630       if(nrefs==1)
3631       {
3632         reftracks = new TObjArray(0);
3633         TString trackname = Form("%sTracks", GetAODObjArrayName().Data());
3634         reftracks->SetName(trackname.Data());
3635         reftracks->SetOwner(kFALSE);
3636       }
3637
3638       reftracks->Add(track);
3639     }// reference track to AOD
3640   }// track loop
3641
3642   //Fill AOD with reference tracks, if not filling histograms
3643   if(fFillAODWithReferences && reftracks)
3644   {
3645     aodParticle->AddObjArray(reftracks);
3646   }
3647
3648 }
3649
3650 //_________________________________________________________________________________________________________
3651 void AliAnaParticleHadronCorrelation::MakeChargedMixCorrelation(AliAODPWG4ParticleCorrelation *aodParticle)
3652 {
3653   // Mix current trigger with tracks in another MB event
3654
3655   AliDebug(1,Form("Make trigger particle - charged hadron mixed event correlation"));
3656
3657   if(GetMixedEvent()) return;  // This is not the mixed event from general mixing frame
3658
3659   // Get the event with similar caracteristics
3660   //printf("MakeChargedMixCorrelation for %s\n",GetInputAODName().Data());
3661
3662   AliAnalysisManager   * manager      = AliAnalysisManager::GetAnalysisManager();
3663
3664   AliInputEventHandler * inputHandler = dynamic_cast<AliInputEventHandler*>(manager->GetInputEventHandler());
3665
3666   if(!inputHandler) return;
3667
3668   if(!(inputHandler->IsEventSelected( ) & GetReader()->GetEventTriggerMask())) return;
3669
3670   // Get the pool, check if it exits
3671   Int_t eventBin = GetEventMixBin();
3672
3673   //Check that the bin exists, if not (bad determination of RP, ntrality or vz bin) do nothing
3674   if(eventBin < 0) return;
3675
3676   fhEventBin->Fill(eventBin);
3677
3678   // get neutral clusters pool?
3679   Bool_t isoCase = OnlyIsolated() && (GetIsolationCut()->GetParticleTypeInCone() != AliIsolationCut::kOnlyCharged);
3680   Bool_t neutralMix = fFillNeutralEventMixPool || isoCase ;
3681
3682   TList * pool     = 0;
3683   TList * poolCalo = 0;
3684   if(fUseMixStoredInReader)
3685   {
3686     pool     = GetReader()->GetListWithMixedEventsForTracks(eventBin);
3687     if(neutralMix) poolCalo = GetReader()->GetListWithMixedEventsForCalo  (eventBin);
3688   }
3689   else
3690   {
3691     pool     = fListMixTrackEvents[eventBin];
3692     if(neutralMix) poolCalo = fListMixCaloEvents [eventBin];
3693   }
3694
3695   if(!pool) return ;
3696
3697   if( neutralMix && !poolCalo )
3698     AliWarning("Careful, cluster pool not available");
3699
3700   Double_t ptTrig  = aodParticle->Pt();
3701   Double_t etaTrig = aodParticle->Eta();
3702   Double_t phiTrig = aodParticle->Phi();
3703   if(phiTrig < 0.) phiTrig+=TMath::TwoPi();
3704
3705   AliDebug(1,Form("Pool bin %d size %d, trigger trigger pt=%f, phi=%f, eta=%f",
3706                   eventBin,pool->GetSize(), ptTrig,phiTrig,etaTrig));
3707
3708   Double_t ptAssoc  = -999.;
3709   Double_t phiAssoc = -999.;
3710   Double_t etaAssoc = -999.;
3711   Double_t deltaPhi = -999.;
3712   Double_t deltaEta = -999.;
3713   Double_t xE       = -999.;
3714
3715   // Start from first event in pool except if in this same event the pool was filled
3716   Int_t ev0 = 0;
3717   if(GetReader()->GetLastTracksMixedEvent() == GetEventNumber()) ev0 = 1;
3718
3719   for(Int_t ev=ev0; ev < pool->GetSize(); ev++)
3720   {
3721     //
3722     // Recover the lists of tracks or clusters
3723     //
3724     TObjArray* bgTracks = static_cast<TObjArray*>(pool->At(ev));
3725     TObjArray* bgCalo   = 0;
3726
3727     // Recover the clusters list if requested
3728     if( neutralMix && poolCalo )
3729     {
3730       if(pool->GetSize()!=poolCalo->GetSize())
3731         AliWarning("Different size of calo and track pools");
3732
3733       bgCalo = static_cast<TObjArray*>(poolCalo->At(ev));
3734
3735       if(!bgCalo) AliDebug(1,Form("Event %d in calo pool not available?",ev));
3736     }
3737
3738     //
3739     // Isolate the trigger in the mixed event with mixed tracks and clusters
3740     //
3741     if( OnlyIsolated() )
3742     {
3743       Int_t   n=0, nfrac = 0;
3744       Bool_t  isolated = kFALSE;
3745       Float_t coneptsum = 0, coneptlead = 0;
3746       GetIsolationCut()->MakeIsolationCut(bgTracks,bgCalo,
3747                                           GetReader(), GetCaloPID(),
3748                                           kFALSE, aodParticle, "",
3749                                           n,nfrac,coneptsum,coneptlead,isolated);
3750
3751       //printf("AliAnaParticleHadronCorrelation::MakeChargedMixCorrelation() - Isolated? %d - cone %f, ptthres %f",
3752       //       isolated,GetIsolationCut()->GetConeSize(),GetIsolationCut()->GetPtThreshold());
3753       //if(bgTracks)printf(" - n track %d", bgTracks->GetEntriesFast());
3754       //printf("\n");
3755
3756       if(!isolated) continue ;
3757     }
3758
3759     //
3760     // Check if the trigger is leading of mixed event
3761     //
3762     Int_t nTracks=bgTracks->GetEntriesFast();
3763
3764     if(fMakeNearSideLeading || fMakeAbsoluteLeading)
3765     {
3766       Bool_t leading = kTRUE;
3767       for(Int_t jlead = 0;jlead < nTracks; jlead++ )
3768       {
3769         AliAODPWG4Particle *track = (AliAODPWG4Particle*) bgTracks->At(jlead) ;
3770
3771         ptAssoc  = track->Pt();
3772         phiAssoc = track->Phi() ;
3773         if(phiAssoc < 0) phiAssoc+=TMath::TwoPi();
3774
3775         if (fMakeNearSideLeading)
3776         {
3777           deltaPhi = phiTrig-phiAssoc;
3778           if(deltaPhi <= -TMath::PiOver2()) deltaPhi+=TMath::TwoPi();
3779           if(deltaPhi > 3*TMath::PiOver2()) deltaPhi-=TMath::TwoPi();
3780
3781           if(ptAssoc > ptTrig && deltaPhi < TMath::PiOver2())
3782           {
3783             leading = kFALSE;
3784             break;
3785           }
3786         }
3787         //jump out this event if there is any other particle with pt larger than trigger
3788         else if(fMakeAbsoluteLeading)
3789         {
3790           if(ptAssoc > ptTrig)
3791           {
3792             leading = kFALSE;
3793             break;
3794           }
3795         }
3796       }
3797
3798       if( !neutralMix && fCheckLeadingWithNeutralClusters )
3799         AliWarning("Leading of clusters requested but no clusters in mixed event");
3800
3801       if(neutralMix && fCheckLeadingWithNeutralClusters && bgCalo)
3802       {
3803         Int_t nClusters=bgCalo->GetEntriesFast();
3804         for(Int_t jlead = 0;jlead <nClusters; jlead++ )
3805         {
3806           AliAODPWG4Particle *cluster= (AliAODPWG4Particle*) bgCalo->At(jlead) ;
3807
3808           ptAssoc  = cluster->Pt();
3809           phiAssoc = cluster->Phi() ;
3810           if(phiAssoc < 0) phiAssoc+=TMath::TwoPi();
3811
3812           if (fMakeNearSideLeading)
3813           {
3814             deltaPhi = phiTrig-phiAssoc;
3815             if(deltaPhi <= -TMath::PiOver2()) deltaPhi+=TMath::TwoPi();
3816             if(deltaPhi > 3*TMath::PiOver2()) deltaPhi-=TMath::TwoPi();
3817
3818             if(ptAssoc > ptTrig && deltaPhi < TMath::PiOver2())
3819             {
3820               leading = kFALSE;
3821               break;
3822             }
3823           }
3824           //jump out this event if there is any other particle with pt larger than trigger
3825           else if(fMakeAbsoluteLeading)
3826           {
3827             if(ptAssoc > ptTrig)
3828             {
3829               leading = kFALSE;
3830               break;
3831             }
3832           }
3833         }
3834       }
3835
3836       if(!leading) continue; // not leading, check the next event in pool
3837     }
3838
3839     //
3840     // Fill histograms for selected triggers
3841     //
3842
3843     fhEventMixBin->Fill(eventBin);
3844
3845     //printf("\t Read Pool event %d, nTracks %d\n",ev,nTracks);
3846
3847     fhPtTriggerMixed   ->Fill(ptTrig);
3848     fhPhiTriggerMixed  ->Fill(ptTrig, phiTrig);
3849     fhEtaTriggerMixed  ->Fill(ptTrig, etaTrig);
3850     fhPtTriggerMixedBin->Fill(ptTrig,eventBin);
3851     if(fCorrelVzBin)fhPtTriggerMixedVzBin->Fill(ptTrig, GetEventVzBin());
3852
3853     //
3854     // Correlation histograms
3855     //
3856     for(Int_t j1 = 0;j1 <nTracks; j1++ )
3857     {
3858       AliAODPWG4Particle *track = (AliAODPWG4Particle*) bgTracks->At(j1) ;
3859
3860       if(!track) continue;
3861
3862       ptAssoc  = track->Pt();
3863       etaAssoc = track->Eta();
3864       phiAssoc = track->Phi() ;
3865       if(phiAssoc < 0) phiAssoc+=TMath::TwoPi();
3866
3867       deltaPhi = phiTrig-phiAssoc;
3868       if(deltaPhi < -TMath::PiOver2())  deltaPhi+=TMath::TwoPi();
3869       if(deltaPhi > 3*TMath::PiOver2()) deltaPhi-=TMath::TwoPi();
3870       deltaEta = etaTrig-etaAssoc;
3871
3872       AliDebug(1,Form("deltaPhi= %f, deltaEta=%f",deltaPhi, deltaEta));
3873
3874       // Angular correlation
3875       fhMixDeltaPhiCharged        ->Fill(ptTrig,   deltaPhi);
3876       fhMixDeltaPhiDeltaEtaCharged->Fill(deltaPhi, deltaEta);
3877
3878       //
3879       // Momentum imbalance
3880       //
3881       if ( (deltaPhi > fDeltaPhiMinCut)   && (deltaPhi < fDeltaPhiMaxCut)   )
3882       {
3883         xE = -ptAssoc/ptTrig*TMath::Cos(deltaPhi); // -(px*pxTrig+py*pyTrig)/(ptTrig*ptTrig);
3884
3885         if(xE < 0.)
3886           AliWarning(Form("Careful!!, negative xE %2.2f for right UE cos(dPhi %2.2f) = %2.2f, check correlation dPhi limits %f to %f",
3887                           xE,deltaPhi*TMath::RadToDeg(),TMath::Cos(deltaPhi),fDeltaPhiMinCut*TMath::RadToDeg(),fDeltaPhiMaxCut*TMath::RadToDeg()));
3888
3889         fhMixXECharged->Fill(ptTrig,xE);
3890         if(xE > 0 ) fhMixHbpXECharged->Fill(ptTrig, TMath::Log(1./xE));
3891       }
3892
3893       //
3894       // Underlying event momentum imbalance
3895       //
3896       if ( (deltaPhi > fUeDeltaPhiMinCut) && (deltaPhi < fUeDeltaPhiMaxCut) )
3897       {
3898         //Underlying event region
3899         Double_t randomphi = gRandom->Uniform(fDeltaPhiMinCut,fDeltaPhiMaxCut);
3900         Double_t uexE = -(ptAssoc/ptTrig)*TMath::Cos(randomphi);
3901
3902         if(uexE < 0.)
3903           AliWarning(Form("Careful!!, negative xE %2.2f for left UE cos(dPhi %2.2f) = %2.2f, check correlation dPhi limits %f to %f",
3904                           uexE,randomphi*TMath::RadToDeg(),TMath::Cos(randomphi),fDeltaPhiMinCut*TMath::RadToDeg(),fDeltaPhiMaxCut*TMath::RadToDeg()));
3905
3906         fhMixXEUeCharged->Fill(ptTrig,uexE);
3907       }
3908
3909       // Set the pt associated bin for the defined bins
3910       Int_t assocBin   = -1;
3911       for(Int_t i = 0 ; i < fNAssocPtBins ; i++)
3912       {
3913         if(ptAssoc > fAssocPtBinLimit[i] && ptAssoc < fAssocPtBinLimit[i+1]) assocBin= i;
3914       }
3915
3916       //
3917       // Assign to the histogram array a bin corresponding to a combination of pTa and vz bins
3918       //
3919       Int_t nz = 1;
3920       Int_t vz = 0;
3921
3922       if(fCorrelVzBin)
3923       {
3924         nz = GetNZvertBin();
3925         vz = GetEventVzBin();
3926       }
3927
3928       Int_t bin = assocBin*nz+vz;
3929
3930       if(bin < 0) continue ; // this pt bin was not considered
3931
3932       fhMixDeltaPhiChargedAssocPtBin        [bin]->Fill(ptTrig,   deltaPhi);
3933       fhMixDeltaPhiDeltaEtaChargedAssocPtBin[bin]->Fill(deltaPhi, deltaEta);
3934
3935       if(fFillEtaGapsHisto)
3936       {
3937         if(TMath::Abs(deltaEta) > 0.8)
3938           fhMixDeltaPhiChargedAssocPtBinDEta08  [bin]->Fill(ptTrig, deltaPhi);
3939         if(TMath::Abs(deltaEta) < 0.01)
3940           fhMixDeltaPhiChargedAssocPtBinDEta0   [bin]->Fill(ptTrig, deltaPhi);
3941       }
3942
3943     } // track loop
3944   } // mixed event loop
3945 }
3946
3947
3948 //_______________________________________________________________________________________________________
3949 void AliAnaParticleHadronCorrelation::MakeNeutralCorrelation(AliAODPWG4ParticleCorrelation * aodParticle)
3950 {
3951   // Neutral Pion Correlation Analysis
3952
3953   TObjArray * pi0list = (TObjArray*) GetAODBranch(fPi0AODBranchName); // For the future, foresee more possible pi0 lists
3954   if(!pi0list) return ;
3955
3956   Int_t npi0 = pi0list->GetEntriesFast();
3957   if(npi0 == 0) return ;
3958
3959   AliDebug(1,Form("Particle - pi0 correlation, %d pi0's",npi0));
3960
3961   Int_t evtIndex11 = 0 ;
3962   Int_t evtIndex12 = 0 ;
3963   if (GetMixedEvent())
3964   {
3965     evtIndex11 = GetMixedEvent()->EventIndexForCaloCluster(aodParticle->GetCaloLabel(0)) ;
3966     evtIndex12 = GetMixedEvent()->EventIndexForCaloCluster(aodParticle->GetCaloLabel(1)) ;
3967   }
3968
3969   Float_t pt   = -100. ;
3970   Float_t zT   = -100. ;
3971   Float_t phi  = -100. ;
3972   Float_t eta  = -100. ;
3973   Float_t xE   = -100. ;
3974   Float_t hbpXE= -100. ;
3975   Float_t hbpZT= -100. ;
3976
3977   Float_t ptTrig  = aodParticle->Pt();
3978   Float_t phiTrig = aodParticle->Phi();
3979   Float_t etaTrig = aodParticle->Eta();
3980   Float_t deltaPhi= -100. ;
3981   Float_t deltaEta= -100. ;
3982
3983   // In case of pi0/eta trigger, we may want to check their decay correlation,
3984   // get their decay children
3985
3986   Bool_t decayFound = kFALSE;
3987   if(fPi0Trigger) decayFound = GetDecayPhotonMomentum(aodParticle->GetCaloLabel(0),aodParticle->GetCaloLabel(1),aodParticle->GetDetectorTag());
3988
3989   TObjArray * refpi0 = 0x0;
3990   Int_t nrefs        = 0;
3991
3992   //Loop on stored AOD pi0
3993
3994   for(Int_t iaod = 0; iaod < npi0 ; iaod++)
3995   {
3996     AliAODPWG4Particle* pi0 =  (AliAODPWG4Particle*) (pi0list->At(iaod));
3997
3998     Int_t evtIndex2 = 0 ;
3999     Int_t evtIndex3 = 0 ;
4000     if (GetMixedEvent())
4001     {
4002       evtIndex2 = GetMixedEvent()->EventIndexForCaloCluster(pi0->GetCaloLabel(0)) ;
4003       evtIndex3 = GetMixedEvent()->EventIndexForCaloCluster(pi0->GetCaloLabel(1)) ;
4004
4005       if (evtIndex11 == evtIndex2 || evtIndex12 == evtIndex2 ||
4006           evtIndex11 == evtIndex3 || evtIndex12 == evtIndex3) // trigger and pi0 are not from different events
4007         continue ;
4008     }
4009
4010     pt  = pi0->Pt();
4011
4012     if(pt < fMinAssocPt || pt > fMaxAssocPt) continue ;
4013
4014     //remove trigger itself for correlation when use charged triggers
4015     if(aodParticle->GetCaloLabel(0) >= 0 &&
4016        (pi0->GetCaloLabel(0) == aodParticle->GetCaloLabel(0) || pi0->GetCaloLabel(1) == aodParticle->GetCaloLabel(0))) continue ;
4017
4018     if( aodParticle->GetCaloLabel(1) >= 0 &&
4019        (pi0->GetCaloLabel(0) == aodParticle->GetCaloLabel(1) || pi0->GetCaloLabel(1) == aodParticle->GetCaloLabel(1))) continue ;
4020
4021     //
4022     // Angular correlations
4023     //
4024     phi      = pi0->Phi() ;
4025     eta      = pi0->Eta() ;
4026     deltaEta = etaTrig-eta;
4027     deltaPhi = phiTrig-phi;
4028     if(deltaPhi <= -TMath::PiOver2()) deltaPhi+=TMath::TwoPi();
4029     if(deltaPhi > 3*TMath::PiOver2()) deltaPhi-=TMath::TwoPi();
4030
4031     fhEtaNeutral     ->Fill(pt    , eta     );
4032     fhPhiNeutral     ->Fill(pt    , phi     );
4033     fhDeltaEtaNeutral->Fill(ptTrig, deltaEta);
4034     fhDeltaPhiNeutral->Fill(ptTrig, deltaPhi);
4035
4036     if(pt > 2 ) fhDeltaPhiDeltaEtaNeutral->Fill(deltaPhi, deltaEta);
4037
4038     //
4039     // Momentum imbalance
4040     //
4041     zT  = pt/ptTrig ;
4042
4043     hbpZT = -100;
4044     hbpXE = -100;
4045
4046     if(zT > 0 ) hbpZT = TMath::Log(1./zT);
4047
4048     //delta phi cut for correlation
4049     if( (deltaPhi > fDeltaPhiMinCut) && ( deltaPhi < fDeltaPhiMaxCut) )
4050     {
4051       xE  =-pt/ptTrig*TMath::Cos(deltaPhi); // -(px*pxTrig+py*pyTrig)/(ptTrig*ptTrig);
4052
4053       if(xE < 0.)
4054         AliWarning(Form("Careful!!, negative xE %2.2f for right UE cos(dPhi %2.2f) = %2.2f, check correlation dPhi limits %f to %f",
4055                         xE,deltaPhi*TMath::RadToDeg(),TMath::Cos(deltaPhi),fDeltaPhiMinCut*TMath::RadToDeg(),fDeltaPhiMaxCut*TMath::RadToDeg()));
4056
4057       if(xE > 0 ) hbpXE = TMath::Log(1./xE);
4058
4059       fhDeltaPhiNeutralPt->Fill(pt,deltaPhi);
4060       fhXENeutral        ->Fill(ptTrig,xE);
4061       fhPtHbpXENeutral   ->Fill(ptTrig,hbpXE);
4062       fhZTNeutral        ->Fill(ptTrig,zT);
4063       fhPtHbpZTNeutral   ->Fill(ptTrig,hbpZT);
4064     }
4065     else if ( (deltaPhi > fUeDeltaPhiMinCut) && (deltaPhi < fUeDeltaPhiMaxCut) )
4066     {
4067       // Randomize angle for xE calculation
4068       Double_t randomphi = gRandom->Uniform(fDeltaPhiMinCut,fDeltaPhiMaxCut);
4069
4070       xE = -(pt/ptTrig)*TMath::Cos(randomphi);
4071       if(xE > 0 ) hbpXE = TMath::Log(1./xE);
4072
4073       fhDeltaPhiUeNeutralPt->Fill(pt,deltaPhi);
4074       fhZTUeNeutral        ->Fill(ptTrig,zT);
4075       fhPtHbpZTUeNeutral   ->Fill(ptTrig,hbpZT);
4076       fhXEUeNeutral        ->Fill(ptTrig,xE);
4077       fhPtHbpXEUeNeutral   ->Fill(ptTrig,hbpXE);
4078     }
4079
4080     // Several UE calculation, not sure it is useful
4081     // with partical calorimter acceptance
4082     if(fMakeSeveralUE) FillNeutralUnderlyingEventSidesHistograms(ptTrig,pt,zT,hbpZT,deltaPhi);
4083
4084     //
4085     // Decay photon correlations
4086     //
4087     if(fPi0Trigger && decayFound)
4088       FillDecayPhotonCorrelationHistograms(pt, phi, kFALSE) ;
4089
4090     if(fFillAODWithReferences)
4091     {
4092       nrefs++;
4093       if(nrefs==1)
4094       {
4095         refpi0 = new TObjArray(0);
4096         refpi0->SetName(GetAODObjArrayName()+"Pi0s");
4097         refpi0->SetOwner(kFALSE);
4098       }
4099       refpi0->Add(pi0);
4100     } // put references in trigger AOD
4101
4102     AliDebug(1,Form("Selected pi0: pt %2.2f, phi %2.2f, eta %2.2f",pt,phi,eta));
4103
4104   }//loop
4105
4106   //Fill AOD with reference tracks, if not filling histograms
4107   if(fFillAODWithReferences && refpi0)
4108   {
4109     aodParticle->AddObjArray(refpi0);
4110   }
4111 }
4112
4113 //__________________________________________________________________________________________________________________
4114 void  AliAnaParticleHadronCorrelation::MakeMCChargedCorrelation(Int_t label, Int_t histoIndex, Bool_t lostDecayPair)
4115 {
4116   // Charged Hadron Correlation Analysis with MC information
4117
4118   AliDebug(1,"Make trigger particle - charged hadron correlation in AOD MC level");
4119
4120   if( label < 0 )
4121   {
4122     AliDebug(1,Form(" *** bad label ***:  label %d", label));
4123     return;
4124   }
4125
4126   // Do MC correlation for a given particle type range.
4127   // Types defined in GetMCTagHistogramIndex:
4128   // 0 direct gamma; 1 pi0; 2 pi0 decay; 3 eta decay; 4 other decay; 5 electron; 6 other (hadron)
4129   if(histoIndex < fMCGenTypeMin || histoIndex > fMCGenTypeMax) return ;
4130
4131   AliStack         * stack        = 0x0 ;
4132   TParticle        * primary      = 0x0 ;
4133   TClonesArray     * mcparticles  = 0x0 ;
4134   AliAODMCParticle * aodprimary   = 0x0 ;
4135
4136   Double_t eprim   = 0 ;
4137   Double_t ptprim  = 0 ;
4138   Double_t phiprim = 0 ;
4139   Double_t etaprim = 0 ;
4140   Int_t    nTracks = 0 ;
4141   Int_t iParticle  = 0 ;
4142
4143   Bool_t leadTrig = kTRUE;
4144
4145   if( GetReader()->ReadStack() )
4146   {
4147     stack =  GetMCStack() ;
4148     if(!stack)
4149     {
4150       AliFatal("Stack not available, is the MC handler called? STOP");
4151       return;
4152     }
4153
4154     //nTracks = stack->GetNtrack() ;
4155     nTracks = stack->GetNprimary();
4156     if( label >=  stack->GetNtrack() )
4157     {
4158       if(GetDebug() > 2)
4159         AliInfo(Form("*** large label ***:  label %d, n tracks %d", label, stack->GetNtrack()));
4160       return ;
4161     }
4162
4163     primary = stack->Particle(label);
4164     if ( !primary )
4165     {
4166       AliInfo(Form(" *** no primary ***:  label %d", label));
4167       return;
4168     }
4169
4170     eprim    = primary->Energy();
4171     ptprim   = primary->Pt();
4172     etaprim  = primary->Eta();
4173     phiprim  = primary->Phi();
4174     if(phiprim < 0) phiprim+=TMath::TwoPi();
4175
4176     if(ptprim < 0.01 || eprim < 0.01) return ;
4177
4178     for (iParticle = 0 ; iParticle <  nTracks ; iParticle++)
4179     {
4180       TParticle * particle = stack->Particle(iParticle);
4181
4182       //keep only final state particles
4183       if( particle->GetStatusCode() != 1 ) continue ;
4184
4185       //---------- Charged particles ----------------------
4186       Int_t pdg    = particle->GetPdgCode();
4187       Int_t charge = (Int_t) TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
4188       if(charge == 0) continue;
4189
4190       particle->Momentum(fMomentum);
4191
4192       //Particles in CTS acceptance, make sure to use the same selection as in the reader
4193       Bool_t inCTS =  GetReader()->GetFiducialCut()->IsInFiducialCut(fMomentum.Eta(),fMomentum.Phi(),kCTS);
4194       //printf("Accepted? %d; pt %2.2f, eta %2.2f, phi %2.2f\n",inCTS,fMomentum.Pt(),fMomentum.Eta(),fMomentum.Phi()*TMath::RadToDeg());
4195       if( !inCTS ) continue;
4196
4197       // Remove conversions
4198       if ( TMath::Abs(pdg) == 11 && stack->Particle(particle->GetFirstMother())->GetPdgCode() == 22 ) continue ;
4199
4200       if ( label == iParticle ) continue; // avoid trigger particle
4201
4202       Float_t phi = particle->Phi();
4203       if(phi < 0) phi+=TMath::TwoPi();
4204
4205       Bool_t lead = FillChargedMCCorrelationHistograms(particle->Pt(),phi,particle->Eta(),ptprim,phiprim,etaprim,histoIndex,lostDecayPair);
4206       if(!lead) leadTrig = kFALSE;
4207       //if ( !lead && (fMakeAbsoluteLeading || fMakeNearSideLeading) ) return;
4208
4209     } //track loop
4210
4211   } //ESD MC
4212
4213   else if( GetReader()->ReadAODMCParticles() )
4214   {
4215     //Get the list of MC particles
4216     mcparticles = GetReader()->GetAODMCParticles();
4217     if( !mcparticles ) return;
4218
4219     nTracks = mcparticles->GetEntriesFast() ;
4220
4221     if( label >= nTracks )
4222     {
4223       if(GetDebug() > 2)
4224         AliInfo(Form(" *** large label ***:  label %d, n tracks %d", label,nTracks));
4225       return;
4226     }
4227
4228     //Get the particle
4229     aodprimary = (AliAODMCParticle*) mcparticles->At(label);
4230     if( !aodprimary )
4231     {
4232       AliInfo(Form(" *** no AOD primary ***:  label %d", label));
4233       return;
4234     }
4235
4236     eprim   = aodprimary->E();
4237     ptprim  = aodprimary->Pt();
4238     etaprim = aodprimary->Eta();
4239     phiprim = aodprimary->Phi();
4240     if(phiprim < 0) phiprim+=TMath::TwoPi();
4241
4242     if(ptprim < 0.01 || eprim < 0.01) return ;
4243
4244     for (iParticle = 0; iParticle < nTracks; iParticle++)
4245     {
4246       AliAODMCParticle *part = (AliAODMCParticle*) mcparticles->At(iParticle);
4247
4248       if (!part->IsPhysicalPrimary() ) continue; // same as part->GetStatus() !=1
4249
4250       if ( part->Charge() == 0 ) continue;
4251
4252       fMomentum.SetPxPyPzE(part->Px(),part->Py(),part->Pz(),part->E());
4253
4254       //Particles in CTS acceptance, make sure to use the same selection as in the reader
4255       Bool_t inCTS =  GetReader()->GetFiducialCut()->IsInFiducialCut(fMomentum.Eta(),fMomentum.Phi(),kCTS);
4256       //printf("Accepted? %d; pt %2.2f, eta %2.2f, phi %2.2f\n",inCTS,fMomentum.Pt(),fMomentum.Eta(),fMomentum.Phi()*TMath::RadToDeg());
4257       if( !inCTS ) continue;
4258
4259       // Remove conversions
4260       Int_t indexmother = part->GetMother();
4261       if ( indexmother > -1 )
4262       {
4263         Int_t pdg = part->GetPdgCode();
4264         Int_t mPdg = ((AliAODMCParticle*) mcparticles->At(indexmother)) ->GetPdgCode();
4265         if (TMath::Abs(pdg) == 11 && mPdg == 22) continue;
4266       }
4267
4268       if ( label == iParticle ) continue; // avoid trigger particle
4269
4270       Float_t phi = part->Phi();
4271       if(phi < 0) phi+=TMath::TwoPi();
4272
4273       Bool_t lead = FillChargedMCCorrelationHistograms(part->Pt(),phi,part->Eta(),ptprim,phiprim,etaprim, histoIndex,lostDecayPair);
4274       if(!lead) leadTrig = kFALSE;
4275       //if ( !lead && (fMakeAbsoluteLeading || fMakeNearSideLeading)) return;
4276
4277     }  //MC particle loop
4278   }// AOD MC
4279
4280   // Trigger MC particle histograms
4281   //if (!lead  && (fMakeAbsoluteLeading || fMakeNearSideLeading)) return;
4282
4283   fhMCPtTrigger [histoIndex]->Fill(ptprim);
4284   fhMCPhiTrigger[histoIndex]->Fill(ptprim,phiprim);
4285   fhMCEtaTrigger[histoIndex]->Fill(ptprim,etaprim);
4286
4287   if(lostDecayPair)
4288   {
4289     // check index in GetMCTagIndexHistogram
4290     if     (histoIndex == 2 && 8 >= fMCGenTypeMin && 8 <= fMCGenTypeMax )
4291     {
4292       // pi0 decay
4293       fhMCPtTrigger [8]->Fill(ptprim);
4294       fhMCPhiTrigger[8]->Fill(ptprim,phiprim);
4295       fhMCEtaTrigger[8]->Fill(ptprim,etaprim);
4296     }
4297     else if(histoIndex == 4 && 9 >= fMCGenTypeMin && 9 <= fMCGenTypeMax )
4298     {
4299       // eta decay
4300       fhMCPtTrigger [9]->Fill(ptprim);
4301       fhMCPhiTrigger[9]->Fill(ptprim,phiprim);
4302       fhMCEtaTrigger[9]->Fill(ptprim,etaprim);
4303     }
4304   }
4305
4306   if(!leadTrig && (fMakeAbsoluteLeading || fMakeNearSideLeading) )
4307   {
4308     AliDebug(1,Form("Not leading primary trigger: pT %2.2f, phi %2.2f, eta %2.2f",
4309                     ptprim,phiprim*TMath::RadToDeg(),etaprim));
4310
4311     fhMCPtTriggerNotLeading [histoIndex]->Fill(ptprim);
4312     fhMCPhiTriggerNotLeading[histoIndex]->Fill(ptprim,phiprim);
4313     fhMCEtaTriggerNotLeading[histoIndex]->Fill(ptprim,etaprim);
4314
4315     if(lostDecayPair)
4316     {
4317       // check index in GetMCTagIndexHistogram
4318       if      (histoIndex == 2  && 8 >= fMCGenTypeMin && 8 <= fMCGenTypeMax )
4319       {
4320         // pi0 decay
4321         fhMCPtTriggerNotLeading [8]->Fill(ptprim);
4322         fhMCPhiTriggerNotLeading[8]->Fill(ptprim,phiprim);
4323         fhMCEtaTriggerNotLeading[8]->Fill(ptprim,etaprim);
4324       }
4325       else  if(histoIndex == 4  && 9 >= fMCGenTypeMin && 9 <= fMCGenTypeMax )
4326       {
4327         // eta decay
4328         fhMCPtTriggerNotLeading [9]->Fill(ptprim);
4329         fhMCPhiTriggerNotLeading[9]->Fill(ptprim,phiprim);
4330         fhMCEtaTriggerNotLeading[9]->Fill(ptprim,etaprim);
4331       }
4332     }
4333   }
4334 }
4335
4336 //_____________________________________________________________________
4337 void AliAnaParticleHadronCorrelation::Print(const Option_t * opt) const
4338 {
4339
4340   //Print some relevant parameters set for the analysis
4341   if(! opt)
4342     return;
4343
4344   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
4345   AliAnaCaloTrackCorrBaseClass::Print(" ");
4346   printf("Pt trigger > %2.2f; < %2.2f\n", GetMinPt() , GetMaxPt()) ;
4347   printf("Pt associa > %2.2f; < %2.2f\n", fMinAssocPt, fMaxAssocPt) ;
4348   printf("Phi trigger particle-Hadron   <  %3.2f\n", fDeltaPhiMaxCut) ;
4349   printf("Phi trigger particle-Hadron   >  %3.2f\n", fDeltaPhiMinCut) ;
4350   printf("Phi trigger particle-UeHadron <  %3.2f\n", fUeDeltaPhiMaxCut) ;
4351   printf("Phi trigger particle-UeHadron >  %3.2f\n", fUeDeltaPhiMinCut) ;
4352   printf("Isolated Trigger?  %d\n"     , fSelectIsolated) ;
4353   printf("Several UE?  %d\n"           , fMakeSeveralUE) ;
4354   printf("Name of AOD Pi0 Branch %s \n", fPi0AODBranchName.Data());
4355   printf("Do Decay-hadron correlation ?  %d\n", fPi0Trigger) ;
4356   printf("Select absolute leading for cluster triggers ?  %d or Near Side %d\n",
4357          fMakeAbsoluteLeading, fMakeNearSideLeading) ;
4358   printf("Trigger pt bins  %d\n", fNAssocPtBins) ;
4359   for (Int_t ibin = 0; ibin<fNAssocPtBins; ibin++) {
4360     printf("\t bin %d = [%2.1f,%2.1f]\n", ibin, fAssocPtBinLimit[ibin], fAssocPtBinLimit[ibin+1]) ;
4361   }
4362
4363 }
4364
4365 //____________________________________________________________
4366 void AliAnaParticleHadronCorrelation::SetNAssocPtBins(Int_t n)
4367 {
4368   // Set number of bins
4369
4370   fNAssocPtBins  = n ;
4371
4372   if(n < 20 && n > 0)
4373   {
4374     fNAssocPtBins  = n ;
4375   }
4376   else
4377   {
4378     AliWarning("n = larger than 19 or too small, set to 19");
4379     fNAssocPtBins = 19;
4380   }
4381 }
4382
4383 //______________________________________________________________________________
4384 void AliAnaParticleHadronCorrelation::SetAssocPtBinLimit(Int_t ibin, Float_t pt)
4385 {
4386   // Set the list of limits for the trigger pt bins
4387
4388   if(ibin <= fNAssocPtBins || ibin >= 0)
4389   {
4390     fAssocPtBinLimit[ibin] = pt  ;
4391   }
4392   else
4393   {
4394     AliWarning(Form("Bin  number too large %d > %d or small, nothing done", ibin, fNAssocPtBins)) ;
4395   }
4396 }
4397