]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGGA/CaloTrackCorrelations/AliAnaParticleHadronCorrelation.cxx
add dphi histograms for deta > 0.8
[u/mrichter/AliRoot.git] / PWGGA / CaloTrackCorrelations / AliAnaParticleHadronCorrelation.cxx
index 2fca91d2ffb118d6586acb43316b38542507e1d8..e25ceb253707c5231887ea0f45636f225e077772 100755 (executable)
 
 // --- ROOT system ---
 //#include "TClonesArray.h"
-#include "TClass.h"
-#include "TMath.h"
-#include "TH3D.h"
-#include "TDatabasePDG.h"
+#include <TClass.h>
+#include <TMath.h>
+#include <TH2F.h>
+#include <TDatabasePDG.h>
 
 //---- ANALYSIS system ----
 #include "AliNeutralMesonSelection.h" 
 #include "AliAnaParticleHadronCorrelation.h" 
 #include "AliCaloTrackReader.h"
-#include "AliCaloPID.h"
 #include "AliAODPWG4ParticleCorrelation.h"
 #include "AliFiducialCut.h"
 #include "AliVTrack.h"
@@ -51,6 +50,9 @@
 #include "AliStack.h"
 #include "AliAODMCParticle.h"
 #include "AliMixedEvent.h"
+#include "AliAnalysisManager.h"
+#include "AliInputEventHandler.h"
+#include "AliEventplane.h"
 
 ClassImp(AliAnaParticleHadronCorrelation)
 
@@ -64,16 +66,20 @@ ClassImp(AliAnaParticleHadronCorrelation)
     fSelectIsolated(0),             fMakeSeveralUE(0),              
     fUeDeltaPhiMaxCut(0.),          fUeDeltaPhiMinCut(0.), 
     fPi0AODBranchName(""),          fNeutralCorr(0),       
-    fPi0Trigger(0),                 
+    fPi0Trigger(0),                 fDecayTrigger(0),
     fMakeAbsoluteLeading(0),        fMakeNearSideLeading(0),      
-    fLeadingTriggerIndex(-1),       
+    fLeadingTriggerIndex(-1),       fHMPIDCorrelation(0),  fFillBradHisto(0),
     fNAssocPtBins(0),               fAssocPtBinLimit(),
+    fListMixEvents(),               fUseMixStoredInReader(0),
     //Histograms
     fhPtLeading(0),                 fhPhiLeading(0),       
-    fhEtaLeading(0),                fhDeltaPhiDeltaEtaCharged(0),
+    fhEtaLeading(0),                
+    fhPtLeadingCentrality(0),       fhPtLeadingEventPlane(0), 
+    fhLeadingEventPlaneCentrality(0),fhDeltaPhiDeltaEtaCharged(0),
     fhPhiCharged(0),                fhEtaCharged(0), 
     fhDeltaPhiCharged(0),           fhDeltaEtaCharged(0), 
     fhDeltaPhiChargedPt(0),         fhDeltaPhiUeChargedPt(0), 
+    fhUePart(0),
     fhXECharged(0),                 fhXEUeCharged(0),
     fhXEPosCharged(0),              fhXENegCharged(0),
     fhPtHbpXECharged(0),            fhPtHbpXEUeCharged(0),
@@ -89,11 +95,11 @@ ClassImp(AliAnaParticleHadronCorrelation)
     fhTrigDeltaPhiCharged(0x0),     fhTrigDeltaEtaCharged(0x0),
     fhTrigXECorr(0x0),              fhTrigXEUeCorr(0x0),
     fhTrigZTCorr(0x0),              fhTrigZTUeCorr(0x0),
-    fhAssocPt(0),                   fhAssocPtBkg(0),  
-    fhDeltaPhiAssocPtBin(0),        fhDeltaPhiAssocPtBinHMPID(0),fhDeltaPhiAssocPtBinHMPIDAcc(0),
+    fhAssocPtBkg(0),  
+    fhDeltaPhiAssocPtBin(0),        fhDeltaPhiAssocPtBinDEta08(0), 
+    fhDeltaPhiAssocPtBinHMPID(0),   fhDeltaPhiAssocPtBinHMPIDAcc(0),
     fhDeltaPhiBradAssocPtBin(0),    fhDeltaPhiBrad(0),
-    fhXEAssocPtBin(0),              fhXE(0),
-    fhZTAssocPtBin(0),              fhZT(0),
+    fhXEAssocPtBin(0),              fhZTAssocPtBin(0),             
     fhDeltaPhiDeltaEtaNeutral(0), 
     fhPhiNeutral(0),                fhEtaNeutral(0), 
     fhDeltaPhiNeutral(0),           fhDeltaEtaNeutral(0),
@@ -110,22 +116,542 @@ ClassImp(AliAnaParticleHadronCorrelation)
     fhPtPi0DecayRatio(0),
     fhDeltaPhiDecayCharged(0),      fhXEDecayCharged(0), fhZTDecayCharged(0), 
     fhDeltaPhiDecayNeutral(0),      fhXEDecayNeutral(0), fhZTDecayNeutral(0),
-    fh2phiLeadingParticle(0x0),
-    fhMCLeadingCount(0),
+    fhDeltaPhiDecayChargedAssocPtBin(0), 
+    fhXEDecayChargedAssocPtBin(0),  fhZTDecayChargedAssocPtBin(0),                
+    fh2phiLeadingParticle(0x0),     fhMCPtLeading(0),
     fhMCEtaCharged(0),              fhMCPhiCharged(0), 
     fhMCDeltaEtaCharged(0),         fhMCDeltaPhiCharged(0x0),
     fhMCDeltaPhiDeltaEtaCharged(0), fhMCDeltaPhiChargedPt(0),
-    fhMCPtXECharged(0),             fhMCPtHbpXECharged(0),          
+    fhMCPtXECharged(0),             fhMCPtXEUeCharged(0),
+    fhMCPtHbpXECharged(0),          fhMCPtHbpXEUeCharged(0),
+    fhMCUePart(0),
     fhMCPtZTCharged(0),             fhMCPtHbpZTCharged(0),
