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