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) ||
2955            GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCISR)          ) return 0;
2956   else if( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0)          ) return 1;
2957   else if( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0Decay)     ) return 2;
2958   else if( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEta)          ) return 3;
2959   else if( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEtaDecay)     ) return 4;
2960   else if( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton)       ) return 5; // other decays
2961   else if(!GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCElectron)     ) return 6;
2962   else                                                                                    return 7;
2963
2964 }
2965
2966 //_________________________________________
2967 void AliAnaParticleHadronCorrelation::Init()
2968 {
2969   //Init
2970   //Do some checks
2971
2972   if(!GetReader()->IsCTSSwitchedOn())
2973     AliFatal("STOP!: You want to use CTS tracks in analysis but not read!! \n!!Check the configuration file!!");
2974 }
2975
2976 //____________________________________________________
2977 void AliAnaParticleHadronCorrelation::InitParameters()
2978 {
2979
2980   //Initialize the parameters of the analysis.
2981   SetInputAODName("Particle");
2982   SetAODObjArrayName("Hadrons");
2983   AddToHistogramsName("AnaHadronCorr_");
2984
2985   SetPtCutRange(0.,300);
2986   fDeltaPhiMinCut       = TMath::DegToRad()*120.;
2987   fDeltaPhiMaxCut       = TMath::DegToRad()*240. ;
2988   fSelectIsolated       = kFALSE;
2989   fMakeSeveralUE        = kFALSE;
2990   fUeDeltaPhiMinCut     = TMath::DegToRad()*60.;
2991   fUeDeltaPhiMaxCut     = TMath::DegToRad()*120 ;
2992
2993   fNeutralCorr          = kFALSE ;
2994   fPi0Trigger           = kFALSE ;
2995   fDecayTrigger         = kFALSE ;
2996   fHMPIDCorrelation     = kFALSE ;
2997
2998   fMakeAbsoluteLeading  = kTRUE;
2999   fMakeNearSideLeading  = kFALSE;
3000
3001   fNAssocPtBins         = 9   ;
3002   fAssocPtBinLimit[0]   = 0.2 ;
3003   fAssocPtBinLimit[1]   = 0.5 ;
3004   fAssocPtBinLimit[2]   = 1.0 ;
3005   fAssocPtBinLimit[3]   = 2.0 ;
3006   fAssocPtBinLimit[4]   = 3.0 ;
3007   fAssocPtBinLimit[5]   = 4.0 ;
3008   fAssocPtBinLimit[6]   = 5.0 ;
3009   fAssocPtBinLimit[7]   = 6.0 ;
3010   fAssocPtBinLimit[8]   = 7.0 ;
3011   fAssocPtBinLimit[9]   = 8.0 ;
3012   fAssocPtBinLimit[10]  = 9.0 ;
3013   fAssocPtBinLimit[11]  = 10.0 ;
3014   fAssocPtBinLimit[12]  = 12.0 ;
3015   fAssocPtBinLimit[13]  = 14.0 ;
3016   fAssocPtBinLimit[14]  = 16.0 ;
3017   fAssocPtBinLimit[15]  = 20.0 ;
3018   fAssocPtBinLimit[16]  = 30.0 ;
3019   fAssocPtBinLimit[17]  = 40.0 ;
3020   fAssocPtBinLimit[18]  = 50.0 ;
3021   fAssocPtBinLimit[19]  = 200.0 ;
3022
3023   fUseMixStoredInReader = kTRUE;
3024
3025   fM02MinCut   = -1 ;
3026   fM02MaxCut   = -1 ;
3027
3028   fSelectLeadingHadronAngle = kFALSE;
3029   fFillLeadHadOppositeHisto = kFALSE;
3030   fMinLeadHadPhi = 150*TMath::DegToRad();
3031   fMaxLeadHadPhi = 210*TMath::DegToRad();
3032
3033   fMinLeadHadPt  = 1;
3034   fMaxLeadHadPt  = 100;
3035
3036   fMCGenTypeMin =  0;
3037   fMCGenTypeMax = 10;
3038
3039   fNDecayBits = 1;
3040   fDecayBits[0] = AliNeutralMesonSelection::kPi0;
3041   fDecayBits[1] = AliNeutralMesonSelection::kEta;
3042   fDecayBits[2] = AliNeutralMesonSelection::kPi0Side;
3043   fDecayBits[3] = AliNeutralMesonSelection::kEtaSide;
3044 }
3045
3046 //_________________________________________________________________________
3047 Bool_t AliAnaParticleHadronCorrelation::IsTriggerTheEventLeadingParticle()
3048 {
3049   // Check if the what of the selected triggers is leading particle comparing
3050   // with all the triggers, all the tracks or all the clusters (if requested for the clusters).
3051
3052   Double_t ptTrig      = GetMinPt();
3053   Double_t phiTrig     = 0 ;
3054   fLeadingTriggerIndex =-1 ;
3055   Int_t index          =-1 ;
3056   AliAODPWG4ParticleCorrelation* pLeading = 0;
3057
3058   // Loop on stored AOD particles, find leading trigger on the selected list, with at least min pT.
3059
3060   for(Int_t iaod = 0; iaod < GetInputAODBranch()->GetEntriesFast() ; iaod++)
3061   {
3062     AliAODPWG4ParticleCorrelation* particle =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
3063     particle->SetLeadingParticle(kFALSE); // set it later
3064
3065     // Vertex cut in case of mixing
3066     Int_t check = CheckMixedEventVertex(particle->GetCaloLabel(0), particle->GetTrackLabel(0));
3067     if(check ==  0) continue;
3068     if(check == -1) return kFALSE; // not sure if it is correct.
3069
3070     // find the leading particles with highest momentum
3071     if (particle->Pt() > ptTrig)
3072     {
3073       ptTrig   = particle->Pt() ;
3074       phiTrig  = particle->Phi();
3075       index    = iaod     ;
3076       pLeading = particle ;
3077     }
3078   }// finish search of leading trigger particle on the AOD branch.
3079
3080   if(index < 0) return kFALSE;
3081
3082   //printf("AOD leading pT %2.2f, ID %d\n", pLeading->Pt(),pLeading->GetCaloLabel(0));
3083
3084   if(phiTrig < 0 ) phiTrig += TMath::TwoPi();
3085
3086   // Compare if it is the leading of all tracks
3087
3088   for(Int_t ipr = 0;ipr < GetCTSTracks()->GetEntriesFast() ; ipr ++ )
3089   {
3090     AliVTrack * track = (AliVTrack *) (GetCTSTracks()->At(ipr)) ;
3091
3092     if(track->GetID() == pLeading->GetTrackLabel(0) || track->GetID() == pLeading->GetTrackLabel(1) ||
3093        track->GetID() == pLeading->GetTrackLabel(2) || track->GetID() == pLeading->GetTrackLabel(3)   ) continue ;
3094
3095     fTrackVector.SetXYZ(track->Px(),track->Py(),track->Pz());
3096     Float_t pt   = fTrackVector.Pt();
3097     Float_t phi  = fTrackVector.Phi() ;
3098     if(phi < 0) phi+=TMath::TwoPi();
3099
3100     //jump out this event if near side associated particle pt larger than trigger
3101     if (fMakeNearSideLeading)
3102     {
3103       Float_t deltaPhi = phiTrig-phi;
3104       if(deltaPhi <= -TMath::PiOver2()) deltaPhi+=TMath::TwoPi();
3105       if(deltaPhi > 3*TMath::PiOver2()) deltaPhi-=TMath::TwoPi();
3106
3107       if(pt > ptTrig && deltaPhi < TMath::PiOver2())  return kFALSE;
3108     }
3109     //jump out this event if there is any other particle with pt larger than trigger
3110     else
3111     {
3112       if(pt > ptTrig)  return kFALSE ;
3113     }
3114    }// track loop
3115
3116   // Compare if it is leading of all calorimeter clusters
3117
3118   if(fCheckLeadingWithNeutralClusters)
3119   {
3120     // Select the calorimeter cluster list
3121     TObjArray * nePl = 0x0;
3122     if      (pLeading->GetDetectorTag() == kPHOS )
3123       nePl = GetPHOSClusters();
3124     else
3125       nePl = GetEMCALClusters();
3126
3127     if(!nePl) return kTRUE; // Do the selection just with the tracks if no calorimeter is available.
3128
3129     for(Int_t ipr = 0;ipr < nePl->GetEntriesFast() ; ipr ++ )
3130     {
3131       AliVCluster * cluster = (AliVCluster *) (nePl->At(ipr)) ;
3132
3133       if(cluster->GetID() == pLeading->GetCaloLabel(0) || cluster->GetID() == pLeading->GetCaloLabel(1) ) continue ;
3134
3135       cluster->GetMomentum(fMomentum,GetVertex(0));
3136
3137       Float_t pt   = fMomentum.Pt();
3138       Float_t phi  = fMomentum.Phi() ;
3139       if(phi < 0) phi+=TMath::TwoPi();
3140
3141       if(IsTrackMatched(cluster,GetReader()->GetInputEvent())) continue ; // avoid charged clusters, already covered by tracks, or cluster merging with track.
3142
3143       //jump out this event if near side associated particle pt larger than trigger
3144       // not really needed for calorimeter, unless DCal is included
3145       if (fMakeNearSideLeading)
3146       {
3147         Float_t deltaPhi = phiTrig-phi;
3148         if(deltaPhi <= -TMath::PiOver2()) deltaPhi+=TMath::TwoPi();
3149         if(deltaPhi > 3*TMath::PiOver2()) deltaPhi-=TMath::TwoPi();
3150
3151         if(pt > ptTrig && deltaPhi < TMath::PiOver2()) return kFALSE ;
3152       }
3153       //jump out this event if there is any other particle with pt larger than trigger
3154       else
3155       {
3156         if(pt > ptTrig)  return kFALSE ;
3157       }
3158     }// cluster loop
3159   } // check neutral clusters
3160
3161   fLeadingTriggerIndex = index ;
3162   pLeading->SetLeadingParticle(kTRUE);
3163
3164   AliDebug(1,Form("\t particle AOD with index %d is leading with pT %2.2f", fLeadingTriggerIndex, pLeading->Pt()));
3165
3166   return kTRUE;
3167
3168 }
3169
3170 //_________________________________________________________________
3171 void  AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms()
3172 {
3173   //Particle-Hadron Correlation Analysis, fill histograms
3174
3175   if(!GetInputAODBranch())
3176   {
3177     AliFatal(Form("No input particles in AOD with name branch < %s >, STOP",GetInputAODName().Data()));
3178     return ; // coverity
3179   }
3180
3181   Int_t naod = GetInputAODBranch()->GetEntriesFast();
3182   if( naod == 0 )
3183   {
3184     AliWarning("No particle AOD found!");
3185     return ; // no trigger particles found.
3186   }
3187
3188   AliDebug(1,Form("Begin hadron correlation analysis, fill histograms"));
3189   AliDebug(1,Form("n particle branch aod entries %d", naod));
3190   AliDebug(1,Form("In CTS aod entries %d",GetCTSTracks()->GetEntriesFast()));
3191
3192   //------------------------------------------------------
3193   // Find leading trigger if analysis request only leading,
3194   // if there is no leading trigger, then skip the event
3195
3196   Int_t iaod = 0 ;
3197   if( fMakeAbsoluteLeading || fMakeNearSideLeading )
3198   {
3199     Bool_t leading = IsTriggerTheEventLeadingParticle();
3200
3201     AliDebug(1,Form("AOD Leading trigger? %d, with index %d",leading,fLeadingTriggerIndex));
3202
3203     if(!leading)
3204     {
3205       AliDebug(1,"Leading was requested and not found");
3206       return ;
3207     }
3208     else
3209     {
3210       // Select only the leading in the trigger AOD loop
3211       naod = fLeadingTriggerIndex+1 ;
3212       iaod = fLeadingTriggerIndex   ;
3213     }
3214   }
3215
3216   //------------------------------------------------------
3217   // Get event multiplicity and bins
3218
3219   Float_t cen         = GetEventCentrality();
3220   Float_t ep          = GetEventPlaneAngle();
3221   if(IsHighMultiplicityAnalysisOn()) fhTriggerEventPlaneCentrality->Fill(cen,ep);
3222
3223   Int_t   mixEventBin = GetEventMixBin();
3224   Int_t   vzbin       = GetEventVzBin();
3225
3226   //------------------------------------------------------
3227   // Loop on trigger AOD
3228
3229   for( iaod = 0; iaod < naod; iaod++ )
3230   {
3231     AliAODPWG4ParticleCorrelation* particle =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
3232
3233     //
3234     // Trigger particle selection criteria:
3235     //
3236     Float_t pt = particle->Pt();
3237
3238     if(pt < GetMinPt() || pt > GetMaxPt() ) continue ;
3239
3240     fhPtTriggerInput->Fill(pt);
3241
3242     //
3243     // check if it was a calorimeter cluster
3244     // and if the shower shape cut was requested apply it.
3245     // Not needed if already done at the particle identification level,
3246     // but for isolation studies, it is preferred not to remove so we do it here
3247     //
3248
3249     AliDebug(1,Form("%s Trigger : min %f, max %f, det %d",
3250                     GetInputAODName().Data(),fM02MinCut,fM02MaxCut,particle->GetDetectorTag()));
3251
3252     if(fM02MaxCut > 0 && fM02MinCut > 0) //clID1 > 0 && clID2 < 0 &&
3253     {
3254 //      Int_t iclus = -1;
3255 //      TObjArray* clusters = 0x0;
3256 //      if     (particle->GetDetectorTag() == kEMCAL) clusters = GetEMCALClusters();
3257 //      else if(particle->GetDetectorTag() == kPHOS ) clusters = GetPHOSClusters();
3258 //
3259 //      if(clusters)
3260 //      {
3261 //        AliVCluster *cluster = FindCluster(clusters,clID1,iclus);
3262 //        Float_t m02 = cluster->GetM02();
3263 //        if(m02 > fM02MaxCut || m02 < fM02MinCut) continue ;
3264 //      }
3265
3266       Float_t m02 = particle->GetM02();
3267
3268       if(m02 > fM02MaxCut || m02 < fM02MinCut) continue ;
3269
3270       fhPtTriggerSSCut->Fill(pt);
3271
3272       AliDebug(1,"Pass the shower shape cut");
3273     }
3274
3275     //
3276     // Check if the particle is isolated or if we want to take the isolation into account
3277     // This bool is set in AliAnaParticleIsolation
3278     //
3279     if(OnlyIsolated())
3280     {
3281       if( !particle->IsIsolated() ) continue;
3282
3283       fhPtTriggerIsoCut->Fill(pt);
3284
3285       AliDebug(1,"Pass the isolation cut");
3286     }
3287
3288     //
3289     // Check if trigger is in fiducial region
3290     //
3291     if(IsFiducialCutOn())
3292     {
3293       Bool_t in = GetFiducialCut()->IsInFiducialCut(particle->Eta(),particle->Phi(),particle->GetDetectorTag()) ;
3294
3295       if(! in ) continue ;
3296
3297       AliDebug(1,"Pass the fiducial cut");
3298     }
3299
3300     fhPtTriggerFidCut->Fill(pt);
3301
3302     //---------------------------------------
3303     // Make correlation
3304
3305     // Find the leading hadron in the opposite hemisphere to the triggeer
3306     // and accept the trigger if leading is in defined window.
3307     Bool_t okLeadHad = kTRUE;
3308     if(fSelectLeadingHadronAngle || fFillLeadHadOppositeHisto)
3309     {
3310       okLeadHad = FindLeadingOppositeHadronInWindow(particle);
3311       if(!okLeadHad && fSelectLeadingHadronAngle) continue;
3312     }
3313
3314     //
3315     // Charged particles correlation
3316     //
3317     MakeChargedCorrelation(particle);
3318
3319     // MC
3320     Int_t mcIndex = -1;
3321     Int_t mcTag   = particle->GetTag();
3322     Bool_t lostDecayPair = kFALSE;
3323     if(IsDataMC())
3324     {
3325       mcIndex = GetMCTagHistogramIndex(mcTag);
3326       lostDecayPair = GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCDecayPairLost);
3327       MakeMCChargedCorrelation(particle->GetLabel(), mcIndex,lostDecayPair);
3328     }
3329
3330     // Do own mixed event with charged,
3331     // add event and remove previous or fill the mixed histograms
3332     if(DoOwnMix())
3333       MakeChargedMixCorrelation(particle);
3334
3335     //
3336     // Neutral particles correlation
3337     //
3338     if(fNeutralCorr)
3339       MakeNeutralCorrelation(particle);
3340
3341     //----------------------------------------------------------------
3342     // Fill trigger pT related histograms if not absolute leading
3343
3344     //
3345     // pT of the trigger, vs trigger origin if MC
3346     //
3347     fhPtTrigger->Fill(pt);
3348     if(IsDataMC() && mcIndex >=0 && mcIndex < fgkNmcTypes)
3349     {
3350       fhPtTriggerMC[mcIndex]->Fill(pt);
3351       if( lostDecayPair )
3352       {
3353         // check index in GetMCTagIndexHistogram
3354         if     ( mcIndex == 2 ) fhPtTriggerMC[8]->Fill(pt); // pi0 decay
3355         else if( mcIndex == 4 ) fhPtTriggerMC[9]->Fill(pt); // eta decay
3356       }
3357     }
3358
3359     if(fDecayTrigger)
3360     {
3361       Int_t decayTag = particle->DecayTag();
3362       if(decayTag < 0) decayTag = 0;
3363
3364       for(Int_t ibit = 0; ibit<fNDecayBits; ibit++)
3365       {
3366         if(GetNeutralMesonSelection()->CheckDecayBit(decayTag,fDecayBits[ibit]))
3367         {
3368           fhPtDecayTrigger[ibit]->Fill(pt);
3369
3370           if(IsDataMC() && mcIndex >=0 && mcIndex < fgkNmcTypes)
3371           {
3372             fhPtDecayTriggerMC[ibit][mcIndex]->Fill(pt);
3373             if( lostDecayPair )
3374             {
3375               // check index in GetMCTagIndexHistogram
3376               if( mcIndex == 2 )  fhPtDecayTriggerMC[ibit][8]->Fill(pt); // pi0 decay
3377               if( mcIndex == 4 )  fhPtDecayTriggerMC[ibit][9]->Fill(pt); // eta decay
3378             }
3379           }
3380         }// check bit
3381       }// bit loop
3382     }
3383
3384     //
3385     // Acceptance of the trigger
3386     //
3387     Float_t phi = particle->Phi();
3388     if( phi < 0 ) phi+=TMath::TwoPi();
3389     fhPhiTrigger->Fill(pt, phi);
3390
3391     fhEtaTrigger->Fill(pt, particle->Eta());
3392     //printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms() - Trigger particle : pt %f, eta %f, phi %f\n",particle->Pt(),particle->Eta(),phi);
3393
3394     //----------------------------------
3395     // Trigger particle pT vs event bins
3396
3397     fhPtTriggerBin->Fill(pt,mixEventBin);
3398     if(fCorrelVzBin)
3399       fhPtTriggerVzBin->Fill(pt,vzbin);
3400
3401     if(IsHighMultiplicityAnalysisOn())
3402     {
3403       fhPtTriggerCentrality->Fill(pt,cen);
3404       fhPtTriggerEventPlane->Fill(pt,ep);
3405     }
3406
3407     //----------------------------------
3408     // Trigger particle pT vs pile-up
3409
3410     if(IsPileUpAnalysisOn())
3411     {
3412       Int_t vtxBC = GetReader()->GetVertexBC();
3413       if(vtxBC == 0 || vtxBC==AliVTrack::kTOFBCNA)     fhPtTriggerVtxBC0->Fill(pt);
3414
3415       if(GetReader()->IsPileUpFromSPD())               fhPtTriggerPileUp[0]->Fill(pt);
3416       if(GetReader()->IsPileUpFromEMCal())             fhPtTriggerPileUp[1]->Fill(pt);
3417       if(GetReader()->IsPileUpFromSPDOrEMCal())        fhPtTriggerPileUp[2]->Fill(pt);
3418       if(GetReader()->IsPileUpFromSPDAndEMCal())       fhPtTriggerPileUp[3]->Fill(pt);
3419       if(GetReader()->IsPileUpFromSPDAndNotEMCal())    fhPtTriggerPileUp[4]->Fill(pt);
3420       if(GetReader()->IsPileUpFromEMCalAndNotSPD())    fhPtTriggerPileUp[5]->Fill(pt);
3421       if(GetReader()->IsPileUpFromNotSPDAndNotEMCal()) fhPtTriggerPileUp[6]->Fill(pt);
3422     }
3423   } // AOD trigger loop
3424
3425   //Reinit for next event
3426   fLeadingTriggerIndex = -1;
3427
3428   AliDebug(1,"End fill histograms");
3429 }
3430
3431 //_______________________________________________________________________________________________________
3432 void  AliAnaParticleHadronCorrelation::MakeChargedCorrelation(AliAODPWG4ParticleCorrelation *aodParticle)
3433 {
3434   // Charged Hadron Correlation Analysis
3435   AliDebug(1,"Make trigger particle - charged hadron correlation");
3436
3437   Float_t phiTrig = aodParticle->Phi();
3438   Float_t etaTrig = aodParticle->Eta();
3439   Float_t ptTrig  = aodParticle->Pt();
3440   Int_t    mcTag  = aodParticle->GetTag();
3441   Double_t bz     = GetReader()->GetInputEvent()->GetMagneticField();
3442
3443
3444   Int_t   decayTag = 0;
3445   if(fDecayTrigger)
3446   {
3447     //decay = aodParticle->IsTagged();
3448     decayTag = aodParticle->DecayTag();
3449     if(decayTag < 0) decayTag = 0;
3450 //    printf("Correlation: pT %2.2f, BTag %d, Tagged %d\n",ptTrig, decayTag, aodParticle->IsTagged());
3451 //    printf("\t check bit Pi0 %d, Eta %d,  Pi0Side %d, EtaSide %d\n",
3452 //           GetNeutralMesonSelection()->CheckDecayBit(decayTag,AliNeutralMesonSelection::kPi0),
3453 //           GetNeutralMesonSelection()->CheckDecayBit(decayTag,AliNeutralMesonSelection::kEta),
3454 //           GetNeutralMesonSelection()->CheckDecayBit(decayTag,AliNeutralMesonSelection::kPi0Side),
3455 //           GetNeutralMesonSelection()->CheckDecayBit(decayTag,AliNeutralMesonSelection::kEtaSide));
3456   }
3457
3458   Float_t pt       = -100. ;
3459   Float_t phi      = -100. ;
3460   Float_t eta      = -100. ;
3461   Float_t deltaPhi = -100. ;
3462
3463   TObjArray * reftracks = 0x0;
3464   Int_t nrefs           = 0;
3465
3466   // Mixed event settings
3467   Int_t evtIndex11   = -1 ; // cluster trigger or pi0 trigger
3468   Int_t evtIndex12   = -1 ; // pi0 trigger
3469   Int_t evtIndex13   = -1 ; // charged trigger
3470
3471   if (GetMixedEvent())
3472   {
3473     evtIndex11 = GetMixedEvent()->EventIndexForCaloCluster(aodParticle->GetCaloLabel(0)) ;
3474     evtIndex12 = GetMixedEvent()->EventIndexForCaloCluster(aodParticle->GetCaloLabel(1)) ;
3475     evtIndex13 = GetMixedEvent()->EventIndex(aodParticle->GetTrackLabel(0)) ;
3476   }
3477
3478   // Track multiplicity or cent bin
3479   Int_t cenbin = 0;
3480   if(IsHighMultiplicityAnalysisOn()) cenbin = GetEventCentralityBin();
3481
3482   //
3483   // In case of pi0/eta trigger, we may want to check their decay correlation,
3484   // get their decay children
3485   //
3486
3487   Bool_t decayFound = kFALSE;
3488   if( fPi0Trigger )
3489   {
3490     decayFound = GetDecayPhotonMomentum(aodParticle->GetCaloLabel(0),aodParticle->GetCaloLabel(1),aodParticle->GetDetectorTag());
3491     if(decayFound)
3492     {
3493       fhPtPi0DecayRatio->Fill(ptTrig, fDecayMom1.Pt()/ptTrig);
3494       fhPtPi0DecayRatio->Fill(ptTrig, fDecayMom2.Pt()/ptTrig);
3495     }
3496   }
3497
3498   //-----------------------------------------------------------------------
3499   // Track loop, select tracks with good pt, phi and fill AODs or histograms
3500   //-----------------------------------------------------------------------
3501
3502   for(Int_t ipr = 0;ipr < GetCTSTracks()->GetEntriesFast() ; ipr ++ )
3503   {
3504     AliVTrack * track = (AliVTrack *) (GetCTSTracks()->At(ipr)) ;
3505
3506     fTrackVector.SetXYZ(track->Px(),track->Py(),track->Pz());
3507     pt   = fTrackVector.Pt();
3508     eta  = fTrackVector.Eta();
3509     phi  = fTrackVector.Phi() ;
3510     if(phi < 0) phi+=TMath::TwoPi();
3511
3512     //Select only hadrons in pt range
3513     if(pt < fMinAssocPt || pt > fMaxAssocPt) continue ;
3514
3515     //remove trigger itself for correlation when use charged triggers
3516     if( track->GetID() == aodParticle->GetTrackLabel(0) || track->GetID() == aodParticle->GetTrackLabel(1) ||
3517         track->GetID() == aodParticle->GetTrackLabel(2) || track->GetID() == aodParticle->GetTrackLabel(3)   )
3518       continue ;
3519
3520     //Only for mixed event frame
3521     Int_t evtIndex2 = 0 ;
3522     if (GetMixedEvent())
3523     {
3524       evtIndex2 = GetMixedEvent()->EventIndex(track->GetID()) ;
3525       if (evtIndex11 == evtIndex2 || evtIndex12 == evtIndex2 || evtIndex13 == evtIndex2 ) // photon and track from different events
3526         continue ;
3527       //vertex cut
3528       if (TMath::Abs(GetVertex(evtIndex2)[2]) > GetZvertexCut())
3529         continue;
3530     }
3531
3532      AliDebug(2,Form("Selected charge for momentum imbalance: pt %2.2f, phi %2.2f, eta %2.2f",pt,phi,eta));
3533
3534     // ------------------------------
3535     // Track type bin or bits setting
3536     //
3537
3538     //
3539     // * Set the pt associated bin for the defined bins *
3540     //
3541     Int_t assocBin   = -1;
3542     for(Int_t i = 0 ; i < fNAssocPtBins ; i++)
3543     {
3544       if(pt > fAssocPtBinLimit[i] && pt < fAssocPtBinLimit[i+1]) assocBin= i;
3545     }
3546
3547     //
3548     // * Assign to the histogram array a bin corresponding
3549     // to a combination of pTa and vz bins *
3550     //
3551     Int_t nz = 1;
3552     Int_t vz = 0;
3553
3554     if(fCorrelVzBin)
3555     {
3556       nz = GetNZvertBin();
3557       vz = GetEventVzBin();
3558     }
3559
3560     Int_t bin = assocBin*nz+vz;
3561
3562     //printf("assoc Bin = %d, vZ bin  = %d, bin = %d \n", assocBin,GetEventVzBin(),bin);
3563
3564     //
3565     // * Get the status of the TOF bit *
3566     //
3567     ULong_t status = track->GetStatus();
3568     Bool_t okTOF = ( (status & AliVTrack::kTOFout) == AliVTrack::kTOFout ) ;
3569     //Double32_t tof = track->GetTOFsignal()*1e-3;
3570     Int_t trackBC = track->GetTOFBunchCrossing(bz);
3571
3572     Int_t outTOF = -1;
3573     if     (okTOF && trackBC!=0) outTOF = 1;
3574     else if(okTOF && trackBC==0) outTOF = 0;
3575
3576     //----------------
3577     // Fill Histograms
3578
3579     //
3580     // Azimuthal Angle histograms
3581     //
3582
3583     deltaPhi = phiTrig-phi;
3584
3585     //
3586     // Calculate deltaPhi shift so that for the particles on the opposite side
3587     // it is defined between 90 and 270 degrees
3588     // Shift [-360,-90]  to [0, 270]
3589     // and [270,360] to [-90,0]
3590     if(deltaPhi <= -TMath::PiOver2()) deltaPhi+=TMath::TwoPi();
3591     if(deltaPhi > 3*TMath::PiOver2()) deltaPhi-=TMath::TwoPi();
3592
3593     FillChargedAngularCorrelationHistograms(pt,  ptTrig,  bin, phi, phiTrig,  deltaPhi,
3594                                             eta, etaTrig, decayTag, track->GetHMPIDsignal(),
3595                                             outTOF, cenbin, mcTag);
3596
3597     //
3598     // Imbalance zT/xE/pOut histograms
3599     //
3600
3601     //
3602     // Delta phi cut for momentum imbalance correlation
3603     //
3604     if  ( (deltaPhi > fDeltaPhiMinCut)   && (deltaPhi < fDeltaPhiMaxCut)   )
3605       FillChargedMomentumImbalanceHistograms(ptTrig, pt, deltaPhi, cenbin, track->Charge(),
3606                                              assocBin, decayTag, outTOF, mcTag);
3607
3608     //
3609     // Underlying event, right side, default case
3610     //
3611     if ( (deltaPhi > fUeDeltaPhiMinCut) && (deltaPhi < fUeDeltaPhiMaxCut) )
3612       FillChargedUnderlyingEventHistograms(ptTrig, pt, deltaPhi, cenbin, outTOF);
3613
3614     //
3615     // Several UE calculation,  in different perpendicular regions, up to 6:
3616     // left, right, upper-left, lower left, upper-right, lower-right
3617     //
3618     if(fMakeSeveralUE)
3619       FillChargedUnderlyingEventSidesHistograms(ptTrig,pt,deltaPhi);
3620
3621     //
3622     if(fPi0Trigger && decayFound)
3623       FillDecayPhotonCorrelationHistograms(pt, phi, kTRUE) ;
3624
3625     //
3626     // Add track reference to array
3627     //
3628     if(fFillAODWithReferences)
3629     {
3630       nrefs++;
3631       if(nrefs==1)
3632       {
3633         reftracks = new TObjArray(0);
3634         TString trackname = Form("%sTracks", GetAODObjArrayName().Data());
3635         reftracks->SetName(trackname.Data());
3636         reftracks->SetOwner(kFALSE);
3637       }
3638
3639       reftracks->Add(track);
3640     }// reference track to AOD
3641   }// track loop
3642
3643   //Fill AOD with reference tracks, if not filling histograms
3644   if(fFillAODWithReferences && reftracks)
3645   {
3646     aodParticle->AddObjArray(reftracks);
3647   }
3648
3649 }
3650
3651 //_________________________________________________________________________________________________________
3652 void AliAnaParticleHadronCorrelation::MakeChargedMixCorrelation(AliAODPWG4ParticleCorrelation *aodParticle)
3653 {
3654   // Mix current trigger with tracks in another MB event
3655
3656   AliDebug(1,Form("Make trigger particle - charged hadron mixed event correlation"));
3657
3658   if(GetMixedEvent()) return;  // This is not the mixed event from general mixing frame
3659
3660   // Get the event with similar caracteristics
3661   //printf("MakeChargedMixCorrelation for %s\n",GetInputAODName().Data());
3662
3663   AliAnalysisManager   * manager      = AliAnalysisManager::GetAnalysisManager();
3664
3665   AliInputEventHandler * inputHandler = dynamic_cast<AliInputEventHandler*>(manager->GetInputEventHandler());
3666
3667   if(!inputHandler) return;
3668
3669   if(!(inputHandler->IsEventSelected( ) & GetReader()->GetEventTriggerMask())) return;
3670
3671   // Get the pool, check if it exits
3672   Int_t eventBin = GetEventMixBin();
3673
3674   //Check that the bin exists, if not (bad determination of RP, ntrality or vz bin) do nothing
3675   if(eventBin < 0) return;
3676
3677   fhEventBin->Fill(eventBin);
3678
3679   // get neutral clusters pool?
3680   Bool_t isoCase = OnlyIsolated() && (GetIsolationCut()->GetParticleTypeInCone() != AliIsolationCut::kOnlyCharged);
3681   Bool_t neutralMix = fFillNeutralEventMixPool || isoCase ;
3682
3683   TList * pool     = 0;
3684   TList * poolCalo = 0;
3685   if(fUseMixStoredInReader)
3686   {
3687     pool     = GetReader()->GetListWithMixedEventsForTracks(eventBin);
3688     if(neutralMix) poolCalo = GetReader()->GetListWithMixedEventsForCalo  (eventBin);
3689   }
3690   else
3691   {
3692     pool     = fListMixTrackEvents[eventBin];
3693     if(neutralMix) poolCalo = fListMixCaloEvents [eventBin];
3694   }
3695
3696   if(!pool) return ;
3697
3698   if( neutralMix && !poolCalo )
3699     AliWarning("Careful, cluster pool not available");
3700
3701   Double_t ptTrig  = aodParticle->Pt();
3702   Double_t etaTrig = aodParticle->Eta();
3703   Double_t phiTrig = aodParticle->Phi();
3704   if(phiTrig < 0.) phiTrig+=TMath::TwoPi();
3705
3706   AliDebug(1,Form("Pool bin %d size %d, trigger trigger pt=%f, phi=%f, eta=%f",
3707                   eventBin,pool->GetSize(), ptTrig,phiTrig,etaTrig));
3708
3709   Double_t ptAssoc  = -999.;
3710   Double_t phiAssoc = -999.;
3711   Double_t etaAssoc = -999.;
3712   Double_t deltaPhi = -999.;
3713   Double_t deltaEta = -999.;
3714   Double_t xE       = -999.;
3715
3716   // Start from first event in pool except if in this same event the pool was filled
3717   Int_t ev0 = 0;
3718   if(GetReader()->GetLastTracksMixedEvent() == GetEventNumber()) ev0 = 1;
3719
3720   for(Int_t ev=ev0; ev < pool->GetSize(); ev++)
3721   {
3722     //
3723     // Recover the lists of tracks or clusters
3724     //
3725     TObjArray* bgTracks = static_cast<TObjArray*>(pool->At(ev));
3726     TObjArray* bgCalo   = 0;
3727
3728     // Recover the clusters list if requested
3729     if( neutralMix && poolCalo )
3730     {
3731       if(pool->GetSize()!=poolCalo->GetSize())
3732         AliWarning("Different size of calo and track pools");
3733
3734       bgCalo = static_cast<TObjArray*>(poolCalo->At(ev));
3735
3736       if(!bgCalo) AliDebug(1,Form("Event %d in calo pool not available?",ev));
3737     }
3738
3739     //
3740     // Isolate the trigger in the mixed event with mixed tracks and clusters
3741     //
3742     if( OnlyIsolated() )
3743     {
3744       Int_t   n=0, nfrac = 0;
3745       Bool_t  isolated = kFALSE;
3746       Float_t coneptsum = 0, coneptlead = 0;
3747       GetIsolationCut()->MakeIsolationCut(bgTracks,bgCalo,
3748                                           GetReader(), GetCaloPID(),
3749                                           kFALSE, aodParticle, "",
3750                                           n,nfrac,coneptsum,coneptlead,isolated);
3751
3752       //printf("AliAnaParticleHadronCorrelation::MakeChargedMixCorrelation() - Isolated? %d - cone %f, ptthres %f",
3753       //       isolated,GetIsolationCut()->GetConeSize(),GetIsolationCut()->GetPtThreshold());
3754       //if(bgTracks)printf(" - n track %d", bgTracks->GetEntriesFast());
3755       //printf("\n");
3756
3757       if(!isolated) continue ;
3758     }
3759
3760     //
3761     // Check if the trigger is leading of mixed event
3762     //
3763     Int_t nTracks=bgTracks->GetEntriesFast();
3764
3765     if(fMakeNearSideLeading || fMakeAbsoluteLeading)
3766     {
3767       Bool_t leading = kTRUE;
3768       for(Int_t jlead = 0;jlead < nTracks; jlead++ )
3769       {
3770         AliAODPWG4Particle *track = (AliAODPWG4Particle*) bgTracks->At(jlead) ;
3771
3772         ptAssoc  = track->Pt();
3773         phiAssoc = track->Phi() ;
3774         if(phiAssoc < 0) phiAssoc+=TMath::TwoPi();
3775
3776         if (fMakeNearSideLeading)
3777         {
3778           deltaPhi = phiTrig-phiAssoc;
3779           if(deltaPhi <= -TMath::PiOver2()) deltaPhi+=TMath::TwoPi();
3780           if(deltaPhi > 3*TMath::PiOver2()) deltaPhi-=TMath::TwoPi();
3781
3782           if(ptAssoc > ptTrig && deltaPhi < TMath::PiOver2())
3783           {
3784             leading = kFALSE;
3785             break;
3786           }
3787         }
3788         //jump out this event if there is any other particle with pt larger than trigger
3789         else if(fMakeAbsoluteLeading)
3790         {
3791           if(ptAssoc > ptTrig)
3792           {
3793             leading = kFALSE;
3794             break;
3795           }
3796         }
3797       }
3798
3799       if( !neutralMix && fCheckLeadingWithNeutralClusters )
3800         AliWarning("Leading of clusters requested but no clusters in mixed event");
3801
3802       if(neutralMix && fCheckLeadingWithNeutralClusters && bgCalo)
3803       {
3804         Int_t nClusters=bgCalo->GetEntriesFast();
3805         for(Int_t jlead = 0;jlead <nClusters; jlead++ )
3806         {
3807           AliAODPWG4Particle *cluster= (AliAODPWG4Particle*) bgCalo->At(jlead) ;
3808
3809           ptAssoc  = cluster->Pt();
3810           phiAssoc = cluster->Phi() ;
3811           if(phiAssoc < 0) phiAssoc+=TMath::TwoPi();
3812
3813           if (fMakeNearSideLeading)
3814           {
3815             deltaPhi = phiTrig-phiAssoc;
3816             if(deltaPhi <= -TMath::PiOver2()) deltaPhi+=TMath::TwoPi();
3817             if(deltaPhi > 3*TMath::PiOver2()) deltaPhi-=TMath::TwoPi();
3818
3819             if(ptAssoc > ptTrig && deltaPhi < TMath::PiOver2())
3820             {
3821               leading = kFALSE;
3822               break;
3823             }
3824           }
3825           //jump out this event if there is any other particle with pt larger than trigger
3826           else if(fMakeAbsoluteLeading)
3827           {
3828             if(ptAssoc > ptTrig)
3829             {
3830               leading = kFALSE;
3831               break;
3832             }
3833           }
3834         }
3835       }
3836
3837       if(!leading) continue; // not leading, check the next event in pool
3838     }
3839
3840     //
3841     // Fill histograms for selected triggers
3842     //
3843
3844     fhEventMixBin->Fill(eventBin);
3845
3846     //printf("\t Read Pool event %d, nTracks %d\n",ev,nTracks);
3847
3848     fhPtTriggerMixed   ->Fill(ptTrig);
3849     fhPhiTriggerMixed  ->Fill(ptTrig, phiTrig);
3850     fhEtaTriggerMixed  ->Fill(ptTrig, etaTrig);
3851     fhPtTriggerMixedBin->Fill(ptTrig,eventBin);
3852     if(fCorrelVzBin)fhPtTriggerMixedVzBin->Fill(ptTrig, GetEventVzBin());
3853
3854     //
3855     // Correlation histograms
3856     //
3857     for(Int_t j1 = 0;j1 <nTracks; j1++ )
3858     {
3859       AliAODPWG4Particle *track = (AliAODPWG4Particle*) bgTracks->At(j1) ;
3860
3861       if(!track) continue;
3862
3863       ptAssoc  = track->Pt();
3864       etaAssoc = track->Eta();
3865       phiAssoc = track->Phi() ;
3866       if(phiAssoc < 0) phiAssoc+=TMath::TwoPi();
3867
3868       deltaPhi = phiTrig-phiAssoc;
3869       if(deltaPhi < -TMath::PiOver2())  deltaPhi+=TMath::TwoPi();
3870       if(deltaPhi > 3*TMath::PiOver2()) deltaPhi-=TMath::TwoPi();
3871       deltaEta = etaTrig-etaAssoc;
3872
3873       AliDebug(1,Form("deltaPhi= %f, deltaEta=%f",deltaPhi, deltaEta));
3874
3875       // Angular correlation
3876       fhMixDeltaPhiCharged        ->Fill(ptTrig,   deltaPhi);
3877       fhMixDeltaPhiDeltaEtaCharged->Fill(deltaPhi, deltaEta);
3878
3879       //
3880       // Momentum imbalance
3881       //
3882       if ( (deltaPhi > fDeltaPhiMinCut)   && (deltaPhi < fDeltaPhiMaxCut)   )
3883       {
3884         xE = -ptAssoc/ptTrig*TMath::Cos(deltaPhi); // -(px*pxTrig+py*pyTrig)/(ptTrig*ptTrig);
3885
3886         if(xE < 0.)
3887           AliWarning(Form("Careful!!, negative xE %2.2f for right UE cos(dPhi %2.2f) = %2.2f, check correlation dPhi limits %f to %f",
3888                           xE,deltaPhi*TMath::RadToDeg(),TMath::Cos(deltaPhi),fDeltaPhiMinCut*TMath::RadToDeg(),fDeltaPhiMaxCut*TMath::RadToDeg()));
3889
3890         fhMixXECharged->Fill(ptTrig,xE);
3891         if(xE > 0 ) fhMixHbpXECharged->Fill(ptTrig, TMath::Log(1./xE));
3892       }
3893
3894       //
3895       // Underlying event momentum imbalance
3896       //
3897       if ( (deltaPhi > fUeDeltaPhiMinCut) && (deltaPhi < fUeDeltaPhiMaxCut) )
3898       {
3899         //Underlying event region
3900         Double_t randomphi = gRandom->Uniform(fDeltaPhiMinCut,fDeltaPhiMaxCut);
3901         Double_t uexE = -(ptAssoc/ptTrig)*TMath::Cos(randomphi);
3902
3903         if(uexE < 0.)
3904           AliWarning(Form("Careful!!, negative xE %2.2f for left UE cos(dPhi %2.2f) = %2.2f, check correlation dPhi limits %f to %f",
3905                           uexE,randomphi*TMath::RadToDeg(),TMath::Cos(randomphi),fDeltaPhiMinCut*TMath::RadToDeg(),fDeltaPhiMaxCut*TMath::RadToDeg()));
3906
3907         fhMixXEUeCharged->Fill(ptTrig,uexE);
3908       }
3909
3910       // Set the pt associated bin for the defined bins
3911       Int_t assocBin   = -1;
3912       for(Int_t i = 0 ; i < fNAssocPtBins ; i++)
3913       {
3914         if(ptAssoc > fAssocPtBinLimit[i] && ptAssoc < fAssocPtBinLimit[i+1]) assocBin= i;
3915       }
3916
3917       //
3918       // Assign to the histogram array a bin corresponding to a combination of pTa and vz bins
3919       //
3920       Int_t nz = 1;
3921       Int_t vz = 0;
3922
3923       if(fCorrelVzBin)
3924       {
3925         nz = GetNZvertBin();
3926         vz = GetEventVzBin();
3927       }
3928
3929       Int_t bin = assocBin*nz+vz;
3930
3931       if(bin < 0) continue ; // this pt bin was not considered
3932
3933       fhMixDeltaPhiChargedAssocPtBin        [bin]->Fill(ptTrig,   deltaPhi);
3934       fhMixDeltaPhiDeltaEtaChargedAssocPtBin[bin]->Fill(deltaPhi, deltaEta);
3935
3936       if(fFillEtaGapsHisto)
3937       {
3938         if(TMath::Abs(deltaEta) > 0.8)
3939           fhMixDeltaPhiChargedAssocPtBinDEta08  [bin]->Fill(ptTrig, deltaPhi);
3940         if(TMath::Abs(deltaEta) < 0.01)
3941           fhMixDeltaPhiChargedAssocPtBinDEta0   [bin]->Fill(ptTrig, deltaPhi);
3942       }
3943
3944     } // track loop
3945   } // mixed event loop
3946 }
3947
3948
3949 //_______________________________________________________________________________________________________
3950 void AliAnaParticleHadronCorrelation::MakeNeutralCorrelation(AliAODPWG4ParticleCorrelation * aodParticle)
3951 {
3952   // Neutral Pion Correlation Analysis
3953
3954   TObjArray * pi0list = (TObjArray*) GetAODBranch(fPi0AODBranchName); // For the future, foresee more possible pi0 lists
3955   if(!pi0list) return ;
3956
3957   Int_t npi0 = pi0list->GetEntriesFast();
3958   if(npi0 == 0) return ;
3959
3960   AliDebug(1,Form("Particle - pi0 correlation, %d pi0's",npi0));
3961
3962   Int_t evtIndex11 = 0 ;
3963   Int_t evtIndex12 = 0 ;
3964   if (GetMixedEvent())
3965   {
3966     evtIndex11 = GetMixedEvent()->EventIndexForCaloCluster(aodParticle->GetCaloLabel(0)) ;
3967     evtIndex12 = GetMixedEvent()->EventIndexForCaloCluster(aodParticle->GetCaloLabel(1)) ;
3968   }
3969
3970   Float_t pt   = -100. ;
3971   Float_t zT   = -100. ;
3972   Float_t phi  = -100. ;
3973   Float_t eta  = -100. ;
3974   Float_t xE   = -100. ;
3975   Float_t hbpXE= -100. ;
3976   Float_t hbpZT= -100. ;
3977
3978   Float_t ptTrig  = aodParticle->Pt();
3979   Float_t phiTrig = aodParticle->Phi();
3980   Float_t etaTrig = aodParticle->Eta();
3981   Float_t deltaPhi= -100. ;
3982   Float_t deltaEta= -100. ;
3983
3984   // In case of pi0/eta trigger, we may want to check their decay correlation,
3985   // get their decay children
3986
3987   Bool_t decayFound = kFALSE;
3988   if(fPi0Trigger) decayFound = GetDecayPhotonMomentum(aodParticle->GetCaloLabel(0),aodParticle->GetCaloLabel(1),aodParticle->GetDetectorTag());
3989
3990   TObjArray * refpi0 = 0x0;
3991   Int_t nrefs        = 0;
3992
3993   //Loop on stored AOD pi0
3994
3995   for(Int_t iaod = 0; iaod < npi0 ; iaod++)
3996   {
3997     AliAODPWG4Particle* pi0 =  (AliAODPWG4Particle*) (pi0list->At(iaod));
3998
3999     Int_t evtIndex2 = 0 ;
4000     Int_t evtIndex3 = 0 ;
4001     if (GetMixedEvent())
4002     {
4003       evtIndex2 = GetMixedEvent()->EventIndexForCaloCluster(pi0->GetCaloLabel(0)) ;
4004       evtIndex3 = GetMixedEvent()->EventIndexForCaloCluster(pi0->GetCaloLabel(1)) ;
4005
4006       if (evtIndex11 == evtIndex2 || evtIndex12 == evtIndex2 ||
4007           evtIndex11 == evtIndex3 || evtIndex12 == evtIndex3) // trigger and pi0 are not from different events
4008         continue ;
4009     }
4010
4011     pt  = pi0->Pt();
4012
4013     if(pt < fMinAssocPt || pt > fMaxAssocPt) continue ;
4014
4015     //remove trigger itself for correlation when use charged triggers
4016     if(aodParticle->GetCaloLabel(0) >= 0 &&
4017        (pi0->GetCaloLabel(0) == aodParticle->GetCaloLabel(0) || pi0->GetCaloLabel(1) == aodParticle->GetCaloLabel(0))) continue ;
4018
4019     if( aodParticle->GetCaloLabel(1) >= 0 &&
4020        (pi0->GetCaloLabel(0) == aodParticle->GetCaloLabel(1) || pi0->GetCaloLabel(1) == aodParticle->GetCaloLabel(1))) continue ;
4021
4022     //
4023     // Angular correlations
4024     //
4025     phi      = pi0->Phi() ;
4026     eta      = pi0->Eta() ;
4027     deltaEta = etaTrig-eta;
4028     deltaPhi = phiTrig-phi;
4029     if(deltaPhi <= -TMath::PiOver2()) deltaPhi+=TMath::TwoPi();
4030     if(deltaPhi > 3*TMath::PiOver2()) deltaPhi-=TMath::TwoPi();
4031
4032     fhEtaNeutral     ->Fill(pt    , eta     );
4033     fhPhiNeutral     ->Fill(pt    , phi     );
4034     fhDeltaEtaNeutral->Fill(ptTrig, deltaEta);
4035     fhDeltaPhiNeutral->Fill(ptTrig, deltaPhi);
4036
4037     if(pt > 2 ) fhDeltaPhiDeltaEtaNeutral->Fill(deltaPhi, deltaEta);
4038
4039     //
4040     // Momentum imbalance
4041     //
4042     zT  = pt/ptTrig ;
4043
4044     hbpZT = -100;
4045     hbpXE = -100;
4046
4047     if(zT > 0 ) hbpZT = TMath::Log(1./zT);
4048
4049     //delta phi cut for correlation
4050     if( (deltaPhi > fDeltaPhiMinCut) && ( deltaPhi < fDeltaPhiMaxCut) )
4051     {
4052       xE  =-pt/ptTrig*TMath::Cos(deltaPhi); // -(px*pxTrig+py*pyTrig)/(ptTrig*ptTrig);
4053
4054       if(xE < 0.)
4055         AliWarning(Form("Careful!!, negative xE %2.2f for right UE cos(dPhi %2.2f) = %2.2f, check correlation dPhi limits %f to %f",
4056                         xE,deltaPhi*TMath::RadToDeg(),TMath::Cos(deltaPhi),fDeltaPhiMinCut*TMath::RadToDeg(),fDeltaPhiMaxCut*TMath::RadToDeg()));
4057
4058       if(xE > 0 ) hbpXE = TMath::Log(1./xE);
4059
4060       fhDeltaPhiNeutralPt->Fill(pt,deltaPhi);
4061       fhXENeutral        ->Fill(ptTrig,xE);
4062       fhPtHbpXENeutral   ->Fill(ptTrig,hbpXE);
4063       fhZTNeutral        ->Fill(ptTrig,zT);
4064       fhPtHbpZTNeutral   ->Fill(ptTrig,hbpZT);
4065     }
4066     else if ( (deltaPhi > fUeDeltaPhiMinCut) && (deltaPhi < fUeDeltaPhiMaxCut) )
4067     {
4068       // Randomize angle for xE calculation
4069       Double_t randomphi = gRandom->Uniform(fDeltaPhiMinCut,fDeltaPhiMaxCut);
4070
4071       xE = -(pt/ptTrig)*TMath::Cos(randomphi);
4072       if(xE > 0 ) hbpXE = TMath::Log(1./xE);
4073
4074       fhDeltaPhiUeNeutralPt->Fill(pt,deltaPhi);
4075       fhZTUeNeutral        ->Fill(ptTrig,zT);
4076       fhPtHbpZTUeNeutral   ->Fill(ptTrig,hbpZT);
4077       fhXEUeNeutral        ->Fill(ptTrig,xE);
4078       fhPtHbpXEUeNeutral   ->Fill(ptTrig,hbpXE);
4079     }
4080
4081     // Several UE calculation, not sure it is useful
4082     // with partical calorimter acceptance
4083     if(fMakeSeveralUE) FillNeutralUnderlyingEventSidesHistograms(ptTrig,pt,zT,hbpZT,deltaPhi);
4084
4085     //
4086     // Decay photon correlations
4087     //
4088     if(fPi0Trigger && decayFound)
4089       FillDecayPhotonCorrelationHistograms(pt, phi, kFALSE) ;
4090
4091     if(fFillAODWithReferences)
4092     {
4093       nrefs++;
4094       if(nrefs==1)
4095       {
4096         refpi0 = new TObjArray(0);
4097         refpi0->SetName(GetAODObjArrayName()+"Pi0s");
4098         refpi0->SetOwner(kFALSE);
4099       }
4100       refpi0->Add(pi0);
4101     } // put references in trigger AOD
4102
4103     AliDebug(1,Form("Selected pi0: pt %2.2f, phi %2.2f, eta %2.2f",pt,phi,eta));
4104
4105   }//loop
4106
4107   //Fill AOD with reference tracks, if not filling histograms
4108   if(fFillAODWithReferences && refpi0)
4109   {
4110     aodParticle->AddObjArray(refpi0);
4111   }
4112 }
4113
4114 //__________________________________________________________________________________________________________________
4115 void  AliAnaParticleHadronCorrelation::MakeMCChargedCorrelation(Int_t label, Int_t histoIndex, Bool_t lostDecayPair)
4116 {
4117   // Charged Hadron Correlation Analysis with MC information
4118
4119   AliDebug(1,"Make trigger particle - charged hadron correlation in AOD MC level");
4120
4121   if( label < 0 )
4122   {
4123     AliDebug(1,Form(" *** bad label ***:  label %d", label));
4124     return;
4125   }
4126
4127   // Do MC correlation for a given particle type range.
4128   // Types defined in GetMCTagHistogramIndex:
4129   // 0 direct gamma; 1 pi0; 2 pi0 decay; 3 eta decay; 4 other decay; 5 electron; 6 other (hadron)
4130   if(histoIndex < fMCGenTypeMin || histoIndex > fMCGenTypeMax) return ;
4131
4132   AliStack         * stack        = 0x0 ;
4133   TParticle        * primary      = 0x0 ;
4134   TClonesArray     * mcparticles  = 0x0 ;
4135   AliAODMCParticle * aodprimary   = 0x0 ;
4136
4137   Double_t eprim   = 0 ;
4138   Double_t ptprim  = 0 ;
4139   Double_t phiprim = 0 ;
4140   Double_t etaprim = 0 ;
4141   Int_t    nTracks = 0 ;
4142   Int_t iParticle  = 0 ;
4143
4144   Bool_t leadTrig = kTRUE;
4145
4146   if( GetReader()->ReadStack() )
4147   {
4148     stack =  GetMCStack() ;
4149     if(!stack)
4150     {
4151       AliFatal("Stack not available, is the MC handler called? STOP");
4152       return;
4153     }
4154
4155     //nTracks = stack->GetNtrack() ;
4156     nTracks = stack->GetNprimary();
4157     if( label >=  stack->GetNtrack() )
4158     {
4159       if(GetDebug() > 2)
4160         AliInfo(Form("*** large label ***:  label %d, n tracks %d", label, stack->GetNtrack()));
4161       return ;
4162     }
4163
4164     primary = stack->Particle(label);
4165     if ( !primary )
4166     {
4167       AliInfo(Form(" *** no primary ***:  label %d", label));
4168       return;
4169     }
4170
4171     eprim    = primary->Energy();
4172     ptprim   = primary->Pt();
4173     etaprim  = primary->Eta();
4174     phiprim  = primary->Phi();
4175     if(phiprim < 0) phiprim+=TMath::TwoPi();
4176
4177     if(ptprim < 0.01 || eprim < 0.01) return ;
4178
4179     for (iParticle = 0 ; iParticle <  nTracks ; iParticle++)
4180     {
4181       TParticle * particle = stack->Particle(iParticle);
4182
4183       //keep only final state particles
4184       if( particle->GetStatusCode() != 1 ) continue ;
4185
4186       //---------- Charged particles ----------------------
4187       Int_t pdg    = particle->GetPdgCode();
4188       Int_t charge = (Int_t) TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
4189       if(charge == 0) continue;
4190
4191       particle->Momentum(fMomentum);
4192
4193       //Particles in CTS acceptance, make sure to use the same selection as in the reader
4194       Bool_t inCTS =  GetReader()->GetFiducialCut()->IsInFiducialCut(fMomentum.Eta(),fMomentum.Phi(),kCTS);
4195       //printf("Accepted? %d; pt %2.2f, eta %2.2f, phi %2.2f\n",inCTS,fMomentum.Pt(),fMomentum.Eta(),fMomentum.Phi()*TMath::RadToDeg());
4196       if( !inCTS ) continue;
4197
4198       // Remove conversions
4199       if ( TMath::Abs(pdg) == 11 && stack->Particle(particle->GetFirstMother())->GetPdgCode() == 22 ) continue ;
4200
4201       if ( label == iParticle ) continue; // avoid trigger particle
4202
4203       Float_t phi = particle->Phi();
4204       if(phi < 0) phi+=TMath::TwoPi();
4205
4206       Bool_t lead = FillChargedMCCorrelationHistograms(particle->Pt(),phi,particle->Eta(),ptprim,phiprim,etaprim,histoIndex,lostDecayPair);
4207       if(!lead) leadTrig = kFALSE;
4208       //if ( !lead && (fMakeAbsoluteLeading || fMakeNearSideLeading) ) return;
4209
4210     } //track loop
4211
4212   } //ESD MC
4213
4214   else if( GetReader()->ReadAODMCParticles() )
4215   {
4216     //Get the list of MC particles
4217     mcparticles = GetReader()->GetAODMCParticles();
4218     if( !mcparticles ) return;
4219
4220     nTracks = mcparticles->GetEntriesFast() ;
4221
4222     if( label >= nTracks )
4223     {
4224       if(GetDebug() > 2)
4225         AliInfo(Form(" *** large label ***:  label %d, n tracks %d", label,nTracks));
4226       return;
4227     }
4228
4229     //Get the particle
4230     aodprimary = (AliAODMCParticle*) mcparticles->At(label);
4231     if( !aodprimary )
4232     {
4233       AliInfo(Form(" *** no AOD primary ***:  label %d", label));
4234       return;
4235     }
4236
4237     eprim   = aodprimary->E();
4238     ptprim  = aodprimary->Pt();
4239     etaprim = aodprimary->Eta();
4240     phiprim = aodprimary->Phi();
4241     if(phiprim < 0) phiprim+=TMath::TwoPi();
4242
4243     if(ptprim < 0.01 || eprim < 0.01) return ;
4244
4245     for (iParticle = 0; iParticle < nTracks; iParticle++)
4246     {
4247       AliAODMCParticle *part = (AliAODMCParticle*) mcparticles->At(iParticle);
4248
4249       if (!part->IsPhysicalPrimary() ) continue; // same as part->GetStatus() !=1
4250
4251       if ( part->Charge() == 0 ) continue;
4252
4253       fMomentum.SetPxPyPzE(part->Px(),part->Py(),part->Pz(),part->E());
4254
4255       //Particles in CTS acceptance, make sure to use the same selection as in the reader
4256       Bool_t inCTS =  GetReader()->GetFiducialCut()->IsInFiducialCut(fMomentum.Eta(),fMomentum.Phi(),kCTS);
4257       //printf("Accepted? %d; pt %2.2f, eta %2.2f, phi %2.2f\n",inCTS,fMomentum.Pt(),fMomentum.Eta(),fMomentum.Phi()*TMath::RadToDeg());
4258       if( !inCTS ) continue;
4259
4260       // Remove conversions
4261       Int_t indexmother = part->GetMother();
4262       if ( indexmother > -1 )
4263       {
4264         Int_t pdg = part->GetPdgCode();
4265         Int_t mPdg = ((AliAODMCParticle*) mcparticles->At(indexmother)) ->GetPdgCode();
4266         if (TMath::Abs(pdg) == 11 && mPdg == 22) continue;
4267       }
4268
4269       if ( label == iParticle ) continue; // avoid trigger particle
4270
4271       Float_t phi = part->Phi();
4272       if(phi < 0) phi+=TMath::TwoPi();
4273
4274       Bool_t lead = FillChargedMCCorrelationHistograms(part->Pt(),phi,part->Eta(),ptprim,phiprim,etaprim, histoIndex,lostDecayPair);
4275       if(!lead) leadTrig = kFALSE;
4276       //if ( !lead && (fMakeAbsoluteLeading || fMakeNearSideLeading)) return;
4277
4278     }  //MC particle loop
4279   }// AOD MC
4280
4281   // Trigger MC particle histograms
4282   //if (!lead  && (fMakeAbsoluteLeading || fMakeNearSideLeading)) return;
4283
4284   fhMCPtTrigger [histoIndex]->Fill(ptprim);
4285   fhMCPhiTrigger[histoIndex]->Fill(ptprim,phiprim);
4286   fhMCEtaTrigger[histoIndex]->Fill(ptprim,etaprim);
4287
4288   if(lostDecayPair)
4289   {
4290     // check index in GetMCTagIndexHistogram
4291     if     (histoIndex == 2 && 8 >= fMCGenTypeMin && 8 <= fMCGenTypeMax )
4292     {
4293       // pi0 decay
4294       fhMCPtTrigger [8]->Fill(ptprim);
4295       fhMCPhiTrigger[8]->Fill(ptprim,phiprim);
4296       fhMCEtaTrigger[8]->Fill(ptprim,etaprim);
4297     }
4298     else if(histoIndex == 4 && 9 >= fMCGenTypeMin && 9 <= fMCGenTypeMax )
4299     {
4300       // eta decay
4301       fhMCPtTrigger [9]->Fill(ptprim);
4302       fhMCPhiTrigger[9]->Fill(ptprim,phiprim);
4303       fhMCEtaTrigger[9]->Fill(ptprim,etaprim);
4304     }
4305   }
4306
4307   if(!leadTrig && (fMakeAbsoluteLeading || fMakeNearSideLeading) )
4308   {
4309     AliDebug(1,Form("Not leading primary trigger: pT %2.2f, phi %2.2f, eta %2.2f",
4310                     ptprim,phiprim*TMath::RadToDeg(),etaprim));
4311
4312     fhMCPtTriggerNotLeading [histoIndex]->Fill(ptprim);
4313     fhMCPhiTriggerNotLeading[histoIndex]->Fill(ptprim,phiprim);
4314     fhMCEtaTriggerNotLeading[histoIndex]->Fill(ptprim,etaprim);
4315
4316     if(lostDecayPair)
4317     {
4318       // check index in GetMCTagIndexHistogram
4319       if      (histoIndex == 2  && 8 >= fMCGenTypeMin && 8 <= fMCGenTypeMax )
4320       {
4321         // pi0 decay
4322         fhMCPtTriggerNotLeading [8]->Fill(ptprim);
4323         fhMCPhiTriggerNotLeading[8]->Fill(ptprim,phiprim);
4324         fhMCEtaTriggerNotLeading[8]->Fill(ptprim,etaprim);
4325       }
4326       else  if(histoIndex == 4  && 9 >= fMCGenTypeMin && 9 <= fMCGenTypeMax )
4327       {
4328         // eta decay
4329         fhMCPtTriggerNotLeading [9]->Fill(ptprim);
4330         fhMCPhiTriggerNotLeading[9]->Fill(ptprim,phiprim);
4331         fhMCEtaTriggerNotLeading[9]->Fill(ptprim,etaprim);
4332       }
4333     }
4334   }
4335 }
4336
4337 //_____________________________________________________________________
4338 void AliAnaParticleHadronCorrelation::Print(const Option_t * opt) const
4339 {
4340
4341   //Print some relevant parameters set for the analysis
4342   if(! opt)
4343     return;
4344
4345   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
4346   AliAnaCaloTrackCorrBaseClass::Print(" ");
4347   printf("Pt trigger > %2.2f; < %2.2f\n", GetMinPt() , GetMaxPt()) ;
4348   printf("Pt associa > %2.2f; < %2.2f\n", fMinAssocPt, fMaxAssocPt) ;
4349   printf("Phi trigger particle-Hadron   <  %3.2f\n", fDeltaPhiMaxCut) ;
4350   printf("Phi trigger particle-Hadron   >  %3.2f\n", fDeltaPhiMinCut) ;
4351   printf("Phi trigger particle-UeHadron <  %3.2f\n", fUeDeltaPhiMaxCut) ;
4352   printf("Phi trigger particle-UeHadron >  %3.2f\n", fUeDeltaPhiMinCut) ;
4353   printf("Isolated Trigger?  %d\n"     , fSelectIsolated) ;
4354   printf("Several UE?  %d\n"           , fMakeSeveralUE) ;
4355   printf("Name of AOD Pi0 Branch %s \n", fPi0AODBranchName.Data());
4356   printf("Do Decay-hadron correlation ?  %d\n", fPi0Trigger) ;
4357   printf("Select absolute leading for cluster triggers ?  %d or Near Side %d\n",
4358          fMakeAbsoluteLeading, fMakeNearSideLeading) ;
4359   printf("Trigger pt bins  %d\n", fNAssocPtBins) ;
4360   for (Int_t ibin = 0; ibin<fNAssocPtBins; ibin++) {
4361     printf("\t bin %d = [%2.1f,%2.1f]\n", ibin, fAssocPtBinLimit[ibin], fAssocPtBinLimit[ibin+1]) ;
4362   }
4363
4364 }
4365
4366 //____________________________________________________________
4367 void AliAnaParticleHadronCorrelation::SetNAssocPtBins(Int_t n)
4368 {
4369   // Set number of bins
4370
4371   fNAssocPtBins  = n ;
4372
4373   if(n < 20 && n > 0)
4374   {
4375     fNAssocPtBins  = n ;
4376   }
4377   else
4378   {
4379     AliWarning("n = larger than 19 or too small, set to 19");
4380     fNAssocPtBins = 19;
4381   }
4382 }
4383
4384 //______________________________________________________________________________
4385 void AliAnaParticleHadronCorrelation::SetAssocPtBinLimit(Int_t ibin, Float_t pt)
4386 {
4387   // Set the list of limits for the trigger pt bins
4388
4389   if(ibin <= fNAssocPtBins || ibin >= 0)
4390   {
4391     fAssocPtBinLimit[ibin] = pt  ;
4392   }
4393   else
4394   {
4395     AliWarning(Form("Bin  number too large %d > %d or small, nothing done", ibin, fNAssocPtBins)) ;
4396   }
4397 }
4398