-    fhMCPtTrigPout(0),
-    fhMCPtAssocDeltaPhi(0)
-{
+    fhMCPtTrigPout(0),              fhMCPtAssocDeltaPhi(0),
+    //Mixing
+    fhNEventsTrigger(0),            
+    fhNtracksAll(0),                fhNtracksTrigger(0),            
+    fhNtracksMB(0),               
+    fhMixDeltaPhiCharged(0),        fhMixDeltaPhiDeltaEtaCharged(0),
+    fhMixDeltaPhiChargedAssocPtBin(),
+    fhMixDeltaPhiChargedAssocPtBinDEta08(),
+    fhMixDeltaPhiDeltaEtaChargedAssocPtBin(),
+    fhEventBin(0),                  fhEventMixBin(0)
+{ 
   //Default Ctor
-  
+    
   //Initialize parameters
   InitParameters();
 }
 
+//_________________________________________________________________
+AliAnaParticleHadronCorrelation::~AliAnaParticleHadronCorrelation() 
+{
+  // Remove event containers
+  
+  if(DoOwnMix() && fListMixEvents)
+  {      
+    for(Int_t iz=0; iz < GetNZvertBin(); iz++)
+    {
+      for(Int_t ic=0; ic < GetNCentrBin(); ic++)
+      {
+         for(Int_t irp=0; irp<GetNRPBin(); irp++)
+         {
+           Int_t bin = GetEventMixBin(ic, iz, irp);
+           fListMixEvents[bin]->Delete() ;
+           delete fListMixEvents[bin] ;
+         }
+      }
+    }
+    
+    delete[] fListMixEvents;
+    
+  }
+}
+
+//______________________________________________________________________________________________________________________________________________________
+void AliAnaParticleHadronCorrelation::FillChargedAngularCorrelationHistograms(const Float_t ptAssoc,  const Float_t ptTrig,   const Int_t   assocBin,
+                                                                              const Float_t phiAssoc, const Float_t phiTrig,  Float_t &     deltaPhi,
+                                                                              const Float_t etaAssoc, const Float_t etaTrig,  
+                                                                              const Bool_t  decay,    const Float_t hmpidSignal, const Int_t nTracks)
+{
+  // Fill angular correlation related histograms
+  
+  Float_t deltaEta    = etaTrig-etaAssoc;
+          deltaPhi    = phiTrig-phiAssoc;
+  Float_t deltaPhiOrg = deltaPhi;
+  
+  if(deltaPhi <= -TMath::PiOver2()) deltaPhi+=TMath::TwoPi();
+  if(deltaPhi > 3*TMath::PiOver2()) deltaPhi-=TMath::TwoPi();
+  
+  fhEtaCharged     ->Fill(ptAssoc,etaAssoc);
+  fhPhiCharged     ->Fill(ptAssoc,phiAssoc);
+  fhDeltaEtaCharged->Fill(ptTrig ,deltaEta);
+  fhDeltaPhiCharged->Fill(ptTrig ,deltaPhi);
+  
+  if(ptAssoc > 2 )           fhDeltaPhiDeltaEtaCharged->Fill(deltaPhi, deltaEta);
+  
+  if(fDecayTrigger && decay) fhDeltaPhiDecayCharged   ->Fill(ptTrig  , deltaPhi);      
+  
+  Double_t  dphiBrad = -100;
+  if(fFillBradHisto)
+  {
+    dphiBrad = atan2(sin(deltaPhiOrg), cos(deltaPhiOrg))/TMath::Pi();//-1 to 1
+    if(TMath::Abs(dphiBrad)>0.325 && TMath::Abs(dphiBrad)<0.475)  //Hardcoded values, BAD, FIXME
+    {
+      fhAssocPtBkg->Fill(ptTrig, ptAssoc);
+    }
+    
+    if(dphiBrad<-1./3) dphiBrad += 2;
+    fhDeltaPhiBrad->Fill(ptTrig, dphiBrad);
+  }
+  
+  // Fill histograms in bins of associated particle pT
+  if(assocBin>=0)
+  {
+      fhDeltaPhiAssocPtBin            [assocBin]->Fill(ptTrig, deltaPhi);
+    
+    if(TMath::Abs(deltaEta)> 0.8) 
+      fhDeltaPhiAssocPtBinDEta08      [assocBin]->Fill(ptTrig, deltaPhi);
+    
+    if (fFillBradHisto)
+      fhDeltaPhiBradAssocPtBin        [assocBin]->Fill(ptTrig, dphiBrad);
+    
+    if(fDecayTrigger && decay)
+      fhDeltaPhiDecayChargedAssocPtBin[assocBin]->Fill(ptTrig, deltaPhi);      
+    
+    if(fHMPIDCorrelation)
+    {
+      if( hmpidSignal > 0 )
+      {
+        //printf("Track pt %f with HMPID signal %f \n",pt,hmpidSignal);
+        fhDeltaPhiAssocPtBinHMPID[assocBin]->Fill(ptTrig, deltaPhi);        
+      }
+      
+      if(phiAssoc > 5*TMath::DegToRad() && phiAssoc < 20*TMath::DegToRad())
+      {
+        //printf("Track pt %f in HMPID acceptance phi %f \n ",pt,phi*TMath::RadToDeg() );
+        fhDeltaPhiAssocPtBinHMPIDAcc[assocBin]->Fill(ptTrig, deltaPhi);        
+      }
+    }
+  }
+  
+  //fill different multiplicity histogram
+  if(DoEventSelect())
+  {
+    for(Int_t im = 0; im<GetMultiBin(); im++)
+    {
+      if(nTracks < ( GetMaxMulti() - GetMinMulti() )/GetMultiBin()*(im+1))
+      {
+        fhTrigDeltaPhiCharged[im]->Fill(ptTrig,deltaPhi);
+        fhTrigDeltaEtaCharged[im]->Fill(ptTrig,deltaEta);
+      }
+    }
+  }  
+}
+
+//____________________________________________________________________________________________________________________________________________________
+Bool_t AliAnaParticleHadronCorrelation::FillChargedMCCorrelationHistograms(const Float_t mcAssocPt,      Float_t mcAssocPhi, const Float_t mcAssocEta,
+                                                                           const Float_t mcTrigPt, const Float_t mcTrigPhi,  const Float_t mcTrigEta)
+{
+  // Fill MC histograms independently of AOD or ESD
+  
+  //Select only hadrons in pt range
+  if(mcAssocPt < fMinAssocPt || mcAssocPt > fMaxAssocPt) return kTRUE ; // exclude but continue
+  
+  if(mcAssocPhi < 0) mcAssocPhi+=TMath::TwoPi();    
+  
+  //remove trigger itself for correlation when use charged triggers 
+  if(TMath::Abs(mcAssocPt -mcTrigPt ) < 1e-6 && 
+     TMath::Abs(mcAssocPhi-mcTrigPhi) < 1e-6 && 
+     TMath::Abs(mcAssocEta-mcTrigEta) < 1e-6)            return kTRUE ; // exclude but continue       
+  
+  // Absolute leading?
+  if( fMakeAbsoluteLeading && mcAssocPt > mcTrigPt )     return kFALSE; // jump event
+  
+  //jump out this event if near side associated partile pt larger than trigger
+  if( fMakeNearSideLeading && mcAssocPt > mcTrigPt && 
+     TMath::Abs(mcAssocPhi-mcTrigPhi)<TMath::PiOver2() ) return kFALSE; // jump event
+  
+  Float_t mcdeltaPhi= mcTrigPhi-mcAssocPhi; 
+  if(mcdeltaPhi <= -TMath::PiOver2()) mcdeltaPhi+=TMath::TwoPi();
+  if(mcdeltaPhi > 3*TMath::PiOver2()) mcdeltaPhi-=TMath::TwoPi();            
+  
+  Float_t mcxE    =-mcAssocPt/mcTrigPt*TMath::Cos(mcdeltaPhi);// -(mcAssocPx*pxprim+mcAssocPy*pyprim)/(mcTrigPt*mcTrigPt);  
+  Float_t mchbpXE =-100 ;
+  if(mcxE > 0 ) mchbpXE = TMath::Log(1./mcxE);
+  
+  Float_t mczT    = mcAssocPt/mcTrigPt ;
+  Float_t mchbpZT =-100 ;
+  if(mczT > 0 ) mchbpZT = TMath::Log(1./mczT);
+  
+  //Selection within angular range
+  if( mcdeltaPhi< -TMath::PiOver2())  mcdeltaPhi+=TMath::TwoPi();
+  if( mcdeltaPhi>3*TMath::PiOver2())  mcdeltaPhi-=TMath::TwoPi();              
+  
+  Double_t mcpout = mcAssocPt*TMath::Sin(mcdeltaPhi) ; 
+  
+  if(GetDebug() > 0 )  
+    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",
+           mcAssocPt,mcAssocPhi, mcTrigPhi,fDeltaPhiMinCut, mcdeltaPhi, fDeltaPhiMaxCut, GetMinPt());              
+  
+  // Fill Histograms
+  fhMCEtaCharged     ->Fill(mcAssocPt, mcAssocEta);
+  fhMCPhiCharged     ->Fill(mcAssocPt, mcAssocPhi);
+  fhMCDeltaEtaCharged->Fill(mcTrigPt,    mcTrigEta-mcAssocEta);
+  fhMCDeltaPhiCharged->Fill(mcTrigPt,    mcdeltaPhi);
+  fhMCPtAssocDeltaPhi->Fill(mcAssocPt, mcdeltaPhi);
+  
+  fhMCDeltaPhiDeltaEtaCharged->Fill(mcdeltaPhi,mcTrigEta-mcAssocEta);
+  
+  //delta phi cut for correlation
+  if( (mcdeltaPhi > fDeltaPhiMinCut) && (mcdeltaPhi < fDeltaPhiMaxCut) ) 
+  {
+    fhMCDeltaPhiChargedPt->Fill(mcAssocPt,mcdeltaPhi);
+    fhMCPtXECharged      ->Fill(mcTrigPt,mcxE); 
+    fhMCPtHbpXECharged   ->Fill(mcTrigPt,mchbpXE);
+    fhMCPtZTCharged      ->Fill(mcTrigPt,mczT); 
+    fhMCPtHbpZTCharged   ->Fill(mcTrigPt,mchbpZT);
+    fhMCPtTrigPout       ->Fill(mcTrigPt, mcpout) ;
+  }
+
+  //underlying event
+  if ( (mcdeltaPhi > fUeDeltaPhiMinCut) && (mcdeltaPhi < fUeDeltaPhiMaxCut) )   
+    {
+       Double_t randomphi = gRandom->Uniform(TMath::Pi()/2,3*TMath::Pi()/2);
+       Double_t mcUexE = -(mcAssocPt/mcTrigPt)*TMath::Cos(randomphi);
+  
+       if(mcUexE < 0.) mcUexE = -mcUexE;
+    
+       fhMCPtXEUeCharged->Fill(mcTrigPt,mcUexE);
+       if(mcUexE > 0) fhMCPtHbpXEUeCharged->Fill(mcTrigPt,TMath::Log(1/mcUexE));
+
+       fhMCUePart->Fill(mcTrigPt);
+    }
+  
+  return kTRUE;
+} 
+
+//___________________________________________________________________________________________________________________
+void AliAnaParticleHadronCorrelation::FillChargedMomentumImbalanceHistograms(const Float_t ptTrig,   const Float_t ptAssoc, 
+                                                                             const Float_t xE,       const Float_t hbpXE, 
+                                                                             const Float_t zT,       const Float_t hbpZT, 
+                                                                             const Float_t pout,     const Float_t deltaPhi,
+                                                                             const Int_t   nTracks,  const Int_t   charge,
+                                                                             const Int_t   assocBin, const Bool_t  decay    )
+
+{
+  // Fill mostly momentum imbalance related histograms
+  
+  fhDeltaPhiChargedPt ->Fill(ptAssoc, deltaPhi);
+  fhXECharged         ->Fill(ptTrig , xE); 
+  fhPtHbpXECharged    ->Fill(ptTrig , hbpXE);
+  fhZTCharged         ->Fill(ptTrig , zT); 
+  fhPtHbpZTCharged    ->Fill(ptTrig , hbpZT);
+  fhPtTrigPout        ->Fill(ptTrig , pout) ;
+  fhPtTrigCharged     ->Fill(ptTrig , ptAssoc) ;
+  
+  if(fDecayTrigger && decay)
+  {          
+    fhXEDecayCharged->Fill(ptTrig,xE); 
+    fhZTDecayCharged->Fill(ptTrig,zT);
+  } // photon decay pi0/eta trigger        
+  
+  if(assocBin >= 0 )//away side 
+  {
+    fhXEAssocPtBin[assocBin]->Fill(ptTrig, xE) ;
+    fhZTAssocPtBin[assocBin]->Fill(ptTrig, zT) ;
+    
+    if(fDecayTrigger && decay)
+    {          
+      fhXEDecayChargedAssocPtBin[assocBin]->Fill(ptTrig, xE); 
+      fhZTDecayChargedAssocPtBin[assocBin]->Fill(ptTrig, zT);
+    }
+  }        
+  
+  if(charge > 0)
+  {
+    fhXEPosCharged->Fill(ptTrig,xE) ;
+    fhZTPosCharged->Fill(ptTrig,zT) ;
+  }
+  else  
+  {
+    fhXENegCharged->Fill(ptTrig,xE) ;
+    fhZTNegCharged->Fill(ptTrig,zT) ;
+  }
+  
+  //fill different multiplicity histogram
+  if(DoEventSelect())
+  {
+    for(Int_t im=0; im<GetMultiBin(); im++)
+    {
+      if(nTracks < ( GetMaxMulti() - GetMinMulti() )/GetMultiBin()*(im+1))
+      {
+        fhTrigXECorr[im]->Fill(ptTrig,xE);
+        fhTrigZTCorr[im]->Fill(ptTrig,zT);
+      }
+    }
+  } //multiplicity events selection
+} 
+
+//_______________________________________________________________________________________________________________________
+void AliAnaParticleHadronCorrelation::FillChargedUnderlyingEventHistograms(const Float_t ptTrig,   const Float_t ptAssoc,
+                                                                           const Float_t deltaPhi, const Int_t nTracks)
+{
+  // Fill underlying event histograms
+  
+  fhDeltaPhiUeChargedPt->Fill(ptAssoc,deltaPhi);
+  
+  Double_t randomphi = gRandom->Uniform(TMath::Pi()/2,3*TMath::Pi()/2);
+  Double_t uexE = -(ptAssoc/ptTrig)*TMath::Cos(randomphi);
+  Double_t uezT =   ptAssoc/ptTrig;
+  
+  if(uexE < 0.) uexE = -uexE;
+    
+  fhXEUeCharged->Fill(ptTrig,uexE);
+  if(uexE > 0) fhPtHbpXEUeCharged->Fill(ptTrig,TMath::Log(1/uexE));
+  
+  fhZTUeCharged->Fill(ptTrig,uezT);
+  if(uexE > 0) fhPtHbpZTUeCharged->Fill(ptTrig,TMath::Log(1/uezT));
+  
+  if(DoEventSelect())
+  {
+    for(Int_t im=0; im<GetMultiBin(); im++)
+    {
+      if(nTracks < ( GetMaxMulti() - GetMinMulti() )/GetMultiBin()*(im+1))
+      {
+        fhTrigXEUeCorr[im]->Fill(ptTrig,uexE); // xE? CHECK
+        fhTrigZTUeCorr[im]->Fill(ptTrig,uezT); // zT? CHECK
+      }
+    }
+  } //multiplicity events selection
+}
+
+//___________________________________________________________________________________________________________________________
+void AliAnaParticleHadronCorrelation::FillChargedUnderlyingEventSidesHistograms(const Float_t ptTrig,   const Float_t ptAssoc, 
+                                                                                const Float_t xE,       const Float_t hbpXE, 
+                                                                                const Float_t zT,       const Float_t hbpZT, 
+                                                                                const Float_t deltaPhi)
+{
+  // Fill underlying event histograms to the left and right of trigger
+  
+  if((deltaPhi<-fUeDeltaPhiMinCut) && (deltaPhi >-fUeDeltaPhiMaxCut))
+  {  
+    fhDeltaPhiUeLeftCharged->Fill(ptAssoc, deltaPhi);
+    fhXEUeLeftCharged      ->Fill(ptTrig , xE);
+    fhPtHbpXEUeLeftCharged ->Fill(ptTrig , hbpXE);
+    fhZTUeLeftCharged      ->Fill(ptTrig , zT);
+    fhPtHbpZTUeLeftCharged ->Fill(ptTrig , hbpZT);
+  }
+  
+  if((deltaPhi > fUeDeltaPhiMinCut) && (deltaPhi < fUeDeltaPhiMaxCut))
+  {  
+    fhDeltaPhiUeRightCharged->Fill(ptAssoc, deltaPhi);
+    fhXEUeRightCharged      ->Fill(ptTrig , xE);
+    fhPtHbpXEUeRightCharged ->Fill(ptTrig , hbpXE);
+    fhZTUeRightCharged      ->Fill(ptTrig , zT);
+    fhPtHbpZTUeRightCharged ->Fill(ptTrig , hbpZT);
+  }
+} 
+
+//______________________________________________________________________________________________________________________________
+void AliAnaParticleHadronCorrelation::FillDecayPhotonCorrelationHistograms(const Float_t ptAssoc,     const Float_t phiAssoc, 
+                                                                           const TLorentzVector mom1, const TLorentzVector mom2,
+                                                                           const Bool_t bChargedOrNeutral)
+{
+  // Do correlation with decay photons of triggered pi0 or eta
+  
+  // Calculate the correlation parameters
+  Float_t ptDecay1 = mom1.Pt();
+  Float_t ptDecay2 = mom2.Pt();
+  
+  Float_t zTDecay1 = -100, zTDecay2 = -100;
+  if(ptDecay1) zTDecay1 = ptAssoc/ptDecay1 ;
+  if(ptDecay2) zTDecay2 = ptAssoc/ptDecay2 ;
+  
+  Float_t deltaPhiDecay1 = mom1.Phi()-phiAssoc;
+  if(deltaPhiDecay1< -TMath::PiOver2()) deltaPhiDecay1+=TMath::TwoPi();
+  if(deltaPhiDecay1>3*TMath::PiOver2()) deltaPhiDecay1-=TMath::TwoPi();
+  
+  Float_t deltaPhiDecay2 = mom2.Phi()-phiAssoc;
+  if(deltaPhiDecay2< -TMath::PiOver2()) deltaPhiDecay2+=TMath::TwoPi();
+  if(deltaPhiDecay2>3*TMath::PiOver2()) deltaPhiDecay2-=TMath::TwoPi();
+  
+  Float_t xEDecay1   =-zTDecay1*TMath::Cos(deltaPhiDecay1); // -(px*pxTrig+py*pyTrig)/(ptTrig*ptTrig);
+  Float_t xEDecay2   =-zTDecay2*TMath::Cos(deltaPhiDecay2); // -(px*pxTrig+py*pyTrig)/(ptTrig*ptTrig);
+  
+  if(bChargedOrNeutral) // correlate with charges
+  {
+    fhDeltaPhiDecayCharged->Fill(ptDecay1, deltaPhiDecay1);
+    fhDeltaPhiDecayCharged->Fill(ptDecay2, deltaPhiDecay2);
+    
+    if(GetDebug() > 1) printf("AliAnaParticleHadronCorrelation::FillDecayPhotonHistograms( Charged corr) - deltaPhoton1 = %f, deltaPhoton2 = %f \n", deltaPhiDecay1, deltaPhiDecay2);
+    
+    if( (deltaPhiDecay1 > fDeltaPhiMinCut) && ( deltaPhiDecay1 < fDeltaPhiMaxCut) )
+    {
+      fhZTDecayCharged->Fill(ptDecay1,zTDecay1); 
+      fhXEDecayCharged->Fill(ptDecay1,xEDecay1); 
+    }
+    if( (deltaPhiDecay2 > fDeltaPhiMinCut) && ( deltaPhiDecay2 < fDeltaPhiMaxCut) )
+    {
+      fhZTDecayCharged->Fill(ptDecay2,zTDecay2);
+      fhXEDecayCharged->Fill(ptDecay2,xEDecay2);
+    }
+  }
+  else // correlate with neutrals
+  {
+    fhDeltaPhiDecayCharged->Fill(ptDecay1, deltaPhiDecay1);
+    fhDeltaPhiDecayCharged->Fill(ptDecay2, deltaPhiDecay2);
+    
+    if(GetDebug() > 1) printf("AliAnaParticleHadronCorrelation::FillDecayPhotonHistograms(Neutral corr) - deltaPhoton1 = %f, deltaPhoton2 = %f \n", deltaPhiDecay1, deltaPhiDecay2);
+    
+    if( (deltaPhiDecay1 > fDeltaPhiMinCut) && ( deltaPhiDecay1 < fDeltaPhiMaxCut) )
+    {
+      fhZTDecayCharged->Fill(ptDecay1,zTDecay1); 
+      fhXEDecayCharged->Fill(ptDecay1,xEDecay1); 
+    }
+    if( (deltaPhiDecay2 > fDeltaPhiMinCut) && ( deltaPhiDecay2 < fDeltaPhiMaxCut) )
+    {
+      fhZTDecayCharged->Fill(ptDecay2,zTDecay2);
+      fhXEDecayCharged->Fill(ptDecay2,xEDecay2);
+    }    
+  }
+}
+
+//______________________________________________________________________________________________________________________________________________________
+void AliAnaParticleHadronCorrelation::FillNeutralAngularCorrelationHistograms(const Float_t ptAssoc,  const Float_t ptTrig,  
+                                                                              const Float_t phiAssoc, const Float_t phiTrig,  Float_t &     deltaPhi,
+                                                                              const Float_t etaAssoc, const Float_t etaTrig)
+{
+  // Fill angular correlation related histograms
+  
+  Float_t deltaEta    = etaTrig-etaAssoc;
+  deltaPhi    = phiTrig-phiAssoc;
+  
+  if(deltaPhi <= -TMath::PiOver2()) deltaPhi+=TMath::TwoPi();
+  if(deltaPhi > 3*TMath::PiOver2()) deltaPhi-=TMath::TwoPi();
+  
+  fhEtaNeutral     ->Fill(ptAssoc,etaAssoc);
+  fhPhiNeutral     ->Fill(ptAssoc,phiAssoc);
+  fhDeltaEtaNeutral->Fill(ptTrig ,deltaEta);
+  fhDeltaPhiNeutral->Fill(ptTrig ,deltaPhi);
+  
+  if(ptAssoc > 2 ) fhDeltaPhiDeltaEtaNeutral->Fill(deltaPhi, deltaEta);
+  
+}
+
+//_____________________________________________________________________________________________________________________________
+void AliAnaParticleHadronCorrelation::FillNeutralUnderlyingEventSidesHistograms(const Float_t ptTrig,   const Float_t ptAssoc, 
+                                                                                const Float_t xE,       const Float_t hbpXE, 
+                                                                                const Float_t zT,       const Float_t hbpZT, 
+                                                                                const Float_t deltaPhi)
+{
+  // Fill underlying event histograms to the left and right of trigger
+  
+  if((deltaPhi<-fUeDeltaPhiMinCut) && (deltaPhi >-fUeDeltaPhiMaxCut))
+  {  
+    fhDeltaPhiUeLeftNeutral->Fill(ptAssoc, deltaPhi);
+    fhXEUeLeftNeutral      ->Fill(ptTrig , xE);
+    fhPtHbpXEUeLeftNeutral ->Fill(ptTrig , hbpXE);
+    fhZTUeLeftNeutral      ->Fill(ptTrig , zT);
+    fhPtHbpZTUeLeftNeutral ->Fill(ptTrig , hbpZT);
+  }
+  
+  if((deltaPhi > fUeDeltaPhiMinCut) && (deltaPhi < fUeDeltaPhiMaxCut))
+  {  
+    fhDeltaPhiUeRightNeutral->Fill(ptAssoc, deltaPhi);
+    fhXEUeRightNeutral      ->Fill(ptTrig , xE);
+    fhPtHbpXEUeRightNeutral ->Fill(ptTrig , hbpXE);
+    fhZTUeRightNeutral      ->Fill(ptTrig , zT);
+    fhPtHbpZTUeRightNeutral ->Fill(ptTrig , hbpZT);
+  }
+} 
+
+//_____________________________________________________________
+void AliAnaParticleHadronCorrelation::FillChargedEventMixPool()
+{
+  // Mixed event init
+    
+  //printf("FillChargedEventMixPool for %s\n",GetInputAODName().Data());
+   
+  Int_t nTracks = GetCTSTracks()->GetEntriesFast();
+  
+  fhNtracksAll->Fill(nTracks);
+  
+  AliAnalysisManager   * manager      = AliAnalysisManager::GetAnalysisManager();
+  AliInputEventHandler * inputHandler = dynamic_cast<AliInputEventHandler*>(manager->GetInputEventHandler());
+  
+  if(!inputHandler) return ;
+  
+  if(inputHandler->IsEventSelected( ) & GetReader()->GetEventTriggerMask())
+  {
+    fhNtracksTrigger->Fill(nTracks);
+  }
+  
+  if(inputHandler->IsEventSelected( ) & AliVEvent::kMB)
+  {
+    
+    fhNtracksMB->Fill(nTracks);
+        
+    if(fUseMixStoredInReader && GetReader()->GetLastTracksMixedEvent() == GetEventNumber())
+    {
+      //printf("%s : Pool already filled for this event !!!\n",GetInputAODName().Data());
+      return ; // pool filled previously for another trigger
+    }
+    Int_t eventBin = GetEventMixBin();
+    
+    //Check that the bin exists, if not (bad determination of RP, centrality or vz bin) do nothing
+    if(eventBin < 0) return;
+    
+    TObjArray * mixEventTracks = new TObjArray;
+    
+    if(fUseMixStoredInReader) 
+    {
+      fListMixEvents[eventBin] = GetReader()->GetListWithMixedEventsForTracks(eventBin);
+    }
+    
+    if(!fListMixEvents[eventBin]) fListMixEvents[eventBin] = new TList();
+    
+    //printf("%s ***** Pool Event bin : %d - nTracks %d\n",GetInputAODName().Data(),eventBin, GetCTSTracks()->GetEntriesFast());
+    
+    TList * pool = fListMixEvents[eventBin];
+    
+    TVector3 p3;  
+    for(Int_t ipr = 0;ipr < GetCTSTracks()->GetEntriesFast() ; ipr ++ )
+    {
+      AliVTrack * track = (AliVTrack *) (GetCTSTracks()->At(ipr)) ;
+      
+      Double_t mom[3] = {track->Px(),track->Py(),track->Pz()};
+      p3.SetXYZ(mom[0],mom[1],mom[2]);
+      Float_t pt   = p3.Pt();
+      
+      //Select only hadrons in pt range
+      if(pt < fMinAssocPt || pt > fMaxAssocPt) continue ;
+      
+      AliAODPWG4Particle * mixedTrack = new AliAODPWG4Particle(mom[0],mom[1],mom[2],0);
+      mixedTrack->SetDetector("CTS");
+      mixedTrack->SetChargedBit(track->Charge()>0);
+      
+      mixEventTracks->Add(mixedTrack);
+    }
+    
+    //Set the event number where the last event was added, to avoid double pool filling
+    GetReader()->SetLastTracksMixedEvent(GetEventNumber());
+    
+    pool->AddFirst(mixEventTracks);
+    mixEventTracks = 0;
+    //printf("Pool size %d, max %d\n",pool->GetSize(), GetNMaxEvMix());
+    
+    if(pool->GetSize() > GetNMaxEvMix())
+    {//Remove last event
+      TClonesArray * tmp = static_cast<TClonesArray*>(pool->Last()) ;
+      pool->RemoveLast() ;
+      delete tmp ;
+    }
+        
+  } // MB event
+  
+}
+
 //____________________________________________________________
 TObjString* AliAnaParticleHadronCorrelation::GetAnalysisCuts()
 {
@@ -150,7 +676,7 @@ TObjString* AliAnaParticleHadronCorrelation::GetAnalysisCuts()
   parList+=onePar ;
   snprintf(onePar,buffersize,"Name of AOD Pi0 Branch %s ", fPi0AODBranchName.Data());
   parList+=onePar ;
-  snprintf(onePar,buffersize,"Do Decay-hadron correlation ?  %d", fPi0Trigger) ;
+  snprintf(onePar,buffersize,"Do Decay-hadron correlation ?  pi0 %d, decay %d", fPi0Trigger, fDecayTrigger) ;
   parList+=onePar ;
   snprintf(onePar,buffersize,"Select absolute leading for cluster triggers ? %d or Near Side Leading %d \n", 
            fMakeAbsoluteLeading, fMakeNearSideLeading) ;
@@ -165,9 +691,6 @@ TObjString* AliAnaParticleHadronCorrelation::GetAnalysisCuts()
   //Get parameters set in base class.
   parList += GetBaseParametersList() ;
   
-  //Get parameters set in PID class.
-  //parList += GetCaloPID()->GetPIDParametersList() ;
-  
   //Get parameters set in FiducialCut class (not available yet)
   //parlist += GetFidCut()->GetFidCutParametersList() 
   
@@ -181,12 +704,13 @@ TList *  AliAnaParticleHadronCorrelation::GetCreateOutputObjects()
   
   // Create histograms to be saved in output file and 
   // store them in fOutputContainer
+  
   TList * outputContainer = new TList() ; 
   outputContainer->SetName("CorrelationHistos") ; 
   
-  Int_t   nptbins = GetHistogramRanges()->GetHistoPtBins(); Int_t  nphibins = GetHistogramRanges()->GetHistoPhiBins(); Int_t   netabins = GetHistogramRanges()->GetHistoEtaBins();
-  Float_t ptmax   = GetHistogramRanges()->GetHistoPtMax();  Float_t phimax  = GetHistogramRanges()->GetHistoPhiMax();  Float_t etamax   = GetHistogramRanges()->GetHistoEtaMax();
-  Float_t ptmin   = GetHistogramRanges()->GetHistoPtMin();  Float_t phimin  = GetHistogramRanges()->GetHistoPhiMin();  Float_t etamin   = GetHistogramRanges()->GetHistoEtaMin();      
+  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();
+  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();
+  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();    
   
   fhPtLeading  = new TH1F ("hPtLeading","p_T distribution of leading particles", nptbins,ptmin,ptmax); 
   fhPtLeading->SetXTitle("p_{T}^{trig} (GeV/c)");
@@ -201,12 +725,27 @@ TList *  AliAnaParticleHadronCorrelation::GetCreateOutputObjects()
   outputContainer->Add(fhPhiLeading);
   outputContainer->Add(fhEtaLeading);
   
+  fhPtLeadingCentrality   = new TH2F("hPtLeadingCentrality","Leading particle p_{T} vs centrality",nptbins,ptmin,ptmax,100,0.,100) ;
+  fhPtLeadingCentrality->SetXTitle("p_{T}^{trig} (GeV/c)");
+  fhPtLeadingCentrality->SetYTitle("Centrality (%)");
+  outputContainer->Add(fhPtLeadingCentrality) ;  
+  
+  fhPtLeadingEventPlane  = new TH2F("hPtLeadingEventPlane","Leading particle p_{T} vs event plane angle",nptbins,ptmin,ptmax, 100,0.,TMath::Pi()) ;
+  fhPtLeadingEventPlane->SetXTitle("p_{T}^{trig} (GeV/c)");
+  fhPtLeadingEventPlane->SetXTitle("EP angle (rad)");
+  outputContainer->Add(fhPtLeadingEventPlane) ;
+  
+  fhLeadingEventPlaneCentrality  = new TH2F("hLeadingEventPlane","Leading particle centrality vs event plane angle",100,0.,100,100,0.,TMath::Pi()) ;
+  fhLeadingEventPlaneCentrality->SetXTitle("Centrality (%)");
+  fhLeadingEventPlaneCentrality->SetYTitle("EP angle (rad)");
+  outputContainer->Add(fhLeadingEventPlaneCentrality) ;
+  
   //Correlation with charged hadrons
   if(GetReader()->IsCTSSwitchedOn()) 
   {
     fhDeltaPhiDeltaEtaCharged  = new TH2F
-    ("hDeltaPhiDeltaEtaCharged","#phi_{trigger} - #phi_{h^{#pm}} vs #eta_{trigger} - #eta_{h^{#pm}}",
-     140,-2.,5.,200,-2,2); 
+    ("hDeltaPhiDeltaEtaCharged","#eta_{trigger} - #eta_{h^{#pm}} vs #phi_{trigger} - #phi_{h^{#pm}}",
+    ndeltaphibins ,deltaphimin,deltaphimax,ndeltaetabins,deltaetamin,deltaetamax); 
     fhDeltaPhiDeltaEtaCharged->SetXTitle("#Delta #phi");
     fhDeltaPhiDeltaEtaCharged->SetYTitle("#Delta #eta");    
     
@@ -224,25 +763,31 @@ TList *  AliAnaParticleHadronCorrelation::GetCreateOutputObjects()
     
     fhDeltaPhiCharged  = new TH2F
     ("hDeltaPhiCharged","#phi_{trigger} - #phi_{h^{#pm}} vs p_{T trigger}",
-     nptbins,ptmin,ptmax,140,-2.,5.); 
+     nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
     fhDeltaPhiCharged->SetYTitle("#Delta #phi");
     fhDeltaPhiCharged->SetXTitle("p_{T trigger} (GeV/c)");
     
     fhDeltaPhiChargedPt  = new TH2F
     ("hDeltaPhiChargedPt","#phi_{trigger} - #phi_{#h^{#pm}} vs p_{T h^{#pm}}",
-     nptbins,ptmin,ptmax,140,-2.,5.);
+     nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
     fhDeltaPhiChargedPt->SetYTitle("#Delta #phi");
     fhDeltaPhiChargedPt->SetXTitle("p_{T h^{#pm}} (GeV/c)");
     
     fhDeltaPhiUeChargedPt  = new TH2F
     ("hDeltaPhiUeChargedPt","#phi_{trigger} - #phi_{#Ueh^{#pm}} vs p_{T Ueh^{#pm}}",
-     nptbins,ptmin,ptmax,140,-2.,5.);
+     nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
     fhDeltaPhiUeChargedPt->SetYTitle("#Delta #phi");
     fhDeltaPhiUeChargedPt->SetXTitle("p_{T h^{#pm}} (GeV/c)");
     
+    fhUePart  =  new TH1F("hUePart","UE particles distribution vs pt trig",
+             nptbins,ptmin,ptmax); 
+    fhUePart->SetYTitle("dNch");
+    fhUePart->SetXTitle("p_{T trigger}");
+    
+    
     fhDeltaEtaCharged  = new TH2F
     ("hDeltaEtaCharged","#eta_{trigger} - #eta_{h^{#pm}} vs p_{T trigger}",
-     nptbins,ptmin,ptmax,200,-2,2); 
+     nptbins,ptmin,ptmax,ndeltaetabins,deltaetamin,deltaetamax);  
     fhDeltaEtaCharged->SetYTitle("#Delta #eta");
     fhDeltaEtaCharged->SetXTitle("p_{T trigger} (GeV/c)");
     
@@ -337,14 +882,15 @@ TList *  AliAnaParticleHadronCorrelation::GetCreateOutputObjects()
     outputContainer->Add(fhDeltaEtaCharged) ;
     outputContainer->Add(fhDeltaPhiChargedPt) ;
     outputContainer->Add(fhDeltaPhiUeChargedPt) ;
-    
+    outputContainer->Add(fhUePart);
+
     outputContainer->Add(fhXECharged) ;
     outputContainer->Add(fhXEPosCharged) ;
     outputContainer->Add(fhXENegCharged) ;
     outputContainer->Add(fhXEUeCharged) ;
     outputContainer->Add(fhPtHbpXECharged) ;
     outputContainer->Add(fhPtHbpXEUeCharged) ;
-    
+
     outputContainer->Add(fhZTCharged) ;
     outputContainer->Add(fhZTPosCharged) ;
     outputContainer->Add(fhZTNegCharged) ;
@@ -368,12 +914,12 @@ TList *  AliAnaParticleHadronCorrelation::GetCreateOutputObjects()
       for(Int_t im=0; im<nMultiBins; im++)
       {
         fhTrigDeltaPhiCharged[im]  = new TH2F 
-        (Form("hTrigDeltaPhiCharged_%d",im),Form("hTrigDeltaPhiCharged_%d",im), nptbins,ptmin,ptmax, 140,-2.,5.); 
+        (Form("hTrigDeltaPhiCharged_%d",im),Form("hTrigDeltaPhiCharged_%d",im), nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax); 
         fhTrigDeltaPhiCharged[im]->SetXTitle("p_{T trigger} (GeV/c)");
         fhTrigDeltaPhiCharged[im]->SetYTitle("#Delta #phi");
         
         fhTrigDeltaEtaCharged[im]  = new TH2F 
-        (Form("hTrigDeltaEtaCharged_%d",im),Form("hTrigDeltaEtaCharged_%d",im), nptbins,ptmin,ptmax, 200,-2,2); 
+        (Form("hTrigDeltaEtaCharged_%d",im),Form("hTrigDeltaEtaCharged_%d",im), nptbins,ptmin,ptmax, ndeltaetabins ,deltaetamin,deltaetamax); 
         fhTrigDeltaEtaCharged[im]->SetXTitle("p_{T trigger} (GeV/c)");
         fhTrigDeltaEtaCharged[im]->SetYTitle("#Delta #eta");
         
@@ -406,69 +952,59 @@ TList *  AliAnaParticleHadronCorrelation::GetCreateOutputObjects()
       }
     }
     
