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