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