-    fhAssocPt           = new TH2F("hAssocPt", " Trigger p_{T} vs associated hadron p_{T}",
-                                   nptbins, ptmin, ptmax,nptbins,ptmin,ptmax);
-    fhAssocPt->SetXTitle("p_{T trigger}");
-    fhAssocPt->SetYTitle("p_{T associated}");
-    outputContainer->Add(fhAssocPt) ;
-    
-    fhAssocPtBkg        = new TH2F("hAssocPtBkg", " Trigger p_{T} vs associated hadron p_{T} from background",
-                                   nptbins, ptmin, ptmax,nptbins,ptmin,ptmax);
-    fhAssocPtBkg->SetXTitle("p_{T trigger}");
-    fhAssocPtBkg->SetYTitle("p_{T associated}");
-    outputContainer->Add(fhAssocPtBkg) ;
-    
-    fhDeltaPhiBrad = new TH2F("hDeltaPhiBrad","atan2(sin(#Delta #phi), cos(#Delta #phi))/#pi vs p_{T trigger} ", 
-                              nptbins, ptmin, ptmax,288, -1.0/3.0, 5.0/3.0);
-    fhDeltaPhiBrad->SetXTitle("p_{T trigger}");
-    fhDeltaPhiBrad->SetYTitle("atan2(sin(#Delta #phi), cos(#Delta #phi))/#pi");
-    outputContainer->Add(fhDeltaPhiBrad) ;
-
-    fhXE       = new TH2F("hXE", "x_{E} vs p_{T trigger}", 
-                          nptbins, ptmin, ptmax,200, 0.0, 2.0);
-    fhXE->SetXTitle("p_{T trigger}");
-    fhXE->SetYTitle("x_{E}");
-    outputContainer->Add(fhXE);
-    
-    fhZT       = new TH2F("hZT", "z_{T} vs p_{T trigger}", 
-                          nptbins, ptmin, ptmax,200, 0.0, 2.0);
-    fhZT->SetXTitle("p_{T trigger}");
-    fhZT->SetYTitle("z_{T}");
-    outputContainer->Add(fhZT);
-    
-    fhDeltaPhiAssocPtBin     = new TH2F*[fNAssocPtBins] ;
-    fhDeltaPhiAssocPtBinHMPID= new TH2F*[fNAssocPtBins] ;
-    fhDeltaPhiAssocPtBinHMPIDAcc= new TH2F*[fNAssocPtBins] ;
-    fhDeltaPhiBradAssocPtBin = new TH2F*[fNAssocPtBins] ;
-    fhXEAssocPtBin           = new TH2F*[fNAssocPtBins] ;
-    fhZTAssocPtBin           = new TH2F*[fNAssocPtBins] ;
+    if(fFillBradHisto)
+    {
+      fhAssocPtBkg        = new TH2F("hAssocPtBkg", " Trigger p_{T} vs associated hadron p_{T} from background",
+                                     nptbins, ptmin, ptmax,nptbins,ptmin,ptmax);
+      fhAssocPtBkg->SetXTitle("p_{T trigger}");
+      fhAssocPtBkg->SetYTitle("p_{T associated}");
+      outputContainer->Add(fhAssocPtBkg) ;
+      
+      fhDeltaPhiBrad = new TH2F("hDeltaPhiBrad","atan2(sin(#Delta #phi), cos(#Delta #phi))/#pi vs p_{T trigger} ", 
+                                nptbins, ptmin, ptmax,288, -1.0/3.0, 5.0/3.0);
+      fhDeltaPhiBrad->SetXTitle("p_{T trigger}");
+      fhDeltaPhiBrad->SetYTitle("atan2(sin(#Delta #phi), cos(#Delta #phi))/#pi");
+      outputContainer->Add(fhDeltaPhiBrad) ;
+    }
+    
+    fhDeltaPhiAssocPtBin       = new TH2F*[fNAssocPtBins] ;
+    fhDeltaPhiAssocPtBinDEta08 = new TH2F*[fNAssocPtBins] ;
+    fhXEAssocPtBin             = new TH2F*[fNAssocPtBins] ;
+    fhZTAssocPtBin             = new TH2F*[fNAssocPtBins] ;
+    
+    if(fFillBradHisto)  
+      fhDeltaPhiBradAssocPtBin = new TH2F*[fNAssocPtBins] ;
+    
+    if(fPi0Trigger || fDecayTrigger)
+    {
+      fhDeltaPhiAssocPtBin       = new TH2F*[fNAssocPtBins] ;
+      fhDeltaPhiAssocPtBinDEta08 = new TH2F*[fNAssocPtBins] ;
+      fhXEAssocPtBin             = new TH2F*[fNAssocPtBins] ;
+      fhZTAssocPtBin             = new TH2F*[fNAssocPtBins] ;
+      fhXEDecayChargedAssocPtBin = new TH2F*[fNAssocPtBins] ;
+      fhZTDecayChargedAssocPtBin = new TH2F*[fNAssocPtBins] ;
+      fhDeltaPhiDecayChargedAssocPtBin = new TH2F*[fNAssocPtBins] ;
+    }
+
+    if(fHMPIDCorrelation)
+    {
+      fhDeltaPhiAssocPtBinHMPID   = new TH2F*[fNAssocPtBins] ;
+      fhDeltaPhiAssocPtBinHMPIDAcc= new TH2F*[fNAssocPtBins] ;
+    }
     
     for(Int_t i = 0 ; i < fNAssocPtBins ; i++)
     {
       fhDeltaPhiAssocPtBin[i] = new TH2F(Form("hDeltaPhiPtAssocPt%2.1f_%2.1f", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]), 
                                          Form("#Delta #phi vs p_{T trigger} for associated p_{T} bin [%2.1f,%2.1f]", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]), 
-                                         nptbins, ptmin, ptmax,140,-2.,5.);
+                                         nptbins, ptmin, ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
       fhDeltaPhiAssocPtBin[i]->SetXTitle("p_{T trigger}");
       fhDeltaPhiAssocPtBin[i]->SetYTitle("#Delta #phi");
-      fhDeltaPhiAssocPtBinHMPID[i] = new TH2F(Form("hDeltaPhiPtAssocPt%2.1f_%2.1fHMPID", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]), 
-                                         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]), 
-                                         nptbins, ptmin, ptmax,140,-2.,5.);
-      fhDeltaPhiAssocPtBinHMPID[i]->SetXTitle("p_{T trigger}");
-      fhDeltaPhiAssocPtBinHMPID[i]->SetYTitle("#Delta #phi");      
-      
-      fhDeltaPhiAssocPtBinHMPIDAcc[i] = new TH2F(Form("hDeltaPhiPtAssocPt%2.1f_%2.1fHMPIDAcc", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]), 
-                                              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]), 
-                                              nptbins, ptmin, ptmax,140,-2.,5.);
-      fhDeltaPhiAssocPtBinHMPIDAcc[i]->SetXTitle("p_{T trigger}");
-      fhDeltaPhiAssocPtBinHMPIDAcc[i]->SetYTitle("#Delta #phi");    
-      
-      fhDeltaPhiBradAssocPtBin[i] = new TH2F(Form("hDeltaPhiBradPtAssocPt%2.1f_%2.1f", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]), 
-                                             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]), 
-                                             nptbins, ptmin, ptmax,288, -1.0/3.0, 5.0/3.0);
-      fhDeltaPhiBradAssocPtBin[i]->SetXTitle("p_{T trigger}");
-      fhDeltaPhiBradAssocPtBin[i]->SetYTitle("atan2(sin(#Delta #phi), cos(#Delta #phi))/#pi");
-      
+       
+      fhDeltaPhiAssocPtBinDEta08[i] = new TH2F(Form("hDeltaPhiDeltaEta0.8PtAssocPt%2.1f_%2.1f", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]), 
+                                               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]), 
+                                               nptbins, ptmin, ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
+      fhDeltaPhiAssocPtBinDEta08[i]->SetXTitle("p_{T trigger}");
+      fhDeltaPhiAssocPtBinDEta08[i]->SetYTitle("#Delta #phi");      
       
       fhXEAssocPtBin[i]       = new TH2F(Form("hXEAssocPtBin%1.f_%1.f", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]), 
                                          Form("x_{E} vs p_{T trigger} for associated p_{T} bin [%2.1f,%2.1f]", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]), 
@@ -483,24 +1019,81 @@ TList *  AliAnaParticleHadronCorrelation::GetCreateOutputObjects()
       fhZTAssocPtBin[i]->SetYTitle("z_{T}");
       
       outputContainer->Add(fhDeltaPhiAssocPtBin[i]) ;
