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