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