-      outputContainer->Add(fhDeltaPhiAssocPtBinHMPID[i]) ;
-      outputContainer->Add(fhDeltaPhiAssocPtBinHMPIDAcc[i]) ;
-      outputContainer->Add(fhDeltaPhiBradAssocPtBin[i]) ;
+      outputContainer->Add(fhDeltaPhiAssocPtBinDEta08[i]) ;
       outputContainer->Add(fhXEAssocPtBin[i]);
       outputContainer->Add(fhZTAssocPtBin[i]);
+      
+      if(fPi0Trigger || fDecayTrigger) 
+      {
+        fhDeltaPhiDecayChargedAssocPtBin[i] = new TH2F(Form("hDeltaPhiPtDecayChargedAssocPt%2.1f_%2.1f", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]), 
+                                           Form("#Delta #phi vs p_{T trigger} tagged as decay for associated p_{T} bin [%2.1f,%2.1f]", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]), 
+                                           nptbins, ptmin, ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
+        fhDeltaPhiDecayChargedAssocPtBin[i]->SetXTitle("p_{T trigger}");
+        fhDeltaPhiDecayChargedAssocPtBin[i]->SetYTitle("#Delta #phi");
+        
+        fhXEDecayChargedAssocPtBin[i]       = new TH2F(Form("hXEDecayChargedAssocPtBin%1.f_%1.f", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]), 
+                                           Form("x_{E} vs p_{T trigger} tagged as decay for associated p_{T} bin [%2.1f,%2.1f]", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]), 
+                                           nptbins, ptmin, ptmax,200, 0.0, 2.0);
+        fhXEDecayChargedAssocPtBin[i]->SetXTitle("p_{T trigger}");
+        fhXEDecayChargedAssocPtBin[i]->SetYTitle("x_{E}");
+        
+        fhZTDecayChargedAssocPtBin[i]       = new TH2F(Form("hZTDecayChargedAssocPtBin%1.f_%1.f", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]), 
+                                           Form("z_{T} vs p_{T trigger} tagged as decay for associated p_{T} bin [%2.1f,%2.1f]", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]), 
+                                           nptbins, ptmin, ptmax,200, 0.0, 2.0);
+        fhZTDecayChargedAssocPtBin[i]->SetXTitle("p_{T trigger}");
+        fhZTDecayChargedAssocPtBin[i]->SetYTitle("z_{T}");
+        
+        outputContainer->Add(fhDeltaPhiDecayChargedAssocPtBin[i]) ;
+        outputContainer->Add(fhXEDecayChargedAssocPtBin[i]);
+        outputContainer->Add(fhZTDecayChargedAssocPtBin[i]);
+        
+      }
+      
+      if(fFillBradHisto) 
+      {
+        fhDeltaPhiBradAssocPtBin[i] = new TH2F(Form("hDeltaPhiBradPtAssocPt%2.1f_%2.1f", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]), 
+                                               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]), 
+                                               nptbins, ptmin, ptmax,288, -1.0/3.0, 5.0/3.0);
+        fhDeltaPhiBradAssocPtBin[i]->SetXTitle("p_{T trigger}");
+        fhDeltaPhiBradAssocPtBin[i]->SetYTitle("atan2(sin(#Delta #phi), cos(#Delta #phi))/#pi");
+        outputContainer->Add(fhDeltaPhiBradAssocPtBin[i]) ;
+      }       
+      
+      if(fHMPIDCorrelation)
+      {
+        fhDeltaPhiAssocPtBinHMPID[i] = new TH2F(Form("hDeltaPhiPtAssocPt%2.1f_%2.1fHMPID", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]), 
+                                                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]), 
+                                                nptbins, ptmin, ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
+        fhDeltaPhiAssocPtBinHMPID[i]->SetXTitle("p_{T trigger}");
+        fhDeltaPhiAssocPtBinHMPID[i]->SetYTitle("#Delta #phi");      
+        
+        fhDeltaPhiAssocPtBinHMPIDAcc[i] = new TH2F(Form("hDeltaPhiPtAssocPt%2.1f_%2.1fHMPIDAcc", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]), 
+                                                   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]), 
+                                                   nptbins, ptmin, ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
+        fhDeltaPhiAssocPtBinHMPIDAcc[i]->SetXTitle("p_{T trigger}");
+        fhDeltaPhiAssocPtBinHMPIDAcc[i]->SetYTitle("#Delta #phi"); 
+        
+        outputContainer->Add(fhDeltaPhiAssocPtBinHMPID[i]) ;
+        outputContainer->Add(fhDeltaPhiAssocPtBinHMPIDAcc[i]) ;
+        
+      }      
     }
     
-    if(fPi0Trigger)
+    if(fPi0Trigger || fDecayTrigger)
     {
-      fhPtPi0DecayRatio  = new TH2F
-      ("hPtPi0DecayRatio","p_T of #pi^{0} and the ratio of pt for two decay", 
-       nptbins,ptmin,ptmax, 100,0.,2.); 
-      fhPtPi0DecayRatio->SetXTitle("p_{T}^{#pi^{0}} (GeV/c)");
-      fhPtPi0DecayRatio->SetYTitle("p_{T}^{Decay}/p_{T}^{#pi^{0}}");
+      if(fPi0Trigger)
+      {
+        fhPtPi0DecayRatio  = new TH2F
+        ("hPtPi0DecayRatio","p_T of #pi^{0} and the ratio of pt for two decay", 
+         nptbins,ptmin,ptmax, 100,0.,2.); 
+        fhPtPi0DecayRatio->SetXTitle("p_{T}^{#pi^{0}} (GeV/c)");
+        fhPtPi0DecayRatio->SetYTitle("p_{T}^{Decay}/p_{T}^{#pi^{0}}");
+        outputContainer->Add(fhPtPi0DecayRatio) ; 
+      }
       
       fhDeltaPhiDecayCharged  = new TH2F
       ("hDeltaPhiDecayCharged","#phi_{Decay} - #phi_{h^{#pm}} vs p_{T Decay}",
-       nptbins,ptmin,ptmax,140,-2.,5.); 
+       nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax); 
       fhDeltaPhiDecayCharged->SetYTitle("#Delta #phi");
       fhDeltaPhiDecayCharged->SetXTitle("p_{T Decay} (GeV/c)");
       
@@ -516,7 +1109,6 @@ TList *  AliAnaParticleHadronCorrelation::GetCreateOutputObjects()
       fhZTDecayCharged->SetYTitle("z_{decay h^{#pm}}");
       fhZTDecayCharged->SetXTitle("p_{T decay}");      
       
-      outputContainer->Add(fhPtPi0DecayRatio) ; 
       outputContainer->Add(fhDeltaPhiDecayCharged) ; 
       outputContainer->Add(fhXEDecayCharged) ;
       outputContainer->Add(fhZTDecayCharged) ;
@@ -526,14 +1118,14 @@ TList *  AliAnaParticleHadronCorrelation::GetCreateOutputObjects()
     { 
       fhDeltaPhiUeLeftCharged  = new TH2F
       ("hDeltaPhiUeLeftChargedPt","#phi_{trigger} - #phi_{#Ueh^{#pm}} vs p_{T Ueh^{#pm}} with UE left side range of trigger particles",
-       nptbins,ptmin,ptmax,140,-2.,5.);
+       nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
       fhDeltaPhiUeLeftCharged->SetYTitle("#Delta #phi");
       fhDeltaPhiUeLeftCharged->SetXTitle("p_{T h^{#pm}} (GeV/c)");
       outputContainer->Add(fhDeltaPhiUeLeftCharged) ;
       
       fhDeltaPhiUeRightCharged  = new TH2F
       ("hDeltaPhiUeRightChargedPt","#phi_{trigger} - #phi_{#Ueh^{#pm}} vs p_{T Ueh^{#pm}} with UE right side range of trigger particles",
-       nptbins,ptmin,ptmax,140,-2.,5.);
+       nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
       fhDeltaPhiUeRightCharged->SetYTitle("#Delta #phi");
       fhDeltaPhiUeRightCharged->SetXTitle("p_{T h^{#pm}} (GeV/c)");
       outputContainer->Add(fhDeltaPhiUeRightCharged) ;
@@ -596,13 +1188,13 @@ TList *  AliAnaParticleHadronCorrelation::GetCreateOutputObjects()
       
     }  
   }  //Correlation with charged hadrons
-  
+
   //Correlation with neutral hadrons
   if(fNeutralCorr)
   {
     fhDeltaPhiDeltaEtaNeutral  = new TH2F
     ("hDeltaPhiDeltaEtaNeutral","#phi_{trigger} - #phi_{h^{0}} vs #eta_{trigger} - #eta_{h^{0}}",
-     140,-2.,5.,200,-2,2); 
+     ndeltaphibins ,deltaphimin,deltaphimax, ndeltaetabins ,deltaetamin,deltaetamax); 
     fhDeltaPhiDeltaEtaNeutral->SetXTitle("#Delta #phi");
     fhDeltaPhiDeltaEtaNeutral->SetYTitle("#Delta #eta");   
          
@@ -626,19 +1218,19 @@ TList *  AliAnaParticleHadronCorrelation::GetCreateOutputObjects()
     
     fhDeltaPhiNeutralPt  = new TH2F
     ("hDeltaPhiNeutralPt","#phi_{trigger} - #phi_{#pi^{0}} vs p_{T #pi^{0}}}",
-     nptbins,ptmin,ptmax,140,-2.,5.); 
+     nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax); 
     fhDeltaPhiNeutralPt->SetYTitle("#Delta #phi");
     fhDeltaPhiNeutralPt->SetXTitle("p_{T h^{0}} (GeV/c)");
     
     fhDeltaPhiUeNeutralPt  = new TH2F
     ("hDeltaPhiUeNeutralPt","#phi_{trigger} - #phi_{#pi^{0}} vs p_{T #pi^{0}}}",
-     nptbins,ptmin,ptmax,140,-2.,5.); 
+     nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax); 
     fhDeltaPhiUeNeutralPt->SetYTitle("#Delta #phi");
     fhDeltaPhiUeNeutralPt->SetXTitle("p_{T h^{0}} (GeV/c)");
     
     fhDeltaEtaNeutral  = new TH2F
     ("hDeltaEtaNeutral","#eta_{trigger} - #eta_{#pi^{0}} vs p_{T trigger}",
-     nptbins,ptmin,ptmax,200,-2,2); 
+     nptbins,ptmin,ptmax, ndeltaetabins ,deltaetamin,deltaetamax);  
     fhDeltaEtaNeutral->SetYTitle("#Delta #eta");
     fhDeltaEtaNeutral->SetXTitle("p_{T trigger} (GeV/c)");
     
@@ -706,11 +1298,11 @@ TList *  AliAnaParticleHadronCorrelation::GetCreateOutputObjects()
     outputContainer->Add(fhPtHbpZTNeutral) ;
     outputContainer->Add(fhPtHbpZTUeNeutral) ;    
     
-    if(fPi0Trigger)
+    if(fPi0Trigger || fDecayTrigger)
     {
       fhDeltaPhiDecayNeutral  = new TH2F
       ("hDeltaPhiDecayNeutral","#phi_{Decay} - #phi_{h^{0}} vs p_{T Decay}",
-       nptbins,ptmin,ptmax,140,-2.,5.); 
+       nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);  
       fhDeltaPhiDecayNeutral->SetYTitle("#Delta #phi");
       fhDeltaPhiDecayNeutral->SetXTitle("p_{T Decay} (GeV/c)");
       
@@ -736,14 +1328,14 @@ TList *  AliAnaParticleHadronCorrelation::GetCreateOutputObjects()
     { 
       fhDeltaPhiUeLeftNeutral  = new TH2F
       ("hDeltaPhiUeLeftNeutralPt","#phi_{trigger} - #phi_{#Ueh^{0}} vs p_{T h^{0}} with neutral UE left side range of trigger particles",
-       nptbins,ptmin,ptmax,140,-2.,5.);
+       nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax); 
       fhDeltaPhiUeLeftNeutral->SetYTitle("#Delta #phi");
       fhDeltaPhiUeLeftNeutral->SetXTitle("p_{T h^{0}} (GeV/c)");
       outputContainer->Add(fhDeltaPhiUeLeftNeutral) ;
       
       fhDeltaPhiUeRightNeutral  = new TH2F
       ("hDeltaPhiUeRightNeutralPt","#phi_{trigger} - #phi_{#Ueh^{0}} vs p_{T Ueh^{0}} with neutral UE right side range of trigger particles",
-       nptbins,ptmin,ptmax,140,-2.,5.);
+       nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax); 
       fhDeltaPhiUeRightNeutral->SetYTitle("#Delta #phi");
       fhDeltaPhiUeRightNeutral->SetXTitle("p_{T h^{0}} (GeV/c)");
       outputContainer->Add(fhDeltaPhiUeRightNeutral) ;
@@ -814,9 +1406,9 @@ TList *  AliAnaParticleHadronCorrelation::GetCreateOutputObjects()
     fh2phiLeadingParticle=new TH2F("h2phiLeadingParticle","#phi resolustion for trigger particles",nptbins,ptmin,ptmax,100,-1,1);
     fh2phiLeadingParticle->GetXaxis()->SetTitle("p_{T gen Leading} (GeV/c)");
     fh2phiLeadingParticle->GetYaxis()->SetTitle("(#phi_{rec}-#phi_{gen})/#phi_{gen}");
-    
-    fhMCLeadingCount=new TH1F("hMCLeadingTriggerCount","MCLeadingTriggerCount",nptbins,ptmin,ptmax);
-    fhMCLeadingCount->SetXTitle("p_{T trig}");
+
+    fhMCPtLeading  = new TH1F ("hMCPtLeading","MC : p_T distribution of leading particles", nptbins,ptmin,ptmax); 
+    fhMCPtLeading->SetXTitle("p_{T}^{trig} (GeV/c)");
     
     fhMCEtaCharged  = new TH2F
     ("hMCEtaCharged","MC #eta_{h^{#pm}}  vs p_{T #pm}",
@@ -844,13 +1436,13 @@ TList *  AliAnaParticleHadronCorrelation::GetCreateOutputObjects()
     
     fhMCDeltaPhiCharged  = new TH2F
     ("hMCDeltaPhiCharged","#phi_{trigger} - #phi_{h^{#pm}} vs p_{T trigger}",
-     nptbins,ptmin,ptmax,140,-2.,5.); 
+     nptbins,ptmin,ptmax,ndeltaphibins ,deltaphimin,deltaphimax); 
     fhMCDeltaPhiCharged->SetYTitle("#Delta #phi");
     fhMCDeltaPhiCharged->SetXTitle("p_{T trigger} (GeV/c)");
     
     fhMCDeltaPhiChargedPt  = new TH2F
     ("hMCDeltaPhiChargedPt","MC #phi_{trigger} - #phi_{#h^{#pm}} vs p_{T h^{#pm}}",
-     nptbins,ptmin,ptmax,140,-2.,5.);
+     nptbins,ptmin,ptmax,ndeltaphibins ,deltaphimin,deltaphimax); 
     fhMCDeltaPhiChargedPt->SetYTitle("#Delta #phi");
     fhMCDeltaPhiChargedPt->SetXTitle("p_{T h^{#pm}} (GeV/c)");
     
@@ -858,13 +1450,31 @@ TList *  AliAnaParticleHadronCorrelation::GetCreateOutputObjects()
     new TH2F("hMCPtXECharged","x_{E}",
              nptbins,ptmin,ptmax,200,0.,2.); 
     fhMCPtXECharged->SetYTitle("x_{E}");
-    fhMCPtXECharged->SetXTitle("p_{T trigger}");  
+    fhMCPtXECharged->SetXTitle("p_{T trigger}");
+
+    fhMCPtXEUeCharged  = 
+    new TH2F("hMCPtXEUeCharged","x_{E}",
+             nptbins,ptmin,ptmax,200,0.,2.); 
+    fhMCPtXEUeCharged->SetYTitle("x_{E}");
+    fhMCPtXEUeCharged->SetXTitle("p_{T trigger}");
     
     fhMCPtHbpXECharged  = 
     new TH2F("hMCHbpXECharged","MC #xi = ln(1/x_{E}) with charged hadrons",
              nptbins,ptmin,ptmax,200,0.,10.); 
     fhMCPtHbpXECharged->SetYTitle("ln(1/x_{E})");
     fhMCPtHbpXECharged->SetXTitle("p_{T trigger}");
+
+    fhMCPtHbpXEUeCharged =
+    new TH2F("hMCPtHbpXEUeCharged","#xi = ln(1/x_{E}) with charged hadrons,Underlying Event",
+             nptbins,ptmin,ptmax,200,0.,10.); 
+    fhMCPtHbpXEUeCharged->SetYTitle("ln(1/x_{E})");
+    fhMCPtHbpXEUeCharged->SetXTitle("p_{T trigger}");
+
+    fhMCUePart  = 
+    new TH1F("hMCUePart","MC UE particles distribution vs pt trig",
+             nptbins,ptmin,ptmax); 
+    fhMCUePart->SetYTitle("dNch");
+    fhMCUePart->SetXTitle("p_{T trigger}");
     
     fhMCPtZTCharged  = 
     new TH2F("hMCPtZTCharged","z_{T}",
@@ -886,12 +1496,12 @@ TList *  AliAnaParticleHadronCorrelation::GetCreateOutputObjects()
     
     fhMCPtAssocDeltaPhi  = 
     new TH2F("hMCPtAssocDeltaPhi","AOD MC delta phi with associated charged hadrons",
-             nptbins,ptmin,ptmax,140,-2.,5.); 
+             nptbins,ptmin,ptmax,ndeltaphibins ,deltaphimin,deltaphimax); 
     fhMCPtAssocDeltaPhi->SetYTitle("#Delta #phi");
     fhMCPtAssocDeltaPhi->SetXTitle("p_{T trigger} (GeV/c)"); 
         
     outputContainer->Add(fh2phiLeadingParticle);
