PWGGA/CaloTrackCorrelations/AliAnaParticleHadronCorrelation.cxx

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