PWGGA/CaloTrackCorrelations/AliAnaParticleHadronCorrelation.cxx

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