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