-    outputContainer->Add(fhMCLeadingCount);
+    outputContainer->Add(fhMCPtLeading);
     outputContainer->Add(fhMCDeltaPhiDeltaEtaCharged);
     outputContainer->Add(fhMCPhiCharged) ;
     outputContainer->Add(fhMCEtaCharged) ;
@@ -900,20 +1510,163 @@ TList *  AliAnaParticleHadronCorrelation::GetCreateOutputObjects()
     
     outputContainer->Add(fhMCDeltaPhiChargedPt) ;
     outputContainer->Add(fhMCPtXECharged) ;
+    outputContainer->Add(fhMCPtXEUeCharged) ;
     outputContainer->Add(fhMCPtZTCharged) ;
     outputContainer->Add(fhMCPtHbpXECharged) ;
+    outputContainer->Add(fhMCPtHbpXEUeCharged);
+    outputContainer->Add(fhMCUePart);
     outputContainer->Add(fhMCPtHbpZTCharged) ;
     outputContainer->Add(fhMCPtTrigPout) ;
     outputContainer->Add(fhMCPtAssocDeltaPhi) ;      
   } //for MC histogram
   
-  return outputContainer;
-  
-}
-
-//____________________________________________________
-void AliAnaParticleHadronCorrelation::InitParameters()
-{
+  if(DoOwnMix())
+  {
+    //create event containers
+    
+    if(!fUseMixStoredInReader || (fUseMixStoredInReader && !GetReader()->ListWithMixedEventsForTracksExists())) 
+    {
+      fListMixEvents= new TList*[GetNCentrBin()*GetNZvertBin()*GetNRPBin()] ;
+      
+      for(Int_t ic=0; ic<GetNCentrBin(); ic++)
+      {
+        for(Int_t iz=0; iz<GetNZvertBin(); iz++)
+        {
+          for(Int_t irp=0; irp<GetNRPBin(); irp++)
+          {
+            //printf("GetCreateOutputObjects - Bins : cent %d, vz %d, RP %d, event %d/%d\n",
+            //       ic,iz, irp, ic*GetNZvertBin()*GetNRPBin()+iz*GetNRPBin()+irp,GetNZvertBin()*GetNRPBin()*GetNCentrBin());
+            fListMixEvents[ic*GetNZvertBin()*GetNRPBin()+iz*GetNRPBin()+irp] = new TList() ;
+            fListMixEvents[ic*GetNZvertBin()*GetNRPBin()+iz*GetNRPBin()+irp]->SetOwner(kFALSE);
+          }
+        }
+      }    
+    }
+    
+    //Init the list in the reader if not done previously
+    if(!GetReader()->ListWithMixedEventsForTracksExists() && fUseMixStoredInReader) 
+    {
+      //printf("%s : Set the list of events \n",GetInputAODName().Data());
+      GetReader()->SetListWithMixedEventsForTracks(fListMixEvents);
+    }
+    
+    fhEventBin=new TH1I("hEventBin","Number of real events per bin(cen,vz,rp)",
+                        GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1,0, 
+                        GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1) ;
+    fhEventBin->SetXTitle("bin");
+    outputContainer->Add(fhEventBin) ;
+    
+    fhEventMixBin=new TH1I("hEventMixBin","Number of events  per bin(cen,vz,rp)",
+                           GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1,0,
+                           GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1) ;
+    fhEventMixBin->SetXTitle("bin");
+    outputContainer->Add(fhEventMixBin) ;
+    
+    fhNtracksAll=new TH1F("hNtracksAll","Number of tracks w/o event trigger",2000,0,2000);
+    outputContainer->Add(fhNtracksAll);
+    
+    fhNtracksTrigger=new TH1F("hNtracksTriggerEvent","Number of tracks w/ event trigger",2000,0,2000);
+    outputContainer->Add(fhNtracksTrigger);
+    
+    fhNtracksMB=new TH1F("hNtracksMBEvent","Number of tracks w/ event trigger kMB",2000,0,2000);
+    outputContainer->Add(fhNtracksMB);
+    
+    fhMixDeltaPhiCharged  = new TH2F
+    ("hMixDeltaPhiCharged","Mixed event : #phi_{trigger} - #phi_{h^{#pm}} vs p_{T trigger}",
+     nptbins,ptmin,ptmax,ndeltaphibins ,deltaphimin,deltaphimax); 
+    fhMixDeltaPhiCharged->SetYTitle("#Delta #phi");
+    fhMixDeltaPhiCharged->SetXTitle("p_{T trigger} (GeV/c)");
+    outputContainer->Add(fhMixDeltaPhiCharged);
+    
+    fhMixDeltaPhiDeltaEtaCharged  = new TH2F
+    ("hMixDeltaPhiDeltaEtaCharged","Mixed event : #phi_{trigger} - #phi_{h^{#pm}} vs #eta_{trigger} - #eta_{h^{#pm}}",
+     ndeltaphibins ,deltaphimin,deltaphimax,ndeltaetabins ,deltaetamin,deltaetamax); 
+    fhMixDeltaPhiDeltaEtaCharged->SetXTitle("#Delta #phi");
+    fhMixDeltaPhiDeltaEtaCharged->SetYTitle("#Delta #eta");
+    outputContainer->Add(fhMixDeltaPhiDeltaEtaCharged);
+    
+    fhMixDeltaPhiChargedAssocPtBin         = new TH2F*[fNAssocPtBins] ;
+    fhMixDeltaPhiChargedAssocPtBinDEta08   = new TH2F*[fNAssocPtBins] ;
+    fhMixDeltaPhiDeltaEtaChargedAssocPtBin = new TH2F*[fNAssocPtBins] ;
+    
+    for(Int_t i = 0 ; i < fNAssocPtBins ; i++)
+    {
+      fhMixDeltaPhiChargedAssocPtBin[i] = new TH2F(Form("hMixDeltaPhiChargedAssocPtBin%2.1f_%2.1f", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]), 
+                                                   Form("Mixed event #Delta #phi vs p_{T trigger} for associated p_{T} bin [%2.1f,%2.1f]", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]), 
+                                                   nptbins, ptmin, ptmax,  ndeltaphibins ,deltaphimin,deltaphimax);
+      fhMixDeltaPhiChargedAssocPtBin[i]->SetXTitle("p_{T trigger}");
+      fhMixDeltaPhiChargedAssocPtBin[i]->SetYTitle("#Delta #phi");
+      
+      fhMixDeltaPhiChargedAssocPtBinDEta08[i] = new TH2F(Form("hMixDeltaPhiDeltaEta0.8ChargedAssocPtBinDEta08%2.1f_%2.1f", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]), 
+                                                         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]), 
+                                                         nptbins, ptmin, ptmax,  ndeltaphibins ,deltaphimin,deltaphimax);
+      fhMixDeltaPhiChargedAssocPtBinDEta08[i]->SetXTitle("p_{T trigger}");
+      fhMixDeltaPhiChargedAssocPtBinDEta08[i]->SetYTitle("#Delta #phi");      
+      
+      fhMixDeltaPhiDeltaEtaChargedAssocPtBin[i] = new TH2F(Form("hMixDeltaPhiDeltaEtaChargedAssocPtBin%2.1f_%2.1f", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]), 
+                                                           Form("Mixed event #Delta #phi vs p_{T trigger} for associated p_{T} bin [%2.1f,%2.1f]", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]), 
+                                                           ndeltaphibins ,deltaphimin,deltaphimax,ndeltaetabins ,deltaetamin,deltaetamax); 
+      fhMixDeltaPhiDeltaEtaChargedAssocPtBin[i]->SetXTitle("#Delta #phi");
+      fhMixDeltaPhiDeltaEtaChargedAssocPtBin[i]->SetYTitle("#Delta #eta");
+      
+      outputContainer->Add(fhMixDeltaPhiChargedAssocPtBin[i]);
+      outputContainer->Add(fhMixDeltaPhiChargedAssocPtBinDEta08[i]);
+      outputContainer->Add(fhMixDeltaPhiDeltaEtaChargedAssocPtBin[i]);
+      
+    }          
+  }
+
+  return outputContainer;
+  
+}
+
+//_________________________________________________________________________________________________
+Bool_t AliAnaParticleHadronCorrelation::GetDecayPhotonMomentum(const AliAODPWG4Particle* trigger, 
+                                                             TLorentzVector & mom1, 
+                                                             TLorentzVector & mom2)
+{
+  // Get the momentum of the pi0/eta assigned decay photons
+  // In case of pi0/eta trigger, we may want to check their decay correlation, 
+  // get their decay children
+  
+  Int_t indexPhoton1 = trigger->GetCaloLabel(0);
+  Int_t indexPhoton2 = trigger->GetCaloLabel(1);
+  Float_t ptTrig     = trigger->Pt();
+  
+  if(indexPhoton1!=-1 || indexPhoton2!=-1) return kFALSE;
+  
+  if(GetDebug() > 1) 
+    printf("AliAnaParticleHadronCorrelation::GetDecayPhotonMomentum() - indexPhoton1 = %d, indexPhoton2 = %d \n", indexPhoton1, indexPhoton2);
+  
+  TObjArray * clusters  = 0x0 ;  
+  if(trigger->GetDetector()=="EMCAL") clusters = GetEMCALClusters() ;
+  else                                clusters = GetPHOSClusters()  ;
+  
+  for(Int_t iclus = 0; iclus < clusters->GetEntriesFast(); iclus++)
+  {
+    AliVCluster * photon =  (AliVCluster*) (clusters->At(iclus));      
+    if(photon->GetID()==indexPhoton1) 
+    {
+      photon->GetMomentum(mom1,GetVertex(0)) ;
+      if(ptTrig) fhPtPi0DecayRatio->Fill(ptTrig, mom1.Pt()/ptTrig);
+    }
+    if(photon->GetID()==indexPhoton2) 
+    {
+      photon->GetMomentum(mom1,GetVertex(0)) ;
+      if(ptTrig > 0) fhPtPi0DecayRatio->Fill(ptTrig, mom2.Pt()/ptTrig);
+    } 
+    
+    if(GetDebug() > 1)printf("AliAnaParticleHadronCorrelation::GetDecayPhotonMomentum() - Photon1 = %f, Photon2 = %f \n", mom1.Pt(), mom2.Pt());
+    
+  } //cluster loop        
+  
+  return kTRUE;
+  
+} 
+
+//____________________________________________________
+void AliAnaParticleHadronCorrelation::InitParameters()
+{
   
   //Initialize the parameters of the analysis.
   SetInputAODName("PWG4Particle");
@@ -927,23 +1680,28 @@ void AliAnaParticleHadronCorrelation::InitParameters()
   fMakeSeveralUE        = kFALSE;
   fUeDeltaPhiMinCut     = 1. ;
   fUeDeltaPhiMaxCut     = 1.5 ;
+  
   fNeutralCorr          = kFALSE ;
   fPi0Trigger           = kFALSE ;
+  fDecayTrigger         = kFALSE ;
+  fHMPIDCorrelation     = kFALSE ;
   
   fMakeAbsoluteLeading  = kTRUE;
   fMakeNearSideLeading  = kFALSE;
 
-  fNAssocPtBins         = 7  ;
-  fAssocPtBinLimit[0]   = 1.5  ; 
-  fAssocPtBinLimit[1]   = 3.  ; 
-  fAssocPtBinLimit[2]   = 5.  ; 
-  fAssocPtBinLimit[3]   = 7.  ; 
-  fAssocPtBinLimit[4]   = 9. ; 
-  fAssocPtBinLimit[5]   = 12. ; 
-  fAssocPtBinLimit[6]   = 15. ;
-  fAssocPtBinLimit[7]   = 20. ;
-  fAssocPtBinLimit[8]   = 100.;
-  fAssocPtBinLimit[9]   = 200.;
+  fNAssocPtBins         = 9   ;
+  fAssocPtBinLimit[0]   = 0.2 ; 
+  fAssocPtBinLimit[1]   = 2.0 ; 
+  fAssocPtBinLimit[2]   = 4.0 ; 
+  fAssocPtBinLimit[3]   = 6.0 ; 
+  fAssocPtBinLimit[4]   = 8.0 ; 
+  fAssocPtBinLimit[5]   = 10. ; 
+  fAssocPtBinLimit[6]   = 12. ;
+  fAssocPtBinLimit[7]   = 15. ;
+  fAssocPtBinLimit[8]   = 25. ;
+  fAssocPtBinLimit[9]   = 50. ;
+  
+  fUseMixStoredInReader = kTRUE;
   
 }
 
@@ -964,7 +1722,8 @@ void  AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD()
     abort();
   }
        
-  if(GetDebug() > 1){
+  if(GetDebug() > 1)
+  {
     printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - Begin hadron correlation analysis, fill AODs \n");
     printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - In particle branch aod entries %d\n", GetInputAODBranch()->GetEntriesFast());
     printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - In CTS aod entries %d\n",   GetCTSTracks()    ->GetEntriesFast());
@@ -977,6 +1736,9 @@ void  AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD()
   GetReader()->GetVertex(v);
   if(!GetMixedEvent() && TMath::Abs(v[2]) > GetZvertexCut()) return ;   
   
+  // Fill the pool with tracks if requested
+  if(DoOwnMix()) FillChargedEventMixPool();
+
   //Loop on stored AOD particles, find leading trigger
   Double_t ptTrig      = fMinTriggerPt ;
   fLeadingTriggerIndex = -1 ;
@@ -1060,7 +1822,7 @@ void  AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms()
       Int_t check = CheckMixedEventVertex(particle->GetCaloLabel(0), particle->GetTrackLabel(0));
       if(check ==  0) continue;
       if(check == -1) return;
-      
+            
       //check if the particle is isolated or if we want to take the isolation into account
       if(OnlyIsolated() && !particle->IsIsolated()) continue;
       
@@ -1084,7 +1846,7 @@ void  AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms()
     {
       Bool_t in = GetFiducialCut()->IsInFiducialCut(*particle->Momentum(),particle->GetDetector()) ;
       if(! in ) return ;
-    }    
+    }
     
     // Check if the particle is isolated or if we want to take the isolation into account
     if(OnlyIsolated() && !particle->IsIsolated()) return;
@@ -1112,11 +1874,23 @@ void  AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms()
     // no problem was found, like not absolute leading, or bad vertex in mixing.
     if(okcharged && okneutral)
     {
-      fhPtLeading->Fill(particle->Pt());
+      Float_t pt = particle->Pt();
+      fhPtLeading->Fill(pt);
+      
       Float_t phi = particle->Phi();
       if(phi<0)phi+=TMath::TwoPi();
-      fhPhiLeading->Fill(particle->Pt(), phi);
-      fhEtaLeading->Fill(particle->Pt(), particle->Eta());
+      fhPhiLeading->Fill(pt, phi);
+      
+      fhEtaLeading->Fill(pt, particle->Eta());
+      //printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms() - Leading particle : pt %f, eta %f, phi %f\n",particle->Pt(),particle->Eta(),phi);
+      
+      Float_t cen = GetEventCentrality();
+      Float_t ep  = GetEventPlaneAngle();
+      
+      fhPtLeadingCentrality        ->Fill(pt,cen);
+      fhPtLeadingEventPlane        ->Fill(pt,ep);
+      fhLeadingEventPlaneCentrality->Fill(cen,ep);
+      
     }//ok charged && neutral
   }//Aod branch loop
   
