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