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