@@ -1131,13 +1905,34 @@ Bool_t  AliAnaParticleHadronCorrelation::MakeChargedCorrelation(AliAODPWG4Partic
                                                                 const TObjArray* pl, const Bool_t bFillHisto)
 {  
   // Charged Hadron Correlation Analysis
-  if(GetDebug() > 1)printf("AliAnaParticleHadronCorrelation::MakeChargedCorrelation() - Make trigger particle - charged hadron correlation \n");
+  if(GetDebug() > 1) 
+    printf("AliAnaParticleHadronCorrelation::MakeChargedCorrelation() - Make trigger particle - charged hadron correlation \n");
+    
+  Float_t phiTrig = aodParticle->Phi();
+  Float_t etaTrig = aodParticle->Eta();
+  Float_t ptTrig  = aodParticle->Pt();  
+  Bool_t   decay  = aodParticle->IsTagged();
+
+  Float_t pt       = -100. ;
+  Float_t zT       = -100. ; 
+  Float_t xE       = -100. ; 
+  Float_t hbpXE    = -100. ; 
+  Float_t hbpZT    = -100. ; 
+  Float_t phi      = -100. ;
+  Float_t eta      = -100. ;
+  Float_t pout     = -100. ;
+  Float_t deltaPhi = -100. ;
+  
+  TVector3 p3;  
+  TLorentzVector photonMom ;   
+  TObjArray * reftracks = 0x0;
+  Int_t nrefs           = 0;
+  Int_t nTracks         = GetCTSTracks()->GetEntriesFast() ;
   
-  Int_t evtIndex11   = -1 ; //cluster trigger or pi0 trigger 
+  // Mixed event settings
+  Int_t evtIndex11   = -1 ; // cluster trigger or pi0 trigger 
   Int_t evtIndex12   = -1 ; // pi0 trigger
   Int_t evtIndex13   = -1 ; // charged trigger
-  Int_t indexPhoton1 = -1 ;
-  Int_t indexPhoton2 = -1 ;  
   
   Double_t v[3]      = {0,0,0}; //vertex ;
   GetReader()->GetVertex(v);
@@ -1149,82 +1944,17 @@ Bool_t  AliAnaParticleHadronCorrelation::MakeChargedCorrelation(AliAODPWG4Partic
     evtIndex13 = GetMixedEvent()->EventIndex(aodParticle->GetTrackLabel(0)) ;
   }
   
-  Double_t phiTrig = aodParticle->Phi();
-  Double_t ptTrig  = aodParticle->Pt();  
-  
-  Double_t pt             = -100. ;
-  Double_t px             = -100. ;
-  Double_t py             = -100. ;
-  Double_t zT             = -100. ; 
-  Double_t xE             = -100. ; 
-  Double_t hbpXE          = -100. ; 
-  Double_t hbpZT          = -100. ; 
-  Double_t phi            = -100. ;
-  Double_t eta            = -100. ;
-  Double_t pout           = -100. ;
-  
-  Double_t ptDecay1       = 0. ;
-  Double_t pxDecay1       = 0. ;
-  Double_t pyDecay1       = 0. ;
-  Double_t phiDecay1      = 0. ;
-  Double_t ptDecay2       = 0. ;
-  Double_t pxDecay2       = 0. ;
-  Double_t pyDecay2       = 0. ;
-  Double_t phiDecay2      = 0. ;
-  
-  Double_t ratDecay1      = -100. ;  
-  Double_t ratDecay2      = -100. ; 
-  Double_t deltaPhi       = -100. ;
-  Double_t deltaPhiOrg    = -100. ;
-  Double_t deltaPhiDecay1 = -100. ;
-  Double_t deltaPhiDecay2 = -100. ;
-  
-  TVector3 p3;  
-  TLorentzVector photonMom ;   
-  TObjArray * clusters  = 0x0 ;  
-  TObjArray * reftracks = 0x0;
-  Int_t nrefs           = 0;
-  Int_t nTracks         = GetCTSTracks()->GetEntriesFast() ;
-  
-  if(fPi0Trigger)
-  {
-    indexPhoton1 = aodParticle->GetCaloLabel (0);
-    indexPhoton2 = aodParticle->GetCaloLabel (1);
-    if(GetDebug() > 1)printf("indexPhoton1 = %d, indexPhoton2 = %d \n", indexPhoton1, indexPhoton2);
-    
-    if(indexPhoton1!=-1 && indexPhoton2!=-1)
-    {
-      if(aodParticle->GetDetector()=="EMCAL") clusters = GetEMCALClusters() ;
-      else                                    clusters = GetPHOSClusters()  ;
-      
-      for(Int_t iclus = 0; iclus < clusters->GetEntriesFast(); iclus++)
-      {
-        AliVCluster * photon =  (AliVCluster*) (clusters->At(iclus));  
-        photon->GetMomentum(photonMom,GetVertex(0)) ;
-        if(photon->GetID()==indexPhoton1) 
-        {
-          ptDecay1  = photonMom.Pt();
-          pxDecay1  = photonMom.Px();
-          pyDecay1  = photonMom.Py();
-          phiDecay1 = photonMom.Phi();
-          if(ptTrig && bFillHisto) fhPtPi0DecayRatio->Fill(ptTrig, ptDecay1/ptTrig);
-        }
-        if(photon->GetID()==indexPhoton2) 
-        {
-          ptDecay2  = photonMom.Pt();
-          pxDecay2  = photonMom.Px();
-          pyDecay2  = photonMom.Py();
-          phiDecay2 = photonMom.Phi();
-          if(ptTrig && bFillHisto) fhPtPi0DecayRatio->Fill(ptTrig, ptDecay2/ptTrig);
-        } 
-        
-        if(GetDebug() > 1)printf("Photon1 = %f, Photon2 = %f \n", ptDecay1, ptDecay2);
-        
-      } //cluster loop        
-    } //index of decay photons found
-  } //make decay-hadron correlation
-  
+  // In case of pi0/eta trigger, we may want to check their decay correlation, 
+  // get their decay children
+  TLorentzVector decayMom1;
+  TLorentzVector decayMom2;
+  Bool_t decayFound = kFALSE;
+  if(fPi0Trigger && bFillHisto) decayFound = GetDecayPhotonMomentum(aodParticle,decayMom1, decayMom2);
+
+  //-----------------------------------------------------------------------
   //Track loop, select tracks with good pt, phi and fill AODs or histograms
+  //-----------------------------------------------------------------------
+
   for(Int_t ipr = 0;ipr < pl->GetEntriesFast() ; ipr ++ )
   {
     AliVTrack * track = (AliVTrack *) (pl->At(ipr)) ;
@@ -1232,8 +1962,6 @@ Bool_t  AliAnaParticleHadronCorrelation::MakeChargedCorrelation(AliAODPWG4Partic
     Double_t mom[3] = {track->Px(),track->Py(),track->Pz()};
     p3.SetXYZ(mom[0],mom[1],mom[2]);
     pt   = p3.Pt();
-    px   = p3.Px();
-    py   = p3.Py();
     eta  = p3.Eta();
     phi  = p3.Phi() ;
     if(phi < 0) phi+=TMath::TwoPi();
@@ -1269,216 +1997,60 @@ Bool_t  AliAnaParticleHadronCorrelation::MakeChargedCorrelation(AliAODPWG4Partic
         return kFALSE;
     }    
     
-    if(fPi0Trigger)
-    {
-      if(indexPhoton1!=-1 && indexPhoton2!=-1)
-      {
-        if(ptDecay1) ratDecay1 = pt/ptDecay1 ;
-        if(ptDecay2) ratDecay2 = pt/ptDecay2 ; 
-        deltaPhiDecay1 = phiDecay1-phi;
-        deltaPhiDecay2 = phiDecay2-phi;
-        if(deltaPhiDecay1< -TMath::PiOver2()) deltaPhiDecay1+=TMath::TwoPi();
-        if(deltaPhiDecay1>3*TMath::PiOver2()) deltaPhiDecay1-=TMath::TwoPi();
-        if(deltaPhiDecay2< -TMath::PiOver2()) deltaPhiDecay2+=TMath::TwoPi();
-        if(deltaPhiDecay2>3*TMath::PiOver2()) deltaPhiDecay2-=TMath::TwoPi();    
-      }
-    } //do decay-hadron correlation    
-    
-    //Selection within angular range
-    deltaPhi    = phiTrig-phi;
-    deltaPhiOrg = deltaPhi;
-    if(deltaPhi <= -TMath::PiOver2()) deltaPhi+=TMath::TwoPi();
-    if(deltaPhi > 3*TMath::PiOver2()) deltaPhi-=TMath::TwoPi();
-    
-    zT   = pt/ptTrig ;
-    hbpZT = -100;
-    if(zT > 0 ) hbpZT = TMath::Log(1./zT); 
-
-    pout = pt*TMath::Sin(deltaPhi) ;
-    
-    xE   =-pt/ptTrig*TMath::Cos(deltaPhi); // -(px*pxTrig+py*pyTrig)/(ptTrig*ptTrig);
-    //if(xE <0.)xE =-xE;
-    hbpXE = -100;
-    if(xE > 0 ) hbpXE = TMath::Log(1./xE); 
-       
-    if(GetDebug() > 2)
-      printf("AliAnaParticleHadronCorrelation::MakeChargedCorrelation() - Charged hadron: pt %f, phi %f, phi trigger %f. Cuts:  delta phi  %2.2f < %2.2f < %2.2f, pT Trig min %2.2f \n",
-             pt,phi, phiTrig,fDeltaPhiMinCut, deltaPhi, fDeltaPhiMaxCut, fMinTriggerPt);
-    
     // Fill Histograms
     if(bFillHisto)
-    {
-      Int_t    assocBin   = -1; 
-      for(Int_t i = 0 ; i < fNAssocPtBins ; i++){
-        if(pt > fAssocPtBinLimit[i] && pt < fAssocPtBinLimit[i+1]) assocBin= i; 
-      }
-      
-      fhAssocPt->Fill(ptTrig,pt);
-      
-      //if(xE >= 1 ) printf("pTTrig = %f, pTHadron = %f, xE = %f\n",ptTrig,pt, xE);
-      fhXE->Fill(ptTrig, xE) ;
-      fhZT->Fill(ptTrig, zT) ;
+    {      
 
-      if(TMath::Cos(deltaPhi) < 0 && assocBin >= 0 )//away side 
+      if(GetDebug() > 2 ) 
+        printf("AliAnaParticleHadronCorrelation::MakeChargedCorrelation() - Selected charge for momentum imbalance: pt %2.2f, phi %2.2f, eta %2.2f \n",pt,phi,eta);
+            
+      // Set the pt associated bin for the defined bins
+      Int_t assocBin   = -1; 
+      for(Int_t i = 0 ; i < fNAssocPtBins ; i++)
       {
-        fhXEAssocPtBin[assocBin]->Fill(ptTrig, xE) ;
-        fhZTAssocPtBin[assocBin]->Fill(ptTrig, zT) ;
-      }
+        if(pt > fAssocPtBinLimit[i] && pt < fAssocPtBinLimit[i+1]) assocBin= i; 
+      }      
       
-      //Hardcoded values, BAD, FIXME
-      Double_t  dphiBrad = atan2(sin(deltaPhiOrg), cos(deltaPhiOrg))/TMath::Pi();//-1 to 1
-      if(TMath::Abs(dphiBrad)>0.325 && TMath::Abs(dphiBrad)<0.475)
-      {
-        fhAssocPtBkg->Fill(ptTrig, pt);
-      }
+      // Azimuthal Angle
+      // calculate deltaPhi for later, shift when needed
+      FillChargedAngularCorrelationHistograms(pt,  ptTrig,  assocBin, phi, phiTrig,  deltaPhi,
+                                              eta, etaTrig, decay, track->GetHMPIDsignal(),nTracks);
       
-      if(dphiBrad<-1./3) dphiBrad += 2;
-      fhDeltaPhiBrad->Fill(ptTrig, dphiBrad);
+      // Imbalance zT/xE/pOut
+      zT = pt/ptTrig ;
+      if(zT > 0 ) hbpZT = TMath::Log(1./zT);
+      else        hbpZT =-100;
       
-      if(assocBin>=0)
-      {
-        fhDeltaPhiBradAssocPtBin[assocBin]->Fill(ptTrig, dphiBrad);
-        fhDeltaPhiAssocPtBin    [assocBin]->Fill(ptTrig, deltaPhi);
-        if(track->GetHMPIDsignal()>0)
-        {
-          //printf("Track pt %f with HMPID signal %f \n",pt,track->GetHMPIDsignal());
-          fhDeltaPhiAssocPtBinHMPID[assocBin]->Fill(ptTrig, deltaPhi);        
-        }
-        
-        if(phi > 5*TMath::DegToRad() && phi < 20*TMath::DegToRad()){
-          //printf("Track pt %f in HMPID acceptance phi %f \n ",pt,phi*TMath::RadToDeg() );
-          fhDeltaPhiAssocPtBinHMPIDAcc[assocBin]->Fill(ptTrig, deltaPhi);        
-        }
-        
-      }
-      
-      fhEtaCharged     ->Fill(pt,eta);
-      fhPhiCharged     ->Fill(pt,phi);
-      fhDeltaEtaCharged->Fill(ptTrig,aodParticle->Eta()-eta);
-      fhDeltaPhiCharged->Fill(ptTrig, deltaPhi);
-      if(pt > 2 ) fhDeltaPhiDeltaEtaCharged->Fill(deltaPhi,aodParticle->Eta()-eta);
-      
-      if(GetDebug() > 2 ) printf("AliAnaParticleHadronCorrelation::MakeChargedCorrelation() - Selected charge for momentum imbalance: pt %2.2f, phi %2.2f, eta %2.2f \n",pt,phi,eta);
-      //fill different multiplicity histogram
-      if(DoEventSelect())
-      {
-        for(Int_t im=0; im<GetMultiBin(); im++){
-          if(nTracks < ( GetMaxMulti() - GetMinMulti() )/GetMultiBin()*(im+1)){
-            fhTrigDeltaPhiCharged[im]->Fill(ptTrig,deltaPhi);
-            fhTrigDeltaEtaCharged[im]->Fill(ptTrig,aodParticle->Eta()-eta);
-          }
-        }
-      }
+      xE   =-pt/ptTrig*TMath::Cos(deltaPhi); // -(px*pxTrig+py*pyTrig)/(ptTrig*ptTrig);
+      //if(xE <0.)xE =-xE;
+      if(xE > 0 ) hbpXE = TMath::Log(1./xE); 
+      else        hbpXE =-100;
+    
+      pout = pt*TMath::Sin(deltaPhi) ;
       
-      //delta phi cut for correlation
+      //delta phi cut for momentum imbalance correlation
       if      ( (deltaPhi > fDeltaPhiMinCut)   && (deltaPhi < fDeltaPhiMaxCut)   ) 
       {
-        fhDeltaPhiChargedPt ->Fill(pt    ,deltaPhi);
-        fhXECharged         ->Fill(ptTrig, xE); 
-        fhPtHbpXECharged    ->Fill(ptTrig, hbpXE);
-        fhZTCharged         ->Fill(ptTrig, zT); 
-        fhPtHbpZTCharged    ->Fill(ptTrig, hbpZT);
-        fhPtTrigPout        ->Fill(ptTrig, pout) ;
-        fhPtTrigCharged     ->Fill(ptTrig, pt) ;
         
-        if(track->Charge() > 0)
-        {
-          fhXEPosCharged->Fill(ptTrig,xE) ;
-          fhZTPosCharged->Fill(ptTrig,zT) ;
-        }
-        else  
-        {
-          fhXENegCharged->Fill(ptTrig,xE) ;
-          fhZTNegCharged->Fill(ptTrig,zT) ;
-        }
+        FillChargedMomentumImbalanceHistograms(ptTrig, pt, xE, hbpXE, zT, hbpZT, pout, deltaPhi, 
+                                               nTracks, track->Charge(), assocBin, decay);
         
-        //fill different multiplicity histogram
-        if(DoEventSelect())
-        {
-          for(Int_t im=0; im<GetMultiBin(); im++)
-          {
-            if(nTracks < ( GetMaxMulti() - GetMinMulti() )/GetMultiBin()*(im+1))
-            {
-              fhTrigXECorr[im]->Fill(ptTrig,xE);
-              fhTrigZTCorr[im]->Fill(ptTrig,zT);
-            }
-          }
-        } //multiplicity events selection
-      } //delta phi cut for correlation
-      else if ( (deltaPhi > fUeDeltaPhiMinCut) && (deltaPhi < fUeDeltaPhiMaxCut) ) 
+      } 
+      if ( (deltaPhi > fUeDeltaPhiMinCut) && (deltaPhi < fUeDeltaPhiMaxCut) ) 
       { //UE study
-        fhDeltaPhiUeChargedPt->Fill(pt,deltaPhi);
-      
-        Double_t randomphi = gRandom->Uniform(TMath::Pi()/2,3*TMath::Pi()/2);
-        Double_t uexE = -(pt/ptTrig)*TMath::Cos(randomphi);
-        Double_t uezT =   pt/ptTrig;
-
-        if(uexE < 0.) uexE = -uexE;
-       
-        if(GetDebug() > 1)printf("AliAnaParticleHadronCorrelation::MakeChargedCorrelation() - xe = %f, uexE = %f \n", xE, uexE);
         
-        fhXEUeCharged->Fill(ptTrig,uexE);
-        if(uexE>0)fhPtHbpXEUeCharged->Fill(ptTrig,TMath::Log(1/uexE));
-        
-        fhZTUeCharged->Fill(ptTrig,uezT);
-        if(uexE>0)fhPtHbpZTUeCharged->Fill(ptTrig,TMath::Log(1/uezT));
-        
-        if(DoEventSelect())
-        {
-          for(Int_t im=0; im<GetMultiBin(); im++)
-          {
-            if(nTracks < ( GetMaxMulti() - GetMinMulti() )/GetMultiBin()*(im+1))
-            {
-              fhTrigXEUeCorr[im]->Fill(ptTrig,xE);
-              fhTrigZTUeCorr[im]->Fill(ptTrig,zT);
-            }
-          }
-        } //multiplicity events selection
+        FillChargedUnderlyingEventHistograms(ptTrig, pt, deltaPhi, nTracks);
+
+       fhUePart->Fill(ptTrig);
         
-      } //UE study
+      }
       
-      if(fPi0Trigger)
-      {
-        if(indexPhoton1!=-1 && indexPhoton2!=-1)
-        {
-          fhDeltaPhiDecayCharged->Fill(ptDecay1, deltaPhiDecay1);
-          fhDeltaPhiDecayCharged->Fill(ptDecay2, deltaPhiDecay2);
-          
-          if(GetDebug() > 1)printf("AliAnaParticleHadronCorrelation::MakeChargedCorrelation() - deltaPhoton1 = %f, deltaPhoton2 = %f \n", deltaPhiDecay1, deltaPhiDecay2);
-          
-          if( (deltaPhiDecay1 > fDeltaPhiMinCut) && ( deltaPhiDecay1 < fDeltaPhiMaxCut) )
-            fhXEDecayCharged->Fill(ptDecay1,ratDecay1); 
-          
-          if( (deltaPhiDecay2 > fDeltaPhiMinCut) && ( deltaPhiDecay2 < fDeltaPhiMaxCut) )
-            fhXEDecayCharged->Fill(ptDecay2,ratDecay2);
-          
-          if(GetDebug() > 1)printf("AliAnaParticleHadronCorrelation::MakeChargedCorrelation() - ratPhoton1 = %f, ratPhoton2 = %f \n", pt/ptDecay1, pt/ptDecay2);
-          
-        } //index of decay photons found
-      } //make decay-hadron correlation          
+      if(fPi0Trigger && decayFound) 
+        FillDecayPhotonCorrelationHistograms(pt, phi, decayMom1,decayMom2, kTRUE) ;
       
       //several UE calculation 
-      if(fMakeSeveralUE)
-      {
-        if((deltaPhi<-fUeDeltaPhiMinCut) && (deltaPhi >-fUeDeltaPhiMaxCut))
-        {  
-          fhDeltaPhiUeLeftCharged->Fill(pt    ,deltaPhi);
-          fhXEUeLeftCharged      ->Fill(ptTrig,xE);
-          fhPtHbpXEUeLeftCharged ->Fill(ptTrig,hbpXE);
-          fhZTUeLeftCharged      ->Fill(ptTrig,zT);
-          fhPtHbpZTUeLeftCharged ->Fill(ptTrig,hbpZT);
-        }
-        
-        if((deltaPhi>fUeDeltaPhiMinCut) && (deltaPhi <fUeDeltaPhiMaxCut))
-        {  
-          fhDeltaPhiUeRightCharged->Fill(pt    ,deltaPhi);
-          fhXEUeRightCharged      ->Fill(ptTrig,xE);
-          fhPtHbpXEUeRightCharged ->Fill(ptTrig,hbpXE);
-          fhZTUeRightCharged      ->Fill(ptTrig,zT);
-          fhPtHbpZTUeRightCharged ->Fill(ptTrig,hbpZT);
-        }
-      } //several UE calculation
-            
+      if(fMakeSeveralUE) FillChargedUnderlyingEventSidesHistograms(ptTrig,pt,xE,hbpXE,zT,hbpZT,deltaPhi);
+      
     } //Fill histogram 
     else
     {
@@ -1488,9 +2060,11 @@ Bool_t  AliAnaParticleHadronCorrelation::MakeChargedCorrelation(AliAODPWG4Partic
         reftracks = new TObjArray(0);
         TString trackname = Form("%s+Tracks", GetAODObjArrayName().Data());
         reftracks->SetName(trackname.Data());
-        reftracks->SetOwner(kFALSE);
+        reftracks->SetOwner(kFALSE);        
       }
+      
       reftracks->Add(track);
+      
     }//aod particle loop
   }// track loop
   
