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