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