PWGGA/CaloTrackCorrelations/AliAnaParticleHadronCorrelation.cxx

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