PWGGA/CaloTrackCorrelations/AliAnaParticleHadronCorrelation.cxx

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