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