@@ -1499,17 +2073,134 @@ Bool_t  AliAnaParticleHadronCorrelation::MakeChargedCorrelation(AliAODPWG4Partic
   {
     aodParticle->AddObjArray(reftracks);
   }
+
+  //Own mixed event, add event and remove previous or fill the mixed histograms
+  if(DoOwnMix() && bFillHisto)
+  {
+      MakeChargedMixCorrelation(aodParticle);
+  }
   
   return kTRUE;
   
 }  
 
+
+//______________________________________________________________________________________________________________________
+void AliAnaParticleHadronCorrelation::MakeChargedMixCorrelation(AliAODPWG4ParticleCorrelation *aodParticle)
+{  
+  // Mix current trigger with tracks in another MB event
+  
+  if(GetDebug() > 1) printf("AliAnaParticleHadronCorrelationNew::MakeChargedMixCorrelation() - Make trigger particle - charged hadron mixed event correlation \n");
+  
+  if(GetMixedEvent()) return;  // This is not the mixed event from general mixing frame
+  
+  // Get the event with similar caracteristics
+  //printf("MakeChargedMixCorrelation for %s\n",GetInputAODName().Data());
+
+  AliAnalysisManager   * manager      = AliAnalysisManager::GetAnalysisManager();
+  
+  AliInputEventHandler * inputHandler = dynamic_cast<AliInputEventHandler*>(manager->GetInputEventHandler());
+  
+  if(!inputHandler) return;
+  
+  if(!(inputHandler->IsEventSelected( ) & GetReader()->GetEventTriggerMask())) return;
+  
+  // Get the pool, check if it exits
+  Int_t eventBin = GetEventMixBin();
+
+  fhEventBin->Fill(eventBin);
+  
+  //Check that the bin exists, if not (bad determination of RP, centrality or vz bin) do nothing
+  if(eventBin < 0) return;
+  
+  TList * pool = 0;
+  if(fUseMixStoredInReader) pool = GetReader()->GetListWithMixedEventsForTracks(eventBin);
+  else                      pool = fListMixEvents[eventBin];
+
+  if(!pool) return ;
+    
+  Double_t ptTrig  = aodParticle->Pt();
+  Double_t etaTrig = aodParticle->Eta();
+  Double_t phiTrig = aodParticle->Phi();
+  if(phiTrig < 0.) phiTrig+=TMath::TwoPi();
+  
+  if(GetDebug() > 1) 
+    printf("AliAnaParticleHadronCorrelationNew::MakeChargedMixCorrelation() - Pool bin %d size %d, leading trigger pt=%f, phi=%f, eta=%f\n",
+           eventBin,pool->GetSize(), ptTrig,phiTrig,etaTrig);
+  
+  Double_t ptAssoc  = -999.;
+  Double_t phiAssoc = -999.;
+  Double_t etaAssoc = -999.;
+  Double_t deltaPhi = -999.;
+  Double_t deltaEta = -999.;
+  
+  //Start from first event in pool except if in this same event the pool was filled
+  Int_t ev0 = 0;
+  if(GetReader()->GetLastTracksMixedEvent() == GetEventNumber()) ev0 = 1;
+  for(Int_t ev=ev0; ev < pool->GetSize(); ev++)
+  {
+    TObjArray* bgTracks = static_cast<TObjArray*>(pool->At(ev));
+    
+    fhEventMixBin->Fill(eventBin);
+    
+    Int_t nTracks=bgTracks->GetEntriesFast();
+    //printf("\t Read Pool event %d, nTracks %d\n",ev,nTracks);
+
+    for(Int_t j1 = 0;j1 <nTracks; j1++ )
+    {
+      AliAODPWG4Particle *track = (AliAODPWG4Particle*) bgTracks->At(j1) ;
+      
+      if(!track) continue;
+      
+      ptAssoc  = track->Pt();
+      etaAssoc = track->Eta();
+      phiAssoc = track->Phi() ;
+      if(phiAssoc < 0) phiAssoc+=TMath::TwoPi();
+            
+      if(IsFiducialCutOn())
+      {
+        Bool_t in = GetFiducialCut()->IsInFiducialCut(*aodParticle->Momentum(),"CTS") ;
+        if(!in) continue ;
+      }
+      
+      deltaPhi = phiTrig-phiAssoc;
+      if(deltaPhi < -TMath::PiOver2())  deltaPhi+=TMath::TwoPi();
+      if(deltaPhi > 3*TMath::PiOver2()) deltaPhi-=TMath::TwoPi();
+      deltaEta = etaTrig-etaAssoc;
+      
+      if(GetDebug()>0)
+        printf("AliAnaParticleHadronCorrelationNew::MakeChargedMixCorrelation(): deltaPhi= %f, deltaEta=%f\n",deltaPhi, deltaEta);
+      
+      // Set the pt associated bin for the defined bins
+      Int_t assocBin   = -1; 
+      for(Int_t i = 0 ; i < fNAssocPtBins ; i++)
+      {
+        if(ptAssoc > fAssocPtBinLimit[i] && ptAssoc < fAssocPtBinLimit[i+1]) assocBin= i; 
+      }      
+
+      fhMixDeltaPhiCharged        ->Fill(ptTrig,  deltaPhi);
+      fhMixDeltaPhiDeltaEtaCharged->Fill(deltaPhi, deltaEta);
+
+      if(assocBin < 0) continue ; // this pt bin was not considered
+      
+      if(TMath::Abs(deltaEta) > 0.8) 
+        fhMixDeltaPhiChargedAssocPtBinDEta08  [assocBin]->Fill(ptTrig,   deltaPhi);
+      
+        fhMixDeltaPhiChargedAssocPtBin        [assocBin]->Fill(ptTrig,   deltaPhi);
+        fhMixDeltaPhiDeltaEtaChargedAssocPtBin[assocBin]->Fill(deltaPhi, deltaEta);
+      
+    } // track loop
+  } // mixed event loop
+}
+  
+
 //________________________________________________________________________________________________________________
 Bool_t  AliAnaParticleHadronCorrelation::MakeNeutralCorrelation(AliAODPWG4ParticleCorrelation * const aodParticle, 
                                                                 const TObjArray* pi0list, const Bool_t bFillHisto)  
 {  
   // Neutral Pion Correlation Analysis
-  if(GetDebug() > 1) printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelation() - Make trigger particle - pi0 correlation, %d pi0's \n",pi0list->GetEntriesFast());
+  if(GetDebug() > 1) printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelation() - Make trigger particle - pi0 correlation, %d pi0's \n",
+                            pi0list->GetEntriesFast());
   
   Int_t evtIndex11 = 0 ; 
   Int_t evtIndex12 = 0 ; 
@@ -1519,85 +2210,30 @@ Bool_t  AliAnaParticleHadronCorrelation::MakeNeutralCorrelation(AliAODPWG4Partic
     evtIndex12 = GetMixedEvent()->EventIndexForCaloCluster(aodParticle->GetCaloLabel(1)) ;    
   }  
   
-  Double_t pt   = -100. ;
-  Double_t px   = -100. ;
-  Double_t py   = -100. ;
-  Double_t zT   = -100. ; 
-  Double_t phi  = -100. ;
-  Double_t eta  = -100. ;
-  Double_t xE   = -100. ; 
-  Double_t hbpXE= -100. ; 
-  Double_t hbpZT= -100. ; 
+  Float_t pt   = -100. ;
+  Float_t zT   = -100. ; 
+  Float_t phi  = -100. ;
+  Float_t eta  = -100. ;
+  Float_t xE   = -100. ; 
+  Float_t hbpXE= -100. ; 
+  Float_t hbpZT= -100. ; 
+
+  Float_t ptTrig  = aodParticle->Pt();
+  Float_t phiTrig = aodParticle->Phi();
+  Float_t etaTrig = aodParticle->Eta();
+  Float_t deltaPhi= -100. ;
 
-  Double_t ptTrig  = aodParticle->Pt();
-  Double_t phiTrig = aodParticle->Phi();
-  Double_t etaTrig = aodParticle->Eta();
-  //Double_t pxTrig  = aodParticle->Px();
-  //Double_t pyTrig  = aodParticle->Py();
-  
-  Int_t indexPhoton1 =-1  ;
-  Int_t indexPhoton2 =-1  ;    
-  Double_t ptDecay1  = 0. ;
-  Double_t pxDecay1  = 0. ;
-  Double_t pyDecay1  = 0. ;
-  Double_t phiDecay1 = 0. ;
-  Double_t ptDecay2  = 0. ;
-  Double_t pxDecay2  = 0. ;
-  Double_t pyDecay2  = 0. ;
-  Double_t phiDecay2 = 0. ;
-  
-  Double_t ratDecay1      = -100. ;  
-  Double_t ratDecay2      = -100. ; 
-  Double_t deltaPhi       = -100. ;
-  Double_t deltaPhiDecay1 = -100. ;
-  Double_t deltaPhiDecay2 = -100. ;
-  
-  TObjArray * clusters = 0x0 ;  
   TLorentzVector photonMom ;
        
-  if(fPi0Trigger){
-    indexPhoton1 = aodParticle->GetCaloLabel (0);
-    indexPhoton2 = aodParticle->GetCaloLabel (1);
-    if(GetDebug() > 1)
-      printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillHistograms() - indexPhoton1 = %d, indexPhoton2 = %d \n", indexPhoton1, indexPhoton2);
-    
-    if(indexPhoton1!=-1 && indexPhoton2!=-1)
-    {
-      if(aodParticle->GetDetector()=="EMCAL") clusters = GetEMCALClusters() ;
-      else                                    clusters = GetPHOSClusters() ;
-      
-      for(Int_t iclus = 0; iclus < clusters->GetEntriesFast(); iclus++)
-      {
-        AliVCluster * photon =  (AliVCluster*) (clusters->At(iclus));  
-        photon->GetMomentum(photonMom,GetVertex(0)) ;
-        if(photon->GetID()==indexPhoton1) 
-        {
-          ptDecay1  = photonMom.Pt();
-          pxDecay1  = photonMom.Px();
-          pyDecay1  = photonMom.Py();
-          phiDecay1 = photonMom.Phi();
-        }
-        
-        if(photon->GetID()==indexPhoton2) 
-        {
-          ptDecay2  = photonMom.Pt();
-          pxDecay2  = photonMom.Px();
-          pyDecay2  = photonMom.Py();
-          phiDecay2 = photonMom.Phi();
-        } 
-        
-        if(GetDebug() > 1)
-          printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillHistograms() - Photon1 = %f, Photon2 = %f \n", ptDecay1, ptDecay2);
-        
-      } //photonAOD loop        
-    } //index of decay photons found
-    
-    if(ptTrig && bFillHisto) fhPtPi0DecayRatio->Fill(ptTrig, ptDecay1/ptTrig, ptDecay2/ptTrig);
-
-  } //make decay-hadron correlation
+  // In case of pi0/eta trigger, we may want to check their decay correlation, 
+  // get their decay children
+  TLorentzVector decayMom1;
+  TLorentzVector decayMom2;
+  Bool_t decayFound = kFALSE;
+  if(fPi0Trigger && bFillHisto) decayFound = GetDecayPhotonMomentum(aodParticle,decayMom1, decayMom2);  
   
-  TObjArray * refpi0    =0x0;
-  Int_t nrefs = 0;
+  TObjArray * refpi0 0x0;
+  Int_t nrefs        = 0;
   
   //Loop on stored AOD pi0
   
@@ -1620,34 +2256,29 @@ Bool_t  AliAnaParticleHadronCorrelation::MakeNeutralCorrelation(AliAODPWG4Partic
           evtIndex11 == evtIndex3 || evtIndex12 == evtIndex3) // trigger and pi0 are not from different events
         continue ; 
     }      
-    
-    //Int_t pdg = pi0->GetPdg();
-    //if(pdg != AliCaloPID::kPi0) continue;  
-    
+
     pt  = pi0->Pt();
-    px  = pi0->Px();
-    py  = pi0->Py();    
-    
+     
     if(pt < fMinAssocPt || pt > fMaxAssocPt) continue ;
     
