PWGGA/CaloTrackCorrelations/AliAnaParticleHadronCorrelation.cxx

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