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