-    //jumped out this event if near side associated particle pt larger than trigger
-    if(pt > ptTrig && TMath::Abs(phi-phiTrig)<TMath::PiOver2())  break ;
-    
-    //Selection within angular range
-    phi = pi0->Phi();
-    //Float_t deltaPhi = TMath::Abs(phiTrig-phi);
-    //if( (deltaPhi < fDeltaPhiMinCut) || ( deltaPhi > fDeltaPhiMaxCut) ) continue ;
+    //jump out this event if near side associated particle pt larger than trigger
+    if (fMakeNearSideLeading)
+    {
+      if(pt > ptTrig && TMath::Abs(phi-phiTrig) < TMath::PiOver2())  return kFALSE;
+    }
+    //jump out this event if there is any other particle with pt larger than trigger
+    else if(fMakeAbsoluteLeading)
+    {
+      if(pt > ptTrig)  return kFALSE;
+    }
     
     if(bFillHisto)
     {
-      
-      deltaPhi = phiTrig-phi;
-      if(deltaPhi<-TMath::PiOver2())  deltaPhi+=TMath::TwoPi();
-      if(deltaPhi>3*TMath::PiOver2()) deltaPhi-=TMath::TwoPi();
-      
       phi = pi0->Phi() ;
       eta = pi0->Eta() ;
-           
+      
+      FillNeutralAngularCorrelationHistograms(pt, ptTrig, phi, phiTrig, deltaPhi, eta, etaTrig);
+      
       zT  = pt/ptTrig ;
       xE  =-pt/ptTrig*TMath::Cos(deltaPhi); // -(px*pxTrig+py*pyTrig)/(ptTrig*ptTrig);
       
@@ -1655,38 +2286,12 @@ Bool_t  AliAnaParticleHadronCorrelation::MakeNeutralCorrelation(AliAODPWG4Partic
       
       hbpXE = -100;
       hbpZT = -100;
-
+      
       if(xE > 0 ) hbpXE = TMath::Log(1./xE); 
       if(zT > 0 ) hbpZT = TMath::Log(1./zT); 
-
-      if(fPi0Trigger)
-      {
-        if(indexPhoton1!=-1 && indexPhoton2!=-1)
-        {
-          if(ptDecay1) ratDecay1 = pt/ptDecay1 ;
-          if(ptDecay2) ratDecay2 = pt/ptDecay2 ; 
-          deltaPhiDecay1 = phiDecay1-phi;
-          deltaPhiDecay2 = phiDecay2-phi;
-          if(deltaPhiDecay1< -TMath::PiOver2()) deltaPhiDecay1+=TMath::TwoPi();
-          if(deltaPhiDecay1>3*TMath::PiOver2()) deltaPhiDecay1-=TMath::TwoPi();
-          if(deltaPhiDecay2< -TMath::PiOver2()) deltaPhiDecay2+=TMath::TwoPi();
-          if(deltaPhiDecay2>3*TMath::PiOver2()) deltaPhiDecay2-=TMath::TwoPi();   
-          fhDeltaPhiDecayNeutral->Fill(ptDecay1, deltaPhiDecay1);
-          fhDeltaPhiDecayNeutral->Fill(ptDecay2, deltaPhiDecay2);
-          if(GetDebug() > 1)printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelation() - deltaPhoton1 = %f, deltaPhoton2 = %f \n", deltaPhiDecay1, deltaPhiDecay2);
-          if( (deltaPhiDecay1 > fDeltaPhiMinCut) && ( deltaPhiDecay1 < fDeltaPhiMaxCut) )
-            fhXEDecayNeutral->Fill(ptDecay1,ratDecay1); 
-          if( (deltaPhiDecay2 > fDeltaPhiMinCut) && ( deltaPhiDecay2 < fDeltaPhiMaxCut) )
-            fhXEDecayNeutral->Fill(ptDecay2,ratDecay2);
-          if(GetDebug() > 1)printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelation() - ratPhoton1 = %f, ratPhoton2 = %f \n", pt/ptDecay1, pt/ptDecay2);
-        }
-      } //do decay-hadron correlation
       
-      fhEtaNeutral             ->Fill(pt,eta);
-      fhPhiNeutral             ->Fill(pt,phi);
-      fhDeltaEtaNeutral        ->Fill(ptTrig,etaTrig-eta);
-      fhDeltaPhiNeutral        ->Fill(ptTrig,deltaPhi);
-      fhDeltaPhiDeltaEtaNeutral->Fill(deltaPhi,etaTrig-eta);
+      if(fPi0Trigger && decayFound)
+        FillDecayPhotonCorrelationHistograms(pt, phi, decayMom1,decayMom2,kFALSE) ;
       
       //delta phi cut for correlation
       if( (deltaPhi > fDeltaPhiMinCut) && ( deltaPhi < fDeltaPhiMaxCut) ) 
@@ -1695,27 +2300,16 @@ Bool_t  AliAnaParticleHadronCorrelation::MakeNeutralCorrelation(AliAODPWG4Partic
         fhXENeutral        ->Fill(ptTrig,xE); 
         fhPtHbpXENeutral   ->Fill(ptTrig,hbpXE); 
       }
-      else 
+      else if ( (deltaPhi > fUeDeltaPhiMinCut) && (deltaPhi < fUeDeltaPhiMaxCut) )      
       {
         fhDeltaPhiUeNeutralPt->Fill(pt,deltaPhi);
         fhXEUeNeutral        ->Fill(ptTrig,xE);
         fhPtHbpXEUeNeutral   ->Fill(ptTrig,hbpXE); 
       }
+      
       //several UE calculation 
-      if(fMakeSeveralUE){
-        if((deltaPhi<-fUeDeltaPhiMinCut) && (deltaPhi >-fUeDeltaPhiMaxCut))
-        {  
-          fhDeltaPhiUeLeftNeutral->Fill(pt,deltaPhi);
-          fhXEUeLeftNeutral      ->Fill(ptTrig,xE);
-          fhPtHbpXEUeLeftNeutral ->Fill(ptTrig,hbpXE);
-        }
-        if((deltaPhi>fUeDeltaPhiMinCut) && (deltaPhi <fUeDeltaPhiMaxCut))
-        {  
-          fhDeltaPhiUeRightNeutral->Fill(pt,deltaPhi);
-          fhXEUeRightNeutral      ->Fill(ptTrig,xE);
-          fhPtHbpXEUeRightNeutral ->Fill(ptTrig,hbpXE);
-        }
-      } //several UE calculation
+      if(fMakeSeveralUE) FillChargedUnderlyingEventSidesHistograms(ptTrig,pt,xE,hbpXE,zT,hbpZT,deltaPhi);
+
          }
     else
     {
@@ -1741,6 +2335,7 @@ Bool_t  AliAnaParticleHadronCorrelation::MakeNeutralCorrelation(AliAODPWG4Partic
 void  AliAnaParticleHadronCorrelation::MakeMCChargedCorrelation(AliAODPWG4ParticleCorrelation *aodParticle)
 {  
   // Charged Hadron Correlation Analysis with MC information
+  
   if(GetDebug()>1)
     printf("AliAnaParticleHadronCorrelation::MakeMCChargedCorrelation() - Make trigger particle - charged hadron correlation in AOD MC level\n");
   
@@ -1754,27 +2349,10 @@ void  AliAnaParticleHadronCorrelation::MakeMCChargedCorrelation(AliAODPWG4Partic
   Double_t ptprim  = 0 ;
   Double_t phiprim = 0 ;
   Double_t etaprim = 0 ;
-  Double_t pxprim  = 0 ;
-  Double_t pyprim  = 0 ;
-  Double_t pzprim  = 0 ;
   Int_t    nTracks = 0 ;  
   Int_t iParticle  = 0 ;
   Double_t charge  = 0.;
-  
-  Double_t mczT    =-100 ;
-  Double_t mcxE    =-100 ;
-  Double_t mchbpXE =-100 ;
-  Double_t mchbpZT =-100 ;
-  Double_t mcdeltaPhi = -100. ;
 
-  //Track loop, select tracks with good pt, phi and fill AODs or histograms
-  //Int_t currentIndex = -1 ; 
-  Double_t mcTrackPt  = 0 ;
-  Double_t mcTrackPhi = 0 ;
-  Double_t mcTrackEta = 0 ;
-  Double_t mcTrackPx  = 0 ;
-  Double_t mcTrackPy  = 0 ;
-  Double_t mcTrackPz  = 0 ;
 
   if(GetReader()->ReadStack())
   {
@@ -1787,15 +2365,18 @@ void  AliAnaParticleHadronCorrelation::MakeMCChargedCorrelation(AliAODPWG4Partic
   //Int_t trackIndex[nTracks];
   
   Int_t label= aodParticle->GetLabel();
-  if(label<0){
-    printf("AliAnaParticleHadronCorrelation::MakeMCChargedCorrelation *** bad label ***:  label %d \n", label);
+  if(label < 0)
+  {
+    if(GetDebug() > 0) printf("AliAnaParticleHadronCorrelation::MakeMCChargedCorrelation *** bad label ***:  label %d \n", label);
     return;
   }  
+
   
   if(GetReader()->ReadStack())
   {
     stack =  GetMCStack() ;
-    if(!stack) {
+    if(!stack) 
+    {
       printf(" AliAnaParticleHadronCorrelation::MakeMCChargedCorrelation- Stack not available, is the MC handler called? STOP\n");
       abort();
     }
@@ -1820,9 +2401,6 @@ void  AliAnaParticleHadronCorrelation::MakeMCChargedCorrelation(AliAODPWG4Partic
       ptprim   = primary->Pt();
       phiprim  = primary->Phi();
       etaprim  = primary->Eta();
-      pxprim   = primary->Px();
-      pyprim   = primary->Py();
-      pzprim   = primary->Pz(); 
       
       if(ptprim < 0.01 || eprim < 0.01) return ;
       
@@ -1850,69 +2428,12 @@ void  AliAnaParticleHadronCorrelation::MakeMCChargedCorrelation(AliAODPWG4Partic
           
           if(inCTS)
           {            
-            mcTrackPt  = particle->Pt();
-            mcTrackPhi = particle->Phi();
-            mcTrackEta = particle->Eta();
-            mcTrackPx  = particle->Px();
-            mcTrackPy  = particle->Py();
-            mcTrackPz  = particle->Pz();              
-            if(mcTrackPhi < 0) mcTrackPhi+=TMath::TwoPi();    
-            
-            //Select only hadrons in pt range
-            if(mcTrackPt < fMinAssocPt || mcTrackPt > fMaxAssocPt) continue ;
-            
-            //remove trigger itself for correlation when use charged triggers 
-            if(label==iParticle && 
-               TMath::Abs(mcTrackPt -ptprim ) < 1e-6 && 
-               TMath::Abs(mcTrackPhi-phiprim) < 1e-6 && 
-               TMath::Abs(mcTrackEta-etaprim) < 1e-6) 
-              continue ;       
-            
-            //jumped out this event if near side associated partile pt larger than trigger
-            if( mcTrackPt > ptprim && TMath::Abs(mcTrackPhi-phiprim) < TMath::PiOver2()) 
-              return ;
-            
-            mcdeltaPhi= phiprim-mcTrackPhi; 
-            if(mcdeltaPhi <= -TMath::PiOver2()) mcdeltaPhi+=TMath::TwoPi();
-            if(mcdeltaPhi > 3*TMath::PiOver2()) mcdeltaPhi-=TMath::TwoPi();            
-            
-            mcdeltaPhi = phiprim-mcTrackPhi; 
-            mczT       = mcTrackPt/ptprim ;
-            mcxE       =-mcTrackPt/ptprim*TMath::Cos(mcdeltaPhi);// -(mcTrackPx*pxprim+mcTrackPy*pyprim)/(ptprim*ptprim);
-            mchbpXE    =-100 ;
-            mchbpZT    =-100 ;
-            if(mcxE > 0 ) mchbpXE = TMath::Log(1./mcxE);
-            if(mczT > 0 ) mchbpZT = TMath::Log(1./mczT);
-
-            // printf("rat = %f, xE = %f, cosi =%f \n", rat, xE, cosi);
-            // printf("phi = %f \n", phi);
-            
-            //Selection within angular range
-            if( mcdeltaPhi< -TMath::PiOver2())  mcdeltaPhi+=TMath::TwoPi();
-            if( mcdeltaPhi>3*TMath::PiOver2())  mcdeltaPhi-=TMath::TwoPi();              
-            Double_t mcpout = mcTrackPt*TMath::Sin(mcdeltaPhi) ;              
-            if(GetDebug()>0 )  
-              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",
-                     mcTrackPt,mcTrackPhi, phiprim,fDeltaPhiMinCut, mcdeltaPhi, fDeltaPhiMaxCut, GetMinPt());              
-            // Fill Histograms
-            fhMCEtaCharged->Fill(mcTrackPt,mcTrackEta);
-            fhMCPhiCharged->Fill(mcTrackPt,mcTrackPhi);
-            fhMCDeltaEtaCharged->Fill(ptprim,etaprim-mcTrackEta);
-            fhMCDeltaPhiCharged->Fill(ptprim,mcdeltaPhi);
-            fhMCPtAssocDeltaPhi->Fill(mcTrackPt, mcdeltaPhi);
-            fhMCDeltaPhiDeltaEtaCharged->Fill(mcdeltaPhi,etaprim-mcTrackEta);
-            
-            //delta phi cut for correlation
-            if( (mcdeltaPhi > fDeltaPhiMinCut) && ( mcdeltaPhi < fDeltaPhiMaxCut) ) {
-              fhMCDeltaPhiChargedPt->Fill(mcTrackPt,mcdeltaPhi);
-              fhMCPtXECharged      ->Fill(ptprim,mcxE); 
-              fhMCPtHbpXECharged   ->Fill(ptprim,mchbpXE);
-              fhMCPtZTCharged      ->Fill(ptprim,mczT); 
-              fhMCPtHbpZTCharged   ->Fill(ptprim,mchbpZT);
-              fhMCPtTrigPout       ->Fill(ptprim, mcpout) ;
-            }//delta phi cut for correlation
-          } //tracks after cuts
-        }//Charged
+            if( label!=iParticle) // avoid trigger particle
+            {
+              if(!FillChargedMCCorrelationHistograms(particle->Pt(),particle->Phi(),particle->Eta(),ptprim,phiprim,etaprim)) return;
+            }
+          }// in CTS acceptance
+        }// charged
       } //track loop
     } //when the leading particles could trace back to MC
   } //ESD MC
@@ -1934,22 +2455,24 @@ void  AliAnaParticleHadronCorrelation::MakeMCChargedCorrelation(AliAODPWG4Partic
       printf("AliAnaParticleHadronCorrelation::MakeMCChargedCorrelation *** no AOD primary ***:  label %d \n", label);   
       return;
     }
-    
+
+   
     if(aodprimary)
     {
       ptprim  = aodprimary->Pt();
       phiprim = aodprimary->Phi();
       etaprim = aodprimary->Eta();
-      pxprim  = aodprimary->Px();
-      pyprim  = aodprimary->Py();
-      pzprim  = aodprimary->Pz();  
+      eprim   = aodprimary->E();
+
+      Bool_t lead = kFALSE;
       
       if(ptprim < 0.01 || eprim < 0.01) return ;
-      
       mcparticles= GetReader()->GetAODMCParticles();
-      for (Int_t i=0; i<nTracks;i++) 
+      for (Int_t i = 0; i < nTracks; i++) 
       {
         AliAODMCParticle *part = (AliAODMCParticle*) mcparticles->At(i);
+        
         if (!part->IsPhysicalPrimary()) continue;        
         Int_t pdg = part->GetPdgCode();        
         charge = TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
@@ -1969,70 +2492,16 @@ void  AliAnaParticleHadronCorrelation::MakeMCChargedCorrelation(AliAODPWG4Partic
             
             if(inCTS)
             {            
-              mcTrackPt  =part->Pt();
-              mcTrackPhi =part->Phi();
-              mcTrackEta =part->Eta();
-              mcTrackPx  =part->Px();
-              mcTrackPy  =part->Py();
-              mcTrackPz  =part->Pz();              
-              if(mcTrackPhi < 0) mcTrackPhi+=TMath::TwoPi();   
-              
-              //Select only hadrons in pt range
-              if(mcTrackPt < GetMinPt() || mcTrackPt > GetMaxPt()) continue ;
-              
-              //remove trigger itself for correlation when use charged triggers 
-              if(label==i && 
-                 TMath::Abs(mcTrackPt -ptprim ) < 1e-6 && 
-                 TMath::Abs(mcTrackPhi-phiprim) < 1e-6 && 
-                 TMath::Abs(mcTrackEta-etaprim) < 1e-6) continue ;  
-              
-              //jumped out this event if near side associated partile pt larger than trigger
-              if( mcTrackPt> ptprim && TMath::Abs(mcTrackPhi-phiprim)<TMath::PiOver2()) 
-                return ;
-              
-              mcdeltaPhi= phiprim-mcTrackPhi; 
-              if(mcdeltaPhi <= -TMath::PiOver2()) mcdeltaPhi+=TMath::TwoPi();
-              if(mcdeltaPhi > 3*TMath::PiOver2()) mcdeltaPhi-=TMath::TwoPi();
-              
-              mczT      = mcTrackPt/ptprim ;
-              mcxE      =-mcTrackPt/ptprim*TMath::Cos(mcdeltaPhi);// -(mcTrackPx*pxprim+mcTrackPy*pyprim)/(ptprim*ptprim);
-              mchbpXE   =-100 ;
-              mchbpZT   =-100 ;
-              if(mcxE > 0 ) mchbpXE = TMath::Log(1./mcxE);
-              if(mczT > 0 ) mchbpZT = TMath::Log(1./mczT);
-
-              //Selection within angular range
-              if( mcdeltaPhi< -TMath::PiOver2())  mcdeltaPhi+=TMath::TwoPi();
-              if( mcdeltaPhi>3*TMath::PiOver2())  mcdeltaPhi-=TMath::TwoPi();              
-              Double_t mcpout = mcTrackPt*TMath::Sin(mcdeltaPhi) ;              
-              if(GetDebug()>0)  
-                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",
-                       mcTrackPt,mcTrackPhi, phiprim,fDeltaPhiMinCut, mcdeltaPhi, fDeltaPhiMaxCut, GetMinPt());              
-
-              // Fill Histograms
-
-              fhMCEtaCharged     ->Fill(mcTrackPt,mcTrackEta);
-              fhMCPhiCharged     ->Fill(mcTrackPt,mcTrackPhi);
-              fhMCDeltaEtaCharged->Fill(ptprim,etaprim-mcTrackEta);
-              fhMCDeltaPhiCharged->Fill(ptprim,mcdeltaPhi);
-              fhMCPtAssocDeltaPhi->Fill(mcTrackPt, mcdeltaPhi);
-              fhMCDeltaPhiDeltaEtaCharged->Fill(mcdeltaPhi,etaprim-mcTrackEta);
-              
-              //delta phi cut for correlation
-              if( (mcdeltaPhi > fDeltaPhiMinCut) && ( mcdeltaPhi < fDeltaPhiMaxCut) ) {
-                fhMCDeltaPhiChargedPt ->Fill(mcTrackPt,mcdeltaPhi);
-                fhMCPtXECharged       ->Fill(ptprim,mcxE); 
-                fhMCPtHbpXECharged    ->Fill(ptprim,mchbpXE);
-                fhMCPtZTCharged       ->Fill(ptprim,mczT); 
-                fhMCPtHbpZTCharged    ->Fill(ptprim,mchbpZT);
-                fhMCPtTrigPout        ->Fill(ptprim, mcpout) ;
-              }//delta phi cut for correlation
-              
-            } //tracks after cuts
-            
-          } //with minimum pt cut
+              if( label!=iParticle) // avoid trigger particle
+              {
+                if(!FillChargedMCCorrelationHistograms(part->Pt(),part->Phi(),part->Eta(),ptprim,phiprim,etaprim)) return;
+               else lead = kTRUE;
+              }
+            } // in acceptance
+          } // min pt cut
         } //only charged particles
-      }  //MC particle loop      
+      }  //MC particle loop    
+      if (lead) fhMCPtLeading->Fill(ptprim);
     } //when the leading particles could trace back to MC
   }// AOD MC
 }
@@ -2095,9 +2564,9 @@ void AliAnaParticleHadronCorrelation::SetAssocPtBinLimit(Int_t ibin, Float_t pt)
   {
     fAssocPtBinLimit[ibin] = pt  ;
   }
-  else {
+  else 
+  {
     printf("AliAnaParticleHadronCorrelation::SetAssocPtBinLimit() - bin  number too large %d > %d or small, nothing done\n", ibin, fNAssocPtBins) ; 
-    
   }
 }