]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWG4/PartCorrDep/AliAnaParticleHadronCorrelation.cxx
Updates from Yaxian
[u/mrichter/AliRoot.git] / PWG4 / PartCorrDep / AliAnaParticleHadronCorrelation.cxx
index bc5c925ea37e8be189c0a46120f6f85ff8122592..224052a571257df94b90f828c29151a44d68c238 100755 (executable)
 // Class for the analysis of particle - hadron correlations
 // Particle (for example direct gamma) must be found in a previous analysis 
 //-- Author: Gustavo Conesa (LNF-INFN) 
+
+//  Modified by Yaxian Mao:
+// 1. add the UE subtraction for corrlation study
+// 2. change the correlation variable
+// 3. Only use leading particle(cluster/track) as trigger for correlation (2010/07/02)
+// 4. Make decay photon-hadron correlations where decay contribute pi0 mass (2010/09/09)
+// 5. fill the pout to extract kt at the end, also to study charge asymmetry(2010/10/06) 
 //////////////////////////////////////////////////////////////////////////////
 
 
 // --- ROOT system ---
-#include "TH2F.h"
+//#include "TClonesArray.h"
+#include "TClass.h"
+#include "TMath.h"
+#include "TH3D.h"
 
 //---- ANALYSIS system ----
-#include "AliLog.h"
 #include "AliNeutralMesonSelection.h" 
 #include "AliAnaParticleHadronCorrelation.h" 
 #include "AliCaloTrackReader.h"
 #include "AliCaloPID.h"
 #include "AliAODPWG4ParticleCorrelation.h"
+#include "AliFiducialCut.h"
+#include "AliAODTrack.h"
+#include "AliVCluster.h"
+#include "AliMCAnalysisUtils.h"
+#include "TParticle.h"
+#include "AliStack.h"
+#include "AliAODMCParticle.h"
+#include "AliMixedEvent.h"
 
 ClassImp(AliAnaParticleHadronCorrelation)
 
 
 //____________________________________________________________________________
-  AliAnaParticleHadronCorrelation::AliAnaParticleHadronCorrelation() 
+  AliAnaParticleHadronCorrelation::AliAnaParticleHadronCorrelation(): 
     AliAnaPartCorrBaseClass(),
-    fDeltaPhiMaxCut(0.), fDeltaPhiMinCut(0.), 
-    fhPhiCharged(0), fhPhiNeutral(0), fhEtaCharged(0), fhEtaNeutral(0), 
-    fhDeltaPhiCharged(0), fhDeltaPhiNeutral(0), 
-    fhDeltaEtaCharged(0), fhDeltaEtaNeutral(0),
-    fhDeltaPhiChargedPt(0), fhDeltaPhiNeutralPt(0), 
-    fhPtImbalanceNeutral(0), fhPtImbalanceCharged(0)
+    fDeltaPhiMaxCut(0.), fDeltaPhiMinCut(0.), fSelectIsolated(0),
+    fMakeSeveralUE(0),  fUeDeltaPhiMaxCut(0.), fUeDeltaPhiMinCut(0.), 
+    fPi0AODBranchName(""),fNeutralCorr(0), fPi0Trigger(0),
+   // fMultiBin(0),fNZvertBin(0),fNrpBin(0),fZvtxCut(0.),
+   // fUseSelectEvent(kFALSE), 
+   // fhNclustersNtracks(0), //fhVertex(0),
+    fhPtLeading(0),fhPhiLeading(0),fhEtaLeading(0), 
+    fhDeltaPhiDeltaEtaCharged(0),
+    fhPhiCharged(0), fhEtaCharged(0), 
+    fhDeltaPhiCharged(0), 
+    fhDeltaEtaCharged(0), 
+    fhDeltaPhiChargedPt(0), 
+    fhDeltaPhiUeChargedPt(0), 
+    fhPtImbalanceCharged(0), 
+    fhPtImbalanceUeCharged(0),
+    fhPtImbalancePosCharged(0),fhPtImbalanceNegCharged(0),
+    fhPtHbpCharged(0), fhPtHbpUeCharged(0),
+    fhDeltaPhiUeLeftCharged(0),fhDeltaPhiUeRightCharged(0),
+    fhPtImbalanceUeLeftCharged(0),fhPtImbalanceUeRightCharged(0),
+    fhPtHbpUeLeftCharged(0),fhPtHbpUeRightCharged(0), 
+    fhPoutPtTrigPtAssoc(0), fhUePoutPtTrigPtAssoc(0),
+    fhPtTrigCharged(0),
+    fhTrigDeltaPhiDeltaEtaCharged(0x0), fhTrigCorr(0x0),fhTrigUeCorr(0x0),
+    fhDeltaPhiDeltaEtaNeutral(0), 
+    fhPhiNeutral(0), fhEtaNeutral(0), 
+    fhDeltaPhiNeutral(0), fhDeltaEtaNeutral(0),
+    fhDeltaPhiNeutralPt(0),fhDeltaPhiUeNeutralPt(0), 
+    fhPtImbalanceNeutral(0),fhPtImbalanceUeNeutral(0),
+    fhPtHbpNeutral(0), fhPtHbpUeNeutral(0),
+    fhDeltaPhiUeLeftNeutral(0),fhDeltaPhiUeRightNeutral(0),
+    fhPtImbalanceUeLeftNeutral(0),fhPtImbalanceUeRightNeutral(0),
+    fhPtHbpUeLeftNeutral(0),fhPtHbpUeRightNeutral(0),
+    fhPtPi0DecayRatio(0),
+    fhDeltaPhiDecayCharged(0),
+    fhPtImbalanceDecayCharged(0), 
+    fhDeltaPhiDecayNeutral(0),
+    fhPtImbalanceDecayNeutral(0)
 {
   //Default Ctor
+  
   //Initialize parameters
   InitParameters();
 }
 
-//____________________________________________________________________________
-AliAnaParticleHadronCorrelation::AliAnaParticleHadronCorrelation(const AliAnaParticleHadronCorrelation & g) :   
-  AliAnaPartCorrBaseClass(g),
-  fDeltaPhiMaxCut(g.fDeltaPhiMaxCut), fDeltaPhiMinCut(g.fDeltaPhiMinCut), 
-  fhPhiCharged(g.fhPhiCharged), fhPhiNeutral(g.fhPhiNeutral), 
-  fhEtaCharged(g.fhEtaCharged), fhEtaNeutral(g.fhEtaNeutral), 
-  fhDeltaPhiCharged(g.fhDeltaPhiCharged),  
-  fhDeltaPhiNeutral(g.fhDeltaPhiNeutral), 
-  fhDeltaEtaCharged(g.fhDeltaEtaCharged), 
-  fhDeltaEtaNeutral(g.fhDeltaEtaNeutral), 
-  fhDeltaPhiChargedPt(g.fhDeltaPhiChargedPt), 
-  fhDeltaPhiNeutralPt(g.fhDeltaPhiNeutralPt), 
-  fhPtImbalanceNeutral(g.fhPtImbalanceNeutral), 
-  fhPtImbalanceCharged(g.fhPtImbalanceCharged)
-{
-  // cpy ctor
-
-}
-
-//_________________________________________________________________________
-AliAnaParticleHadronCorrelation & AliAnaParticleHadronCorrelation::operator = (const AliAnaParticleHadronCorrelation & source)
-{
-  // assignment operator
-
-  if(this == &source)return *this;
-  ((AliAnaPartCorrBaseClass *)this)->operator=(source);
+//________________________________________________________________________
+TList *  AliAnaParticleHadronCorrelation::GetCreateOutputObjects()
+{  
   
-  fDeltaPhiMaxCut = source.fDeltaPhiMaxCut ; 
-  fDeltaPhiMinCut = source.fDeltaPhiMinCut ; 
-
-  fhPhiCharged = source.fhPhiCharged ; fhPhiNeutral = source.fhPhiNeutral ; 
-  fhEtaCharged = source.fhEtaCharged ; fhEtaNeutral = source.fhEtaNeutral ; 
-  fhDeltaPhiCharged = source.fhDeltaPhiCharged ;  
-  fhDeltaPhiNeutral = source.fhDeltaPhiNeutral ; 
-  fhDeltaPhiNeutralPt = source.fhDeltaPhiNeutralPt ; 
-  fhDeltaEtaCharged = source.fhDeltaEtaCharged ; 
-  fhDeltaEtaNeutral = source.fhDeltaEtaNeutral ; 
-  fhDeltaPhiChargedPt = source.fhDeltaPhiChargedPt ;
-
-  fhPtImbalanceNeutral = source.fhPtImbalanceNeutral ; 
-  fhPtImbalanceCharged = source.fhPtImbalanceCharged ; 
+  // Create histograms to be saved in output file and 
+  // store them in fOutputContainer
+  TList * outputContainer = new TList() ; 
+  outputContainer->SetName("CorrelationHistos") ; 
+  
+  Int_t nptbins  = GetHistoPtBins();
+  Int_t nphibins = GetHistoPhiBins();
+  Int_t netabins = GetHistoEtaBins();
+  Float_t ptmax  = GetHistoPtMax();
+  Float_t phimax = GetHistoPhiMax();
+  Float_t etamax = GetHistoEtaMax();
+  Float_t ptmin  = GetHistoPtMin();
+  Float_t phimin = GetHistoPhiMin();
+  Float_t etamin = GetHistoEtaMin();   
 
-  return *this;
+  
+//  fhNclustersNtracks  = new TH2F ("Multiplicity","Neutral cluster and charged track multiplicity",1000, 0., 1000.,1000, 0., 1000.); 
+//  fhNclustersNtracks->SetYTitle("# of tracks");
+//  fhNclustersNtracks->SetXTitle("# of clusters");
 
-}
+  fhPtLeading  = new TH1F ("hPtLeading","p_T distribution of leading particles", nptbins,ptmin,ptmax); 
+  fhPtLeading->SetXTitle("p_{T}^{trig} (GeV/c)");
+  
+  fhPhiLeading  = new TH2F ("hPhiLeading","#phi distribution of leading Particles",nptbins,ptmin,ptmax, nphibins,phimin,phimax); 
+  fhPhiLeading->SetYTitle("#phi (rad)");
+  
+  fhEtaLeading  = new TH2F ("hEtaLeading","#eta distribution of leading",nptbins,ptmin,ptmax, netabins,etamin,etamax); 
+  fhEtaLeading->SetYTitle("#eta ");  
+ // outputContainer->Add(fhNclustersNtracks);
+  outputContainer->Add(fhPtLeading);
+  outputContainer->Add(fhPhiLeading);
+  outputContainer->Add(fhEtaLeading);
+  
+  //Correlation with charged hadrons
+  if(GetReader()->IsCTSSwitchedOn()) {
+    fhDeltaPhiDeltaEtaCharged  = new TH2F
+    ("DeltaPhiDeltaEtaCharged","#phi_{trigger} - #phi_{h^{#pm}} vs #eta_{trigger} - #eta_{h^{#pm}}",
+     140,-2.,5.,200,-2,2); 
+    fhDeltaPhiDeltaEtaCharged->SetXTitle("#Delta #phi");
+    fhDeltaPhiDeltaEtaCharged->SetYTitle("#Delta #eta");    
+    
+    fhPhiCharged  = new TH2F
+    ("PhiCharged","#phi_{h^{#pm}}  vs p_{T #pm}",
+     nptbins,ptmin,ptmax,nphibins,phimin,phimax); 
+    fhPhiCharged->SetYTitle("#phi_{h^{#pm}} (rad)");
+    fhPhiCharged->SetXTitle("p_{T #pm} (GeV/c)");
+    
+    fhEtaCharged  = new TH2F
+    ("EtaCharged","#eta_{h^{#pm}}  vs p_{T #pm}",
+     nptbins,ptmin,ptmax,netabins,etamin,etamax); 
+    fhEtaCharged->SetYTitle("#eta_{h^{#pm}} (rad)");
+    fhEtaCharged->SetXTitle("p_{T #pm} (GeV/c)");
+    
+    fhDeltaPhiCharged  = new TH2F
+    ("DeltaPhiCharged","#phi_{trigger} - #phi_{h^{#pm}} vs p_{T trigger}",
+     nptbins,ptmin,ptmax,140,-2.,5.); 
+    fhDeltaPhiCharged->SetYTitle("#Delta #phi");
+    fhDeltaPhiCharged->SetXTitle("p_{T trigger} (GeV/c)");
+    
+    fhDeltaPhiChargedPt  = new TH2F
+    ("DeltaPhiChargedPt","#phi_{trigger} - #phi_{#h^{#pm}} vs p_{T h^{#pm}}",
+     nptbins,ptmin,ptmax,140,-2.,5.);
+    fhDeltaPhiChargedPt->SetYTitle("#Delta #phi");
+    fhDeltaPhiChargedPt->SetXTitle("p_{T h^{#pm}} (GeV/c)");
+    
+    fhDeltaPhiUeChargedPt  = new TH2F
+    ("DeltaPhiUeChargedPt","#phi_{trigger} - #phi_{#Ueh^{#pm}} vs p_{T Ueh^{#pm}}",
+     nptbins,ptmin,ptmax,140,-2.,5.);
+    fhDeltaPhiUeChargedPt->SetYTitle("#Delta #phi");
+    fhDeltaPhiUeChargedPt->SetXTitle("p_{T h^{#pm}} (GeV/c)");
+    
+    fhDeltaEtaCharged  = new TH2F
+    ("DeltaEtaCharged","#eta_{trigger} - #eta_{h^{#pm}} vs p_{T trigger}",
+     nptbins,ptmin,ptmax,200,-2,2); 
+    fhDeltaEtaCharged->SetYTitle("#Delta #eta");
+    fhDeltaEtaCharged->SetXTitle("p_{T trigger} (GeV/c)");
+    
+    fhPtImbalanceCharged  = 
+    new TH2F("CorrelationCharged","z_{trigger h^{#pm}} = p_{T h^{#pm}} / p_{T trigger}",
+             nptbins,ptmin,ptmax,200,0.,2.); 
+    fhPtImbalanceCharged->SetYTitle("z_{trigger h^{#pm}}");
+    fhPtImbalanceCharged->SetXTitle("p_{T trigger}");
+    
+    fhPtImbalanceUeCharged  = 
+    new TH2F("CorrelationUeCharged","z_{trigger h^{#pm}} = p_{T Ueh^{#pm}} / p_{T trigger}",
+             nptbins,ptmin,ptmax,200,0.,2.); 
+    fhPtImbalanceUeCharged->SetYTitle("z_{trigger Ueh^{#pm}}");
+    fhPtImbalanceUeCharged->SetXTitle("p_{T trigger}");
+    
+    fhPtImbalancePosCharged  = 
+    new TH2F("CorrelationPositiveCharged","z_{trigger h^{+}} = p_{T h^{+}} / p_{T trigger}",
+             nptbins,ptmin,ptmax,200,0.,2.); 
+    fhPtImbalancePosCharged->SetYTitle("z_{trigger h^{+}}");
+    fhPtImbalancePosCharged->SetXTitle("p_{T trigger}");
+    
+    fhPtImbalanceNegCharged  = 
+    new TH2F("CorrelationNegativeCharged","z_{trigger h^{-}} = p_{T h^{-}} / p_{T trigger}",
+             nptbins,ptmin,ptmax,200,0.,2.); 
+    fhPtImbalanceNegCharged->SetYTitle("z_{trigger h^{-}}");
+    fhPtImbalanceNegCharged->SetXTitle("p_{T trigger}");
+    
+    fhPtHbpCharged  = 
+    new TH2F("HbpCharged","#xi = ln(1/x_{E}) with charged hadrons",
+             nptbins,ptmin,ptmax,200,0.,10.); 
+    fhPtHbpCharged->SetYTitle("ln(1/x_{E})");
+    fhPtHbpCharged->SetXTitle("p_{T trigger}");
+    
+    fhPtHbpUeCharged  = 
+    new TH2F("HbpUeCharged","#xi = ln(1/x_{E}) with charged hadrons",
+             nptbins,ptmin,ptmax,200,0.,10.); 
+    fhPtHbpUeCharged->SetYTitle("ln(1/x_{E})");
+    fhPtHbpUeCharged->SetXTitle("p_{T trigger}");
+    
+    fhPoutPtTrigPtAssoc  = 
+    new TH3D("PoutPtTrigPtAssoc","Pout with charged hadrons",
+             nptbins,ptmin,ptmax,nptbins,ptmin,ptmax,2*nptbins,-ptmax,ptmax); 
+    fhPoutPtTrigPtAssoc->SetZTitle("p_{out} (GeV/c)");
+    fhPoutPtTrigPtAssoc->SetYTitle("p_{T associated} (GeV/c)");  
+    fhPoutPtTrigPtAssoc->SetXTitle("p_{T trigger} (GeV/c)"); 
+    
+    fhUePoutPtTrigPtAssoc  = 
+    new TH3D("UePoutPtTrigPtAssoc"," UE Pout with charged hadrons",
+             nptbins,ptmin,ptmax,nptbins,ptmin,ptmax,2*nptbins,-ptmax,ptmax); 
+    fhUePoutPtTrigPtAssoc->SetZTitle("p_{out} (GeV/c)");
+    fhUePoutPtTrigPtAssoc->SetYTitle("p_{T associated} (GeV/c)");  
+    fhUePoutPtTrigPtAssoc->SetXTitle("p_{T trigger} (GeV/c)"); 
+    
+    fhPtTrigCharged  = 
+    new TH2F("PtTrigCharged","trgger and charged tracks pt distribution",
+             nptbins,ptmin,ptmax,nptbins,ptmin,ptmax); 
+    fhPtTrigCharged->SetYTitle("p_{T asso} (GeV/c)");
+    fhPtTrigCharged->SetXTitle("p_{T trigger} (GeV/c)");    
+         
+    outputContainer->Add(fhDeltaPhiDeltaEtaCharged);
+    outputContainer->Add(fhPhiCharged) ;
+    outputContainer->Add(fhEtaCharged) ;
+    outputContainer->Add(fhDeltaPhiCharged) ; 
+    outputContainer->Add(fhDeltaEtaCharged) ;
+    outputContainer->Add(fhDeltaPhiChargedPt) ;
+    outputContainer->Add(fhDeltaPhiUeChargedPt) ;
+    outputContainer->Add(fhPtImbalanceCharged) ;
+    outputContainer->Add(fhPtImbalancePosCharged) ;
+    outputContainer->Add(fhPtImbalanceNegCharged) ;
+    outputContainer->Add(fhPtImbalanceUeCharged) ;
+    outputContainer->Add(fhPtHbpCharged) ;
+    outputContainer->Add(fhPtHbpUeCharged) ;
+    outputContainer->Add(fhPoutPtTrigPtAssoc) ;
+    outputContainer->Add(fhUePoutPtTrigPtAssoc) ;
+    outputContainer->Add(fhPtTrigCharged) ;
+    
+    if(DoEventSelect()){ 
+      Int_t nMultiBins = GetMultiBin();
+      fhTrigDeltaPhiDeltaEtaCharged = new TH3D*[nMultiBins] ;
+      fhTrigCorr = new TH2F*[nMultiBins];
+      fhTrigUeCorr = new TH2F*[nMultiBins];
+      for(Int_t im=0; im<nMultiBins; im++){
+        fhTrigDeltaPhiDeltaEtaCharged[im]  = new TH3D 
+        (Form("fhTrigDeltaPhiDeltaEtaCharged_%d",im),Form("fhTrigDeltaPhiDeltaEtaCharged_%d",im), nptbins,ptmin,ptmax, 140,-2.,5., 200,-2,2); 
+        fhTrigDeltaPhiDeltaEtaCharged[im]->SetXTitle("p_{T trigger} (GeV/c)");
+        fhTrigDeltaPhiDeltaEtaCharged[im]->SetYTitle("#Delta #phi");
+        fhTrigDeltaPhiDeltaEtaCharged[im]->SetZTitle("#Delta #eta");
+        fhTrigCorr[im]  = new TH2F
+        (Form("fhTrigPtCorr_%d",im),Form("fhTrigPtCorr_%d",im), nptbins,ptmin,ptmax,200,0.,2.); 
+        fhTrigCorr[im]->SetYTitle("z_{trigger h^{#pm}}");
+        fhTrigCorr[im]->SetXTitle("p_{T trigger}");
+        fhTrigUeCorr[im]  = new TH2F
+        (Form("fhTrigPtUeCorr_%d",im),Form("fhTrigPtUeCorr_%d",im), nptbins,ptmin,ptmax,200,0.,2.); 
+        fhTrigUeCorr[im]->SetYTitle("z_{trigger h^{#pm}}");
+        fhTrigUeCorr[im]->SetXTitle("p_{T trigger}");       
 
-//________________________________________________________________________
-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  = GetHistoNPtBins();
-       Int_t nphibins = GetHistoNPhiBins();
-       Int_t netabins = GetHistoNEtaBins();
-       Float_t ptmax  = GetHistoPtMax();
-       Float_t phimax = GetHistoPhiMax();
-       Float_t etamax = GetHistoEtaMax();
-       Float_t ptmin  = GetHistoPtMin();
-       Float_t phimin = GetHistoPhiMin();
-       Float_t etamin = GetHistoEtaMin();      
-       
-       //Correlation with charged hadrons
-       if(GetReader()->IsCTSSwitchedOn()) {
-               fhPhiCharged  = new TH2F
-               ("PhiCharged","#phi_{h^{#pm}}  vs p_{T trigger}",
-                nptbins,ptmin,ptmax,nphibins,phimin,phimax); 
-               fhPhiCharged->SetYTitle("#phi_{h^{#pm}} (rad)");
-               fhPhiCharged->SetXTitle("p_{T trigger} (GeV/c)");
-               
-               fhEtaCharged  = new TH2F
-               ("EtaCharged","#eta_{h^{#pm}}  vs p_{T trigger}",
-                nptbins,ptmin,ptmax,netabins,etamin,etamax); 
-               fhEtaCharged->SetYTitle("#eta_{h^{#pm}} (rad)");
-               fhEtaCharged->SetXTitle("p_{T trigger} (GeV/c)");
-               
-               fhDeltaPhiCharged  = new TH2F
-               ("DeltaPhiCharged","#phi_{trigger} - #phi_{h^{#pm}} vs p_{T trigger}",
-               nptbins,ptmin,ptmax,200,0,6.4); 
-               fhDeltaPhiCharged->SetYTitle("#Delta #phi");
-               fhDeltaPhiCharged->SetXTitle("p_{T trigger} (GeV/c)");
-               
-               fhDeltaPhiChargedPt  = new TH2F
-               ("DeltaPhiChargedPt","#phi_{trigger} - #phi_{#p^{#pm}i} vs p_{T h^{#pm}}",
-                nptbins,ptmin,ptmax,200,0,6.4);
-               fhDeltaPhiChargedPt->SetYTitle("#Delta #phi");
-               fhDeltaPhiChargedPt->SetXTitle("p_{T h^{#pm}} (GeV/c)");
-               
-               fhDeltaEtaCharged  = new TH2F
-               ("DeltaEtaCharged","#eta_{trigger} - #eta_{h^{#pm}} vs p_{T trigger}",
-               nptbins,ptmin,ptmax,200,-2,2); 
-               fhDeltaEtaCharged->SetYTitle("#Delta #eta");
-               fhDeltaEtaCharged->SetXTitle("p_{T trigger} (GeV/c)");
-               
-               fhPtImbalanceCharged  = 
-               new TH2F("CorrelationCharged","z_{trigger h^{#pm}} = p_{T h^{#pm}} / p_{T trigger}",
-                               nptbins,ptmin,ptmax,1000,0.,1.2); 
-               fhPtImbalanceCharged->SetYTitle("z_{trigger h^{#pm}}");
-               fhPtImbalanceCharged->SetXTitle("p_{T trigger}");
-               
-               outputContainer->Add(fhPhiCharged) ;
-               outputContainer->Add(fhEtaCharged) ;
-               outputContainer->Add(fhDeltaPhiCharged) ; 
-               outputContainer->Add(fhDeltaEtaCharged) ;
-               outputContainer->Add(fhPtImbalanceCharged) ;
-               outputContainer->Add(fhDeltaPhiChargedPt) ;
-       }  //Correlation with charged hadrons
-       
-       //Correlation with neutral hadrons
-       if(GetReader()->IsEMCALSwitchedOn() || GetReader()->IsPHOSSwitchedOn()){
-               
-               fhPhiNeutral  = new TH2F
-               ("PhiNeutral","#phi_{#pi^{0}}  vs p_{T trigger}",
-               nptbins,ptmin,ptmax,nphibins,phimin,phimax); 
-               fhPhiNeutral->SetYTitle("#phi_{#pi^{0}} (rad)");
-               fhPhiNeutral->SetXTitle("p_{T trigger} (GeV/c)");
-               
-               fhEtaNeutral  = new TH2F
-               ("EtaNeutral","#eta_{#pi^{0}}  vs p_{T trigger}",
-                nptbins,ptmin,ptmax,netabins,etamin,etamax); 
-               fhEtaNeutral->SetYTitle("#eta_{#pi^{0}} (rad)");
-               fhEtaNeutral->SetXTitle("p_{T trigger} (GeV/c)");
-               
-               fhDeltaPhiNeutral  = new TH2F
-               ("DeltaPhiNeutral","#phi_{trigger} - #phi_{#pi^{0}} vs p_{T trigger}",
-               nptbins,ptmin,ptmax,nphibins,phimin,phimax); 
-               fhDeltaPhiNeutral->SetYTitle("#Delta #phi");
-               fhDeltaPhiNeutral->SetXTitle("p_{T trigger} (GeV/c)");
-               
-               fhDeltaPhiNeutralPt  = new TH2F
-               ("DeltaPhiNeutralPt","#phi_{trigger} - #phi_{#pi^{0}} vs p_{T #pi^{0}}}",
-                nptbins,ptmin,ptmax,200,0,6.4); 
-               fhDeltaPhiNeutralPt->SetYTitle("#Delta #phi");
-               fhDeltaPhiNeutralPt->SetXTitle("p_{T trigger} (GeV/c)");
-               
-               fhDeltaEtaNeutral  = new TH2F
-               ("DeltaEtaNeutral","#eta_{trigger} - #eta_{#pi^{0}} vs p_{T trigger}",
-                nptbins,ptmin,ptmax,200,-2,2); 
-               fhDeltaEtaNeutral->SetYTitle("#Delta #eta");
-               fhDeltaEtaNeutral->SetXTitle("p_{T trigger} (GeV/c)");
-               
-               fhPtImbalanceNeutral  = 
-               new TH2F("CorrelationNeutral","z_{trigger #pi} = p_{T #pi^{0}} / p_{T trigger}",
-                                nptbins,ptmin,ptmax,1000,0.,1.2); 
-               fhPtImbalanceNeutral->SetYTitle("z_{trigger #pi^{0}}");
-               fhPtImbalanceNeutral->SetXTitle("p_{T trigger}");
-               
-               outputContainer->Add(fhPhiNeutral) ;  
-               outputContainer->Add(fhEtaNeutral) ;   
-               outputContainer->Add(fhDeltaPhiNeutral) ; 
-               outputContainer->Add(fhDeltaEtaNeutral) ; 
-               outputContainer->Add(fhPtImbalanceNeutral) ;
-               
-               
-               //Keep neutral meson selection histograms if requiered
-               //Setting done in AliNeutralMesonSelection
-               if(GetNeutralMesonSelection()){
-                       TList * nmsHistos = GetNeutralMesonSelection()->GetCreateOutputObjects() ;
-                       if(GetNeutralMesonSelection()->AreNeutralMesonSelectionHistosKept())
-                               for(Int_t i = 0; i < nmsHistos->GetEntries(); i++) outputContainer->Add(nmsHistos->At(i)) ;
-               }
-               
-       }//Correlation with neutral hadrons
-       
-       return outputContainer;
+        outputContainer->Add(fhTrigDeltaPhiDeltaEtaCharged[im]) ;
+        outputContainer->Add(fhTrigCorr[im]);
+        outputContainer->Add(fhTrigUeCorr[im]);
+        
+        }
+    }
+   
+    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}}");
+      
+      fhDeltaPhiDecayCharged  = new TH2F
+      ("DeltaPhiDecayCharged","#phi_{Decay} - #phi_{h^{#pm}} vs p_{T Decay}",
+       nptbins,ptmin,ptmax,140,-2.,5.); 
+      fhDeltaPhiDecayCharged->SetYTitle("#Delta #phi");
+      fhDeltaPhiDecayCharged->SetXTitle("p_{T Decay} (GeV/c)");
+      
+      fhPtImbalanceDecayCharged  = 
+      new TH2F("CorrelationDecayCharged","z_{trigger h^{#pm}} = p_{T h^{#pm}} / p_{T Decay}",
+               nptbins,ptmin,ptmax,200,0.,2.); 
+      fhPtImbalanceDecayCharged->SetYTitle("z_{decay h^{#pm}}");
+      fhPtImbalanceDecayCharged->SetXTitle("p_{T decay}");
+      
+      outputContainer->Add(fhPtPi0DecayRatio) ; 
+      outputContainer->Add(fhDeltaPhiDecayCharged) ; 
+      outputContainer->Add(fhPtImbalanceDecayCharged) ;
+    }    
+    
+    
+    if(fMakeSeveralUE){ 
+      fhDeltaPhiUeLeftCharged  = new TH2F
+      ("DeltaPhiUeLeftChargedPt","#phi_{trigger} - #phi_{#Ueh^{#pm}} vs p_{T Ueh^{#pm}} with UE left side range of trigger particles",
+       nptbins,ptmin,ptmax,140,-2.,5.);
+      fhDeltaPhiUeLeftCharged->SetYTitle("#Delta #phi");
+      fhDeltaPhiUeLeftCharged->SetXTitle("p_{T h^{#pm}} (GeV/c)");
+      outputContainer->Add(fhDeltaPhiUeLeftCharged) ;
+      
+      fhDeltaPhiUeRightCharged  = new TH2F
+      ("DeltaPhiUeRightChargedPt","#phi_{trigger} - #phi_{#Ueh^{#pm}} vs p_{T Ueh^{#pm}} with UE right side range of trigger particles",
+       nptbins,ptmin,ptmax,140,-2.,5.);
+      fhDeltaPhiUeRightCharged->SetYTitle("#Delta #phi");
+      fhDeltaPhiUeRightCharged->SetXTitle("p_{T h^{#pm}} (GeV/c)");
+      outputContainer->Add(fhDeltaPhiUeRightCharged) ;
+      
+      fhPtImbalanceUeLeftCharged  = 
+      new TH2F("CorrelationUeChargedLeft","z_{trigger h^{#pm}} = p_{T Ueh^{#pm}} / p_{T trigger} with UE left side of trigger",
+               nptbins,ptmin,ptmax,200,0.,2.); 
+      fhPtImbalanceUeLeftCharged->SetYTitle("z_{trigger Ueh^{#pm}}");
+      fhPtImbalanceUeLeftCharged->SetXTitle("p_{T trigger}");
+      outputContainer->Add(fhPtImbalanceUeLeftCharged) ;
+      
+      fhPtImbalanceUeRightCharged  = 
+      new TH2F("CorrelationUeChargedRight","z_{trigger h^{#pm}} = p_{T Ueh^{#pm}} / p_{T trigger} with UE right side of trigger",
+               nptbins,ptmin,ptmax,200,0.,2.); 
+      fhPtImbalanceUeRightCharged->SetYTitle("z_{trigger Ueh^{#pm}}");
+      fhPtImbalanceUeRightCharged->SetXTitle("p_{T trigger}");
+      outputContainer->Add(fhPtImbalanceUeRightCharged) ;
+      
+      fhPtHbpUeLeftCharged  = 
+      new TH2F("HbpUeChargedLeft","#xi = ln(1/x_{E}) with charged UE left side of trigger",
+               nptbins,ptmin,ptmax,200,0.,10.); 
+      fhPtHbpUeLeftCharged->SetYTitle("ln(1/x_{E})");
+      fhPtHbpUeLeftCharged->SetXTitle("p_{T trigger}");
+      outputContainer->Add(fhPtHbpUeLeftCharged) ;
+      
+      fhPtHbpUeRightCharged  = 
+      new TH2F("HbpUeChargedRight","#xi = ln(1/x_{E}) with charged UE right side of trigger",
+               nptbins,ptmin,ptmax,200,0.,10.); 
+      fhPtHbpUeRightCharged->SetYTitle("ln(1/x_{E})");
+      fhPtHbpUeRightCharged->SetXTitle("p_{T trigger}");
+      outputContainer->Add(fhPtHbpUeRightCharged) ;
+      
+    }  
+  }  //Correlation with charged hadrons
+  
+  //Correlation with neutral hadrons
+  if(fNeutralCorr){
+    
+    fhDeltaPhiDeltaEtaNeutral  = new TH2F
+    ("DeltaPhiDeltaEtaNeutral","#phi_{trigger} - #phi_{h^{0}} vs #eta_{trigger} - #eta_{h^{0}}",
+     140,-2.,5.,200,-2,2); 
+    fhDeltaPhiDeltaEtaNeutral->SetXTitle("#Delta #phi");
+    fhDeltaPhiDeltaEtaNeutral->SetYTitle("#Delta #eta");   
+         
+    fhPhiNeutral  = new TH2F
+    ("PhiNeutral","#phi_{#pi^{0}}  vs p_{T #pi^{0}}",
+     nptbins,ptmin,ptmax,nphibins,phimin,phimax); 
+    fhPhiNeutral->SetYTitle("#phi_{#pi^{0}} (rad)");
+    fhPhiNeutral->SetXTitle("p_{T #pi^{0}} (GeV/c)");
+    
+    fhEtaNeutral  = new TH2F
+    ("EtaNeutral","#eta_{#pi^{0}}  vs p_{T #pi^{0}}",
+     nptbins,ptmin,ptmax,netabins,etamin,etamax); 
+    fhEtaNeutral->SetYTitle("#eta_{#pi^{0}} (rad)");
+    fhEtaNeutral->SetXTitle("p_{T #pi^{0}} (GeV/c)");
+    
+    fhDeltaPhiNeutral  = new TH2F
+    ("DeltaPhiNeutral","#phi_{trigger} - #phi_{#pi^{0}} vs p_{T trigger}",
+     nptbins,ptmin,ptmax,nphibins,phimin,phimax); 
+    fhDeltaPhiNeutral->SetYTitle("#Delta #phi");
+    fhDeltaPhiNeutral->SetXTitle("p_{T trigger} (GeV/c)");
+    
+    fhDeltaPhiNeutralPt  = new TH2F
+    ("DeltaPhiNeutralPt","#phi_{trigger} - #phi_{#pi^{0}} vs p_{T #pi^{0}}}",
+     nptbins,ptmin,ptmax,140,-2.,5.); 
+    fhDeltaPhiNeutralPt->SetYTitle("#Delta #phi");
+    fhDeltaPhiNeutralPt->SetXTitle("p_{T h^{0}} (GeV/c)");
+    
+    fhDeltaPhiUeNeutralPt  = new TH2F
+    ("DeltaPhiUeNeutralPt","#phi_{trigger} - #phi_{#pi^{0}} vs p_{T #pi^{0}}}",
+     nptbins,ptmin,ptmax,140,-2.,5.); 
+    fhDeltaPhiUeNeutralPt->SetYTitle("#Delta #phi");
+    fhDeltaPhiUeNeutralPt->SetXTitle("p_{T h^{0}} (GeV/c)");
+    
+    fhDeltaEtaNeutral  = new TH2F
+    ("DeltaEtaNeutral","#eta_{trigger} - #eta_{#pi^{0}} vs p_{T trigger}",
+     nptbins,ptmin,ptmax,200,-2,2); 
+    fhDeltaEtaNeutral->SetYTitle("#Delta #eta");
+    fhDeltaEtaNeutral->SetXTitle("p_{T trigger} (GeV/c)");
+    
+    fhPtImbalanceNeutral  = 
+    new TH2F("CorrelationNeutral","z_{trigger #pi} = p_{T #pi^{0}} / p_{T trigger}",
+             nptbins,ptmin,ptmax,200,0.,2.); 
+    fhPtImbalanceNeutral->SetYTitle("z_{trigger #pi^{0}}");
+    fhPtImbalanceNeutral->SetXTitle("p_{T trigger}");
+    
+    fhPtImbalanceUeNeutral  = 
+    new TH2F("CorrelationUeNeutral","z_{trigger #pi} = p_{T #pi^{0}} / p_{T trigger}",
+             nptbins,ptmin,ptmax,200,0.,2.); 
+    fhPtImbalanceUeNeutral->SetYTitle("z_{trigger #pi^{0}}");
+    fhPtImbalanceUeNeutral->SetXTitle("p_{T trigger}");
+    
+    fhPtHbpNeutral  = 
+    new TH2F("HbpNeutral","#xi = ln(1/x_{E}) with neutral particles",
+             nptbins,ptmin,ptmax,200,0.,10.); 
+    fhPtHbpNeutral->SetYTitle("ln(1/x_{E})");
+    fhPtHbpNeutral->SetXTitle("p_{T trigger}");
+    
+    fhPtHbpUeNeutral  = 
+    new TH2F("HbpUeNeutral","#xi = ln(1/x_{E}) with neutral particles",
+             nptbins,ptmin,ptmax,200,0.,10.); 
+    fhPtHbpUeNeutral->SetYTitle("ln(1/x_{E})");
+    fhPtHbpUeNeutral->SetXTitle("p_{T trigger}");
+    
+    
+    outputContainer->Add(fhDeltaPhiDeltaEtaNeutral); 
+    outputContainer->Add(fhPhiNeutral) ;  
+    outputContainer->Add(fhEtaNeutral) ;   
+    outputContainer->Add(fhDeltaPhiNeutral) ; 
+    outputContainer->Add(fhDeltaPhiNeutralPt) ; 
+    outputContainer->Add(fhDeltaPhiUeNeutralPt) ; 
+    outputContainer->Add(fhDeltaEtaNeutral) ; 
+    outputContainer->Add(fhPtImbalanceNeutral) ;
+    outputContainer->Add(fhPtImbalanceUeNeutral) ;  
+    outputContainer->Add(fhPtHbpNeutral) ;
+    outputContainer->Add(fhPtHbpUeNeutral) ;    
+    
+    if(fPi0Trigger){
+      fhDeltaPhiDecayNeutral  = new TH2F
+      ("DeltaPhiDecayNeutral","#phi_{Decay} - #phi_{h^{0}} vs p_{T Decay}",
+       nptbins,ptmin,ptmax,140,-2.,5.); 
+      fhDeltaPhiDecayNeutral->SetYTitle("#Delta #phi");
+      fhDeltaPhiDecayNeutral->SetXTitle("p_{T Decay} (GeV/c)");
+      
+      fhPtImbalanceDecayNeutral  = 
+      new TH2F("CorrelationDecayNeutral","z_{trigger h^{0}} = p_{T h^{0}} / p_{T Decay}",
+               nptbins,ptmin,ptmax,200,0.,2.); 
+      fhPtImbalanceDecayNeutral->SetYTitle("z_{decay h^{0}}");
+      fhPtImbalanceDecayNeutral->SetXTitle("p_{T decay}");
+      
+      outputContainer->Add(fhDeltaPhiDecayNeutral) ; 
+      outputContainer->Add(fhPtImbalanceDecayNeutral) ;
+    }
+    
+    
+    if(fMakeSeveralUE){ 
+      fhDeltaPhiUeLeftNeutral  = new TH2F
+      ("DeltaPhiUeLeftNeutralPt","#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.);
+      fhDeltaPhiUeLeftNeutral->SetYTitle("#Delta #phi");
+      fhDeltaPhiUeLeftNeutral->SetXTitle("p_{T h^{0}} (GeV/c)");
+      outputContainer->Add(fhDeltaPhiUeLeftNeutral) ;
+      
+      fhDeltaPhiUeRightNeutral  = new TH2F
+      ("DeltaPhiUeRightNeutralPt","#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.);
+      fhDeltaPhiUeRightNeutral->SetYTitle("#Delta #phi");
+      fhDeltaPhiUeRightNeutral->SetXTitle("p_{T h^{0}} (GeV/c)");
+      outputContainer->Add(fhDeltaPhiUeRightNeutral) ;
+      
+      fhPtImbalanceUeLeftNeutral  = 
+      new TH2F("CorrelationUeNeutralLeft","z_{trigger h^{0}} = p_{T Ueh^{0}} / p_{T trigger} with neutral UE left side of trigger",
+               nptbins,ptmin,ptmax,140,0.,2.); 
+      fhPtImbalanceUeLeftNeutral->SetYTitle("z_{trigger Ueh^{0}}");
+      fhPtImbalanceUeLeftNeutral->SetXTitle("p_{T trigger}");
+      outputContainer->Add(fhPtImbalanceUeLeftNeutral) ;
+      
+      fhPtImbalanceUeRightNeutral  = 
+      new TH2F("CorrelationUeNeutralRight","z_{trigger h^{0}} = p_{T Ueh^{0}} / p_{T trigger} with neutral UE right side of trigger",
+               nptbins,ptmin,ptmax,200,0.,2.); 
+      fhPtImbalanceUeRightNeutral->SetYTitle("z_{trigger Ueh^{0}}");
+      fhPtImbalanceUeRightNeutral->SetXTitle("p_{T trigger}");
+      outputContainer->Add(fhPtImbalanceUeRightNeutral) ;
+      
+      fhPtHbpUeLeftNeutral  = 
+      new TH2F("HbpUeNeutralLeft","#xi = ln(1/x_{E}) with neutral UE left side of trigger",
+               nptbins,ptmin,ptmax,200,0.,10.); 
+      fhPtHbpUeLeftNeutral->SetYTitle("ln(1/x_{E})");
+      fhPtHbpUeLeftNeutral->SetXTitle("p_{T trigger}");
+      outputContainer->Add(fhPtHbpUeLeftNeutral) ;
+      
+      fhPtHbpUeRightNeutral  = 
+      new TH2F("HbpUeNeutralRight","#xi = ln(1/x_{E}) with neutral UE right side of trigger",
+               nptbins,ptmin,ptmax,200,0.,10.); 
+      fhPtHbpUeRightNeutral->SetYTitle("ln(1/x_{E})");
+      fhPtHbpUeRightNeutral->SetXTitle("p_{T trigger}");
+      outputContainer->Add(fhPtHbpUeRightNeutral) ;
+      
+    }  
+    
+    //Keep neutral meson selection histograms if requiered
+    //Setting done in AliNeutralMesonSelection
+    if(GetNeutralMesonSelection()){
+      TList * nmsHistos = GetNeutralMesonSelection()->GetCreateOutputObjects() ;
+      if(GetNeutralMesonSelection()->AreNeutralMesonSelectionHistosKept())
+        for(Int_t i = 0; i < nmsHistos->GetEntries(); i++) outputContainer->Add(nmsHistos->At(i)) ;
+      delete nmsHistos;
+    }
+    
+  }//Correlation with neutral hadrons
+  
+  return outputContainer;
+  
 }
 
- //____________________________________________________________________________
+//____________________________________________________________________________
 void AliAnaParticleHadronCorrelation::InitParameters()
 {
+  
   //Initialize the parameters of the analysis.
-  SetInputAODName("photons");
-  SetPtCutRange(2,300);
+  SetInputAODName("PWG4Particle");
+  SetAODObjArrayName("Hadrons");  
+  AddToHistogramsName("AnaHadronCorr_");
+
+  SetPtCutRange(0.,300);
   fDeltaPhiMinCut = 1.5 ;
   fDeltaPhiMaxCut = 4.5 ;
+  fSelectIsolated = kFALSE;
+  fMakeSeveralUE = kFALSE;
+  fUeDeltaPhiMinCut = 1. ;
+  fUeDeltaPhiMaxCut = 1.5 ;
+  fNeutralCorr = kFALSE ;
+  fPi0Trigger = kFALSE ;
 
 }
 
@@ -241,302 +540,902 @@ void AliAnaParticleHadronCorrelation::Print(const Option_t * opt) const
   if(! opt)
     return;
   
-  Info("*** Print *** ", "%s %s", GetName(), GetTitle() ) ;
-  printf("pT Hadron       >    %2.2f\n", GetMinPt()) ; 
-  printf("pT Hadron       <    %2.2f\n", GetMaxPt()) ; 
+  printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
+  AliAnaPartCorrBaseClass::Print(" ");
+
   printf("Phi trigger particle-Hadron      <     %3.2f\n", fDeltaPhiMaxCut) ; 
   printf("Phi trigger particle-Hadron      >     %3.2f\n", fDeltaPhiMinCut) ;
+  printf("Isolated Trigger?  %d\n", fSelectIsolated) ;
+  printf("Phi trigger particle-UeHadron    <    %3.2f\n", fUeDeltaPhiMaxCut) ; 
+  printf("Phi trigger particle-UeHadron    >    %3.2f\n", fUeDeltaPhiMinCut) ;
+  printf("Several UE?  %d\n", fMakeSeveralUE) ;
+  printf("Name of AOD Pi0 Branch %s \n",fPi0AODBranchName.Data());
+  printf("Do Decay-hadron correlation ?  %d\n", fPi0Trigger) ;
+
+  
+} 
+
+//__________________________________________________________________
+TObjString* AliAnaParticleHadronCorrelation::GetAnalysisCuts()
+{
+  //Save parameters used for analysis
+  TString parList ; //this will be list of parameters used for this analysis.
+  const Int_t buffersize = 255;
+  char onePar[buffersize] ;
+  snprintf(onePar,buffersize,"--- AliAnaPaticleHadronCorrelation ---\n") ;
+  parList+=onePar ;    
+  snprintf(onePar,buffersize,"Phi trigger particle-Hadron      <     %3.2f ", fDeltaPhiMaxCut) ; 
+  parList+=onePar ;
+  snprintf(onePar,buffersize,"Phi trigger particle-Hadron      >     %3.2f ", fDeltaPhiMinCut) ;
+  parList+=onePar ;
+  snprintf(onePar,buffersize,"Isolated Trigger?  %d\n", fSelectIsolated) ;
+  parList+=onePar ;
+  snprintf(onePar,buffersize,"Phi trigger particle-UeHadron    <    %3.2f ", fUeDeltaPhiMaxCut) ; 
+  parList+=onePar ;
+  snprintf(onePar,buffersize,"Phi trigger particle-UeHadron    >    %3.2f ", fUeDeltaPhiMinCut) ;
+  parList+=onePar ;
+  snprintf(onePar,buffersize,"Several UE?  %d\n", fMakeSeveralUE) ;
+  parList+=onePar ;
+  snprintf(onePar,buffersize,"Name of AOD Pi0 Branch %s ",fPi0AODBranchName.Data());
+  parList+=onePar ;
+  snprintf(onePar,buffersize,"Do Decay-hadron correlation ?  %d", fPi0Trigger) ;
+  parList+=onePar ;
+
+  //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() 
+  
+  return new TObjString(parList) ;  
 
 } 
 
+
 //____________________________________________________________________________
 void  AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD()  
 {  
-       //Particle-Hadron Correlation Analysis, fill AODs
-       
-       if(!GetInputAODBranch())
-               AliFatal(Form("ParticleHadronCorrelation::FillAOD: No input particles in AOD with name branch < %s > \n",GetInputAODName().Data()));
-       
-       if(GetDebug() > 1){
-               printf("Begin hadron correlation analysis, fill AODs \n");
-               printf("In particle branch aod entries %d\n", GetInputAODBranch()->GetEntriesFast());
-               printf("In CTS aod entries %d\n", GetAODCTS()->GetEntriesFast());
-               printf("In EMCAL aod entries %d\n", GetAODEMCAL()->GetEntriesFast());
-               printf("In PHOS aod entries %d\n", GetAODPHOS()->GetEntriesFast());
-       }
+  //Particle-Hadron Correlation Analysis, fill AODs
+  
+  if(!GetInputAODBranch()){
+    printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - No input particles in AOD with name branch < %s >, STOP \n",GetInputAODName().Data());
+    abort();
+  }
        
-       //Loop on stored AOD particles, trigger
-       Int_t naod = GetInputAODBranch()->GetEntriesFast();
-       for(Int_t iaod = 0; iaod < naod ; iaod++){
-               AliAODPWG4ParticleCorrelation* particle =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
-               //Make correlation with charged hadrons
-               if(GetReader()->IsCTSSwitchedOn() )
-                       MakeChargedCorrelation(particle, (TSeqCollection*)GetAODCTS(),kFALSE);
-               
-               //Make correlation with neutral pions
-               //Trigger particle in PHOS, correlation with EMCAL
-               if(particle->GetDetector()=="PHOS" && GetReader()->IsEMCALSwitchedOn() && GetAODEMCAL()->GetEntriesFast() > 0)
-                       MakeNeutralCorrelation(particle,(TSeqCollection*)GetAODEMCAL(),kFALSE);
-               //Trigger particle in EMCAL, correlation with PHOS
-               else if(particle->GetDetector()=="EMCAL" && GetReader()->IsPHOSSwitchedOn() && GetAODPHOS()->GetEntriesFast() > 0)
-                       MakeNeutralCorrelation(particle,(TSeqCollection*)GetAODPHOS(),kFALSE);
-               //Trigger particle in CTS, correlation with PHOS, EMCAL and CTS
-               else if(particle->GetDetector()=="CTS" ){
-                       if(GetReader()->IsPHOSSwitchedOn() && GetAODPHOS()->GetEntriesFast() > 0) 
-                               MakeNeutralCorrelation(particle,(TSeqCollection*)GetAODPHOS(),kFALSE);
-                       if(GetReader()->IsEMCALSwitchedOn() && GetAODEMCAL()->GetEntriesFast() > 0) 
-                               MakeNeutralCorrelation(particle,(TSeqCollection*)GetAODEMCAL(),kFALSE);
-               }
-               
-       }//Aod branch loop
+  if(strcmp(GetInputAODBranch()->GetClass()->GetName(), "AliAODPWG4ParticleCorrelation")){
+    printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - Wrong type of AOD object, change AOD class name in input AOD: It should be <AliAODPWG4ParticleCorrelation> and not <%s> \n",GetInputAODBranch()->GetClass()->GetName());
+    abort();
+  }
        
-       if(GetDebug() > 1) printf("End hadron correlation analysis, fill AODs \n");
+  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", GetAODCTS()->GetEntriesFast());
+    printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - In EMCAL aod entries %d\n", GetAODEMCAL()->GetEntriesFast());
+    printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - In PHOS aod entries %d\n", GetAODPHOS()->GetEntriesFast());
+  }
+  
+  //Loop on stored AOD particles, trigger
+  Double_t ptTrig    = 0.;
+  Int_t    trigIndex = -1;
+  Int_t naod = GetInputAODBranch()->GetEntriesFast();
+  //fhNclustersNtracks->Fill(naod, GetAODCTS()->GetEntriesFast());
+  for(Int_t iaod = 0; iaod < naod ; iaod++){
+    AliAODPWG4ParticleCorrelation* particle =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
+    //find the leading particles with highest momentum
+    if (particle->Pt()>ptTrig) {
+      ptTrig = particle->Pt() ;
+      trigIndex = iaod ;
+    }
+  }//Aod branch loop
        
+  //Do correlation with leading particle
+  if(trigIndex!=-1){
+         
+    AliAODPWG4ParticleCorrelation* particle =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(trigIndex));
+    //Make correlation with charged hadrons
+    if(GetReader()->IsCTSSwitchedOn() )
+      MakeChargedCorrelation(particle, GetAODCTS(),kFALSE);
+    
+    TObjArray * pi0list = (TObjArray*) GetAODBranch(fPi0AODBranchName); //For the future, foresee more possible pi0 lists
+    if(fNeutralCorr && pi0list && pi0list->GetEntriesFast() > 0)
+      MakeNeutralCorrelation(particle, pi0list,kFALSE);
+    
+  }//Correlate leading
+  
+  if(GetDebug() > 1) printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - End fill AODs \n");
+  
 }
 
 //____________________________________________________________________________
 void  AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms()  
 {  
-       //Particle-Hadron Correlation Analysis, fill histograms
-       
-       if(!GetInputAODBranch())
-               AliFatal(Form("ParticleHadronCorrelation::FillHistos: No input particles in AOD with name branch < %s > \n",GetInputAODName().Data()));
-       
-       if(GetDebug() > 1){
-               printf("Begin hadron correlation analysis, fill histograms \n");
-               printf("In particle branch aod entries %d\n", GetInputAODBranch()->GetEntriesFast());
-       }
-       
-       //Loop on stored AOD particles
-       Int_t naod = GetInputAODBranch()->GetEntriesFast();
-       for(Int_t iaod = 0; iaod < naod ; iaod++){
-               AliAODPWG4ParticleCorrelation* particle =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
-               
-               if(GetDebug() > 1){
-                       printf("Particle %d, In Track Refs  entries %d\n", iaod, (particle->GetRefTracks())->GetEntriesFast());
-                       printf("Particle %d, In Cluster Refs entries %d\n",iaod, (particle->GetRefClusters())->GetEntriesFast());
-               }
-               
-               //Make correlation with charged hadrons
-               if((particle->GetRefTracks())->GetEntriesFast() > 0)
-                       MakeChargedCorrelation(particle, (TSeqCollection*) (particle->GetRefTracks()),kTRUE);
-               
-               //Make correlation with neutral pions
-               if((particle->GetRefClusters())->GetEntriesFast() > 0)
-                       MakeNeutralCorrelation(particle,  (TSeqCollection*) (particle->GetRefClusters()), kTRUE);
-               
-       }//Aod branch loop
-       
-       if(GetDebug() > 1) printf("End hadron correlation analysis, fill histograms \n");
-       
+  //Particle-Hadron Correlation Analysis, fill histograms
+  
+  if(!GetInputAODBranch()){
+    printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms() - No input particles in AOD with name branch < %s >, STOP \n",GetInputAODName().Data());
+    abort();
+  }
+  
+  if(GetDebug() > 1){
+    printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms() - Begin hadron correlation analysis, fill histograms \n");
+    printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms() - In particle branch aod entries %d\n", GetInputAODBranch()->GetEntriesFast());
+  }
+  
+  //Get the vertex and check it is not too large in z
+  Double_t v[3] = {0,0,0}; //vertex ;
+  GetReader()->GetVertex(v);
+  if(!GetMixedEvent() && TMath::Abs(v[2]) > GetZvertexCut()) return ;  
+  
+  //Loop on stored AOD particles, find leading
+  Int_t naod = GetInputAODBranch()->GetEntriesFast();
+ // if(naod!=0)fhVertex->Fill(v[0],v[1],v[2]);
+  Double_t ptTrig = 0.;
+  Int_t trigIndex = -1 ;
+  for(Int_t iaod = 0; iaod < naod ; iaod++){    //loop on input trigger AOD file 
+    AliAODPWG4ParticleCorrelation* particle =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
+    //vertex cut in case of mixing
+    if (GetMixedEvent()) {
+      Int_t evt=-1;
+      Int_t id =-1;
+      if     (particle->GetCaloLabel(0)  >= 0 ){
+        id=particle->GetCaloLabel(0); 
+        if(id >= 0 )evt= GetMixedEvent()-> EventIndexForCaloCluster(id) ;
+      }
+      else if(particle->GetTrackLabel(0) >= 0 ){
+        id=particle->GetTrackLabel(0);
+        if(id >= 0 )evt= GetMixedEvent()->EventIndex(id) ;
+      }
+      else continue;
+      
+      if (TMath::Abs(GetVertex(evt)[2]) > GetZvertexCut()) 
+        return ;
+    }
+
+    //check if the particle is isolated or if we want to take the isolation into account
+    if(OnlyIsolated() && !particle->IsIsolated()) continue;
+    //check if inside the vertex cut
+    //find the leading particles with highest momentum
+    if (particle->Pt()>ptTrig) {
+      ptTrig = particle->Pt() ;
+                 trigIndex = iaod ;
+    }
+  }//finish searching for leading trigger particle
+  if(trigIndex!=-1){ //using trigger partilce to do correlations
+    AliAODPWG4ParticleCorrelation* particle =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(trigIndex));
+         
+    if (GetMixedEvent()) {
+      Int_t evt=-1;
+      Int_t id = 0;
+      if     (particle->GetCaloLabel(0)  >= 0 ){
+        id=particle->GetCaloLabel(0); 
+        if(id >= 0 )evt= GetMixedEvent()-> EventIndexForCaloCluster(id) ;
+      }
+      else if(particle->GetTrackLabel(0) >= 0 ){
+        id=particle->GetTrackLabel(0);
+        if(id >= 0 )evt= GetMixedEvent()->EventIndex(id) ;
+      }
+      else return;
+      
+      if (TMath::Abs(GetVertex(evt)[2]) > GetZvertexCut()) 
+        return ;
+    }
+        
+    //Fill leading particle histogram   
+    fhPtLeading->Fill(particle->Pt());
+    Float_t phi = particle->Phi();
+    if(phi<0)phi+=TMath::TwoPi();
+    fhPhiLeading->Fill(particle->Pt(), phi);
+    fhEtaLeading->Fill(particle->Pt(), particle->Eta());
+    
+    //Make correlation with charged hadrons
+    if(GetReader()->IsCTSSwitchedOn() )
+      MakeChargedCorrelation(particle, GetAODCTS(),kTRUE);
+    
+    TObjArray * pi0list = (TObjArray*) GetAODBranch(fPi0AODBranchName); //For the future, foresee more possible pi0 lists
+    if(fNeutralCorr && pi0list && pi0list->GetEntriesFast() > 0)
+      MakeNeutralCorrelation(particle, pi0list,kTRUE);
+    
+  }//Aod branch loop
+  
+  if(GetDebug() > 1) printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms() - End fill histograms \n");
+  
 }
 
 //____________________________________________________________________________
-void  AliAnaParticleHadronCorrelation::MakeChargedCorrelation(AliAODPWG4ParticleCorrelation *aodParticle, TSeqCollection* pl, const Bool_t bFillHisto)
+void  AliAnaParticleHadronCorrelation::MakeChargedCorrelation(AliAODPWG4ParticleCorrelation *aodParticle, TObjArray* const pl, const Bool_t bFillHisto)
 {  
   // Charged Hadron Correlation Analysis
-  if(GetDebug() > 1)printf("Make trigger particle - charged hadron correlation \n");
+  if(GetDebug() > 1)printf("AliAnaParticleHadronCorrelation::MakeChargedCorrelation() - Make trigger particle - charged hadron correlation \n");
+  
+  Int_t evtIndex11 = 0 ; 
+  Int_t evtIndex12 = 0 ;
+  Int_t indexPhoton1 = -1 ;
+  Int_t indexPhoton2 = -1 ;  
+  Double_t v[3] = {0,0,0}; //vertex ;
+  GetReader()->GetVertex(v);
+  if(!GetMixedEvent() && TMath::Abs(v[2]) > GetZvertexCut()) return ;  
+
+  Int_t nTracks = GetAODCTS()->GetEntriesFast() ;
+  
+  if (GetMixedEvent()) {
+    evtIndex11 = GetMixedEvent()->EventIndexForCaloCluster(aodParticle->GetCaloLabel(0)) ;
+    evtIndex12 = GetMixedEvent()->EventIndexForCaloCluster(aodParticle->GetCaloLabel(1)) ;    
+  }
   
   Double_t ptTrig  = aodParticle->Pt();
+  Double_t pxTrig  = aodParticle->Px();
+  Double_t pyTrig  = aodParticle->Py();
+  
   Double_t phiTrig = aodParticle->Phi();
+  Double_t etaTrig = aodParticle->Eta();
+  
   Double_t pt   = -100.;
+  Double_t px   = -100.;
+  Double_t py   = -100.;
   Double_t rat  = -100.; 
+  Double_t xE   = -100.; 
+  Double_t cosi = -100.; 
   Double_t phi  = -100. ;
   Double_t eta  = -100. ;
-  Double_t p[3];
+  TVector3 p3;  
+       
+  TObjArray * reftracks = 0x0;
+  Int_t nrefs = 0;
+       
+  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.; 
+  Float_t deltaphi    = -100. ;
+  Float_t deltaphiDecay1 = -100. ;
+  Float_t deltaphiDecay2 = -100. ;
+  TObjArray * clusters   = 0x0 ;  
+  TLorentzVector photonMom ;   
+  
+  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 = GetAODEMCAL() ;
+      else  clusters = GetAODPHOS() ;
+      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
   
   //Track loop, select tracks with good pt, phi and fill AODs or histograms
-  for(Int_t ipr = 0;ipr < pl->GetEntries() ; ipr ++ ){
+  //Int_t currentIndex = -1 ; 
+  for(Int_t ipr = 0;ipr < pl->GetEntriesFast() ; ipr ++ ){
     AliAODTrack * track = (AliAODTrack *) (pl->At(ipr)) ;
-    track->GetPxPyPz(p) ;
-    TLorentzVector mom(p[0],p[1],p[2],0);
-    pt    = mom.Pt();
-    eta  = mom.Eta();
-    phi  = mom.Phi() ;
-    if(phi<0) phi+=TMath::TwoPi();
-    rat   = pt/ptTrig ;
 
-    if(IsFidutialCutOn()){
-      Bool_t in = GetFidutialCut()->IsInFidutialCut(mom,"CTS") ;
-      if(! in ) continue ;
-    }    
+    //check if inside the vertex cut
+    //printf("charge = %d\n", track->Charge());
+    Int_t evtIndex2 = 0 ; 
+    if (GetMixedEvent()) {
+      evtIndex2 = GetMixedEvent()->EventIndex(track->GetID()) ;
+      if (evtIndex11 == evtIndex2 || evtIndex12 == evtIndex2) // photon and track from different events
+        continue ; 
+      //vertex cut
+      if (TMath::Abs(GetVertex(evtIndex2)[2]) > GetZvertexCut()) 
+         return ;
+      //      if(currentIndex == evtIndex2) // tracks from different event 
+      //        continue ;
+      //      currentIndex = evtIndex2 ;
+    }
+    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();
 
     //Select only hadrons in pt range
     if(pt < GetMinPt() || pt > GetMaxPt()) continue ;
+    //remove trigger itself for correlation when use charged triggers    
+    if(track->GetID()==aodParticle->GetTrackLabel(0) && pt==ptTrig && phi==phiTrig && eta==etaTrig) 
+      continue ;    
+    if(IsFiducialCutOn()){
+      Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,"CTS") ;
+      if(! in ) continue ;
+    }    
+    //jumped out this event if near side associated partile pt larger than trigger
+    if(pt > ptTrig && TMath::Abs(phi-phiTrig)<TMath::PiOver2())  break ;
+    rat   = pt/ptTrig ;
+    xE    = -(px*pxTrig+py*pyTrig)/(ptTrig*ptTrig);
+    if(xE <0.)xE =-xE;
+    cosi = TMath::Log(1/xE);
+    // printf("rat = %f, xE = %f, cosi =%f \n", rat, xE, cosi);
+    // printf("phi = %f \n", phi);
     
-    if(GetDebug() > 2)
-      printf("charged hadron: pt %f, phi %f, phi trigger %f. Cuts:  delta phi min %2.2f,  max%2.2f, pT min %2.2f \n",
-            pt,phi,phiTrig,fDeltaPhiMinCut, fDeltaPhiMaxCut, GetMinPt());
+     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;
+    if(deltaphi< -TMath::PiOver2()) deltaphi+=TMath::TwoPi();
+    if(deltaphi>3*TMath::PiOver2()) deltaphi-=TMath::TwoPi();
 
+    Double_t pout = pt*TMath::Sin(deltaphi) ;
+    
+    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 min %2.2f \n",
+             pt,phi, phiTrig,fDeltaPhiMinCut, deltaphi, fDeltaPhiMaxCut, GetMinPt());
+    
     if(bFillHisto){
       // Fill Histograms
-      fhEtaCharged->Fill(ptTrig,eta);
-      fhPhiCharged->Fill(ptTrig,phi);
+      fhEtaCharged->Fill(pt,eta);
+      fhPhiCharged->Fill(pt,phi);
       fhDeltaEtaCharged->Fill(ptTrig,aodParticle->Eta()-eta);
-      fhDeltaPhiCharged->Fill(ptTrig,phiTrig-phi);
-      fhDeltaPhiChargedPt->Fill(pt,phiTrig-phi);
-      //Selection within angular range
-      if(((phiTrig-phi)> fDeltaPhiMinCut) && ((phiTrig-phi)<fDeltaPhiMaxCut) ){
-       if(GetDebug() > 2 ) printf("Selected charge for momentum imbalance: pt %2.2f, phi %2.2f, eta %2.2f ",pt,phi,eta);
-       fhPtImbalanceCharged->Fill(ptTrig,rat);
-      } 
-    }
-    else{
-      //Fill AODs
-      aodParticle->AddTrack(track);
+      fhDeltaPhiCharged->Fill(ptTrig, deltaphi);
+      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))fhTrigDeltaPhiDeltaEtaCharged[im]->Fill(ptTrig,deltaphi,aodParticle->Eta()-eta);
+        }
+      }
+      //delta phi cut for correlation
+      if( (deltaphi > fDeltaPhiMinCut) && ( deltaphi < fDeltaPhiMaxCut) ) {
+        fhDeltaPhiChargedPt->Fill(pt,deltaphi);
+        fhPtImbalanceCharged->Fill(ptTrig,xE); 
+        fhPtHbpCharged->Fill(ptTrig,cosi);
+        fhPoutPtTrigPtAssoc->Fill(ptTrig, pt, pout) ;
+        fhPtTrigCharged->Fill(ptTrig, pt) ;
+        if(track->Charge()>0) fhPtImbalancePosCharged->Fill(ptTrig,xE) ;
+        else fhPtImbalanceNegCharged->Fill(ptTrig,xE) ;
+        //fill different multiplicity histogram
+        if(DoEventSelect()){
+          for(Int_t im=0; im<GetMultiBin(); im++){
+            if(nTracks < ( GetMaxMulti() - GetMinMulti() )/GetMultiBin()*(im+1))
+              fhTrigCorr[im]->Fill(ptTrig,xE);
+          }
+        } //multiplicity events selection
+      } //delta phi cut for correlation
+      else { //UE study
+        fhDeltaPhiUeChargedPt->Fill(pt,deltaphi);
+        fhUePoutPtTrigPtAssoc->Fill(ptTrig, pt, pout) ;
+        fhPtImbalanceUeCharged->Fill(ptTrig,xE);
+        fhPtHbpUeCharged->Fill(ptTrig,cosi);
+        if(DoEventSelect()){
+          for(Int_t im=0; im<GetMultiBin(); im++){
+            if(nTracks < ( GetMaxMulti() - GetMinMulti() )/GetMultiBin()*(im+1))
+              fhTrigUeCorr[im]->Fill(ptTrig,xE);
+          }
+        } //multiplicity events selection
+        
+      } //UE study
       
+      if(fPi0Trigger){
+        if(indexPhoton1!=-1 && indexPhoton2!=-1){
+          fhDeltaPhiDecayCharged->Fill(ptDecay1, deltaphiDecay1);
+          fhDeltaPhiDecayCharged->Fill(ptDecay2, deltaphiDecay2);
+          if(GetDebug() > 1)printf("deltaPhoton1 = %f, deltaPhoton2 = %f \n", deltaphiDecay1, deltaphiDecay2);
+          if( (deltaphiDecay1 > fDeltaPhiMinCut) && ( deltaphiDecay1 < fDeltaPhiMaxCut) )
+            fhPtImbalanceDecayCharged->Fill(ptDecay1,ratDecay1); 
+          if( (deltaphiDecay2 > fDeltaPhiMinCut) && ( deltaphiDecay2 < fDeltaPhiMaxCut) )
+            fhPtImbalanceDecayCharged->Fill(ptDecay2,ratDecay2);
+          if(GetDebug() > 1)printf("ratPhoton1 = %f, ratPhoton2 = %f \n", pt/ptDecay1, pt/ptDecay2);
+        } //index of decay photons found
+      } //make decay-hadron correlation          
+      
+      //several UE calculation 
+      if(fMakeSeveralUE){
+        if((deltaphi<-fUeDeltaPhiMinCut) && (deltaphi >-fUeDeltaPhiMaxCut)){  
+          fhDeltaPhiUeLeftCharged->Fill(pt,deltaphi);
+          fhPtImbalanceUeLeftCharged->Fill(ptTrig,rat);
+          fhPtHbpUeLeftCharged->Fill(ptTrig,cosi);
+        }
+        if((deltaphi>fUeDeltaPhiMinCut) && (deltaphi <fUeDeltaPhiMaxCut)){  
+          fhDeltaPhiUeRightCharged->Fill(pt,deltaphi);
+          fhPtImbalanceUeRightCharged->Fill(ptTrig,rat);
+          fhPtHbpUeRightCharged->Fill(ptTrig,cosi);
+          
+        }
+      } //several UE calculation
+      
+    } //Fill histogram 
+    else{
+      nrefs++;
+      if(nrefs==1){
+        reftracks = new TObjArray(0);
+        reftracks->SetName(GetAODObjArrayName()+"Tracks");
+        reftracks->SetOwner(kFALSE);
+      }
+      reftracks->Add(track);
     }//aod particle loop
   }// track loop
-
+  
+  //Fill AOD with reference tracks, if not filling histograms
+  if(!bFillHisto && reftracks) {
+    aodParticle->AddObjArray(reftracks);
+  }
+  
+  //delete reftracks;
+  
 }  
+
+
+
 //____________________________________________________________________________
-void  AliAnaParticleHadronCorrelation::MakeNeutralCorrelation(AliAODPWG4ParticleCorrelation * aodParticle,TSeqCollection* pl, const Bool_t bFillHisto)  
+//void  AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillAOD(AliAODPWG4ParticleCorrelation* const aodParticle,TObjArray* const pl, TString detector)  
+//{  
+//    // Neutral Pion Correlation Analysis, find pi0, put them in new output aod, if correlation cuts passed
+//  if(GetDebug() > 1) printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillAOD() - Make trigger particle - neutral hadron correlation \n");
+//  
+//  if(!NewOutputAOD()){
+//    printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillAOD() - Output aod not created, set AOD class name and branch name in the configuration file, STOP! \n");
+//    abort();
+//  }
+//  
+//  Double_t phiTrig = aodParticle->Phi();
+//  Int_t      tag = 0;
+//  TLorentzVector gammai;
+//  TLorentzVector gammaj;
+//  
+//  //Get vertex for photon momentum calculation
+//  
+//  if(GetReader()->GetDataType()!= AliCaloTrackReader::kMC) 
+//  {
+//    for (Int_t iev = 0; iev < GetNMixedEvent(); iev++) {
+//      if (!GetMixedEvent()) 
+//        GetReader()->GetVertex(GetVertex(iev));
+//      else 
+//        GetMixedEvent()->GetVertexOfEvent(iev)->GetXYZ(GetVertex(iev)); 
+//    }
+//  }
+//  Double_t vertex2[] = {0.0,0.0,0.0} ; //vertex of second input aod
+//  if(GetReader()->GetDataType()!= AliCaloTrackReader::kMC) 
+//  {
+//     if(GetReader()->GetSecondInputAODTree()) GetReader()->GetSecondInputAODVertex(vertex2);
+//  }
+//     
+//    //Cluster loop, select pairs with good pt, phi and fill AODs or histograms
+//    //Int_t iEvent= GetReader()->GetEventNumber() ;
+//  Int_t nclus = pl->GetEntriesFast();
+//  for(Int_t iclus = 0;iclus < nclus ; iclus ++ ){
+//    AliVCluster * calo = (AliVCluster *) (pl->At(iclus)) ;
+//    
+//    Int_t evtIndex1 = 0 ; 
+//    if (GetMixedEvent()) {
+//      evtIndex1=GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ; 
+//    }
+//    
+//
+//      //Input from second AOD?
+//    Int_t inputi = 0;
+//    if     (aodParticle->GetDetector() == "EMCAL" && GetReader()->GetAODEMCALNormalInputEntries() <= iclus) 
+//      inputi = 1 ;
+//    else if(aodParticle->GetDetector() == "PHOS"  && GetReader()->GetAODPHOSNormalInputEntries()  <= iclus) 
+//      inputi = 1;
+//    
+//      //Cluster selection, not charged, with photon or pi0 id and in fiducial cut
+//    //FIXME
+//    Int_t pdg=0;
+//    //if     (inputi == 0 && !SelectCluster(calo, GetVertex(evtIndex1),  gammai, pdg))  
+//      continue ;
+//    //MEFIX
+//    else if(inputi == 1 && !SelectCluster(calo, vertex2, gammai, pdg))  
+//      continue ;
+//    
+//    if(GetDebug() > 2)
+//      printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillAOD() - Neutral cluster in %s: pt %f, phi %f, phi trigger %f. Cuts:  delta phi min %2.2f,  max %2.2f, pT min %2.2f \n",
+//             detector.Data(), gammai.Pt(),gammai.Phi(),phiTrig,fDeltaPhiMinCut, fDeltaPhiMaxCut, GetMinPt());
+//    
+//      //2 gamma overlapped, found with PID
+//    if(pdg == AliCaloPID::kPi0){ 
+//      
+//        //Select only hadrons in pt range
+//      if(gammai.Pt() < GetMinPt() || gammai.Pt() > GetMaxPt()) 
+//        continue ;
+//      
+//        //Selection within angular range
+//      Float_t phi = gammai.Phi();
+//      if(phi < 0) phi+=TMath::TwoPi();
+//        //Float_t deltaphi = TMath::Abs(phiTrig-phi);
+//        //if( (deltaphi < fDeltaPhiMinCut) || ( deltaphi > fDeltaPhiMaxCut) ) continue ;
+//      
+//      AliAODPWG4Particle pi0 = AliAODPWG4Particle(gammai);
+//        //pi0.SetLabel(calo->GetLabel());
+//      pi0.SetPdg(AliCaloPID::kPi0);
+//      pi0.SetDetector(detector);
+//      
+//      if(IsDataMC()){
+//        pi0.SetTag(GetMCAnalysisUtils()->CheckOrigin(calo->GetLabel(),GetReader(),inputi));
+//        if(GetDebug() > 0) printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillAOD() - Origin of candidate %d\n",pi0.GetTag());
+//      }//Work with stack also 
+//       //Set the indeces of the original caloclusters  
+//      pi0.SetCaloLabel(calo->GetID(),-1);
+//      AddAODParticle(pi0);
+//      
+//      if(GetDebug() > 2) 
+//        printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillAOD() - Correlated with selected pi0 (pid): pt %f, phi %f\n",pi0.Pt(),pi0.Phi());
+//      
+//    }// pdg = 111
+//    
+//      //Make invariant mass analysis
+//    else if(pdg == AliCaloPID::kPhoton){     
+//        //Search the photon companion in case it comes from  a Pi0 decay
+//        //Apply several cuts to select the good pair;
+//      for(Int_t jclus = iclus+1; jclus < pl->GetEntries() ; jclus ++ ){
+//        AliVCluster * calo2 = (AliVCluster *) (pl->At(jclus)) ;
+//        Int_t evtIndex2 = 0 ; 
+//        if (GetMixedEvent()) {
+//          evtIndex2=GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ; 
+//        }
+//        if (GetMixedEvent() && (evtIndex1 == evtIndex2))
+//          continue ;
+//        
+//          //Input from second AOD?
+//        Int_t inputj = 0;
+//        if     (aodParticle->GetDetector() == "EMCAL" && GetReader()->GetAODEMCALNormalInputEntries() <= jclus) 
+//          inputj = 1;
+//        else if(aodParticle->GetDetector() == "PHOS"  && GetReader()->GetAODPHOSNormalInputEntries()  <= jclus) 
+//          inputj = 1;
+//        
+//          //Cluster selection, not charged with photon or pi0 id and in fiducial cut
+//        Int_t pdgj=0;
+//        //FIXME
+//        //if     (inputj == 0 && !SelectCluster(calo2, GetVertex(evtIndex2),  gammaj, pdgj))  
+//          continue ;
+//        //MEFIX
+//
+//        else if(inputj == 1 && !SelectCluster(calo2, vertex2, gammaj, pdgj))  
+//          continue ;
+//        //FIXME
+//        //if(!SelectCluster(calo2,GetVertex(evtIndex2), gammaj, pdgj)) 
+//        //MEFIX
+//          continue ;
+//        
+//        if(pdgj == AliCaloPID::kPhoton ){
+//          
+//          if((gammai+gammaj).Pt() < GetMinPt() || (gammai+gammaj).Pt() > GetMaxPt()) 
+//            continue ;
+//          
+//            //Selection within angular range
+//          Float_t phi = (gammai+gammaj).Phi();
+//          if(phi < 0) phi+=TMath::TwoPi();
+//            //Float_t deltaphi = TMath::Abs(phiTrig-phi);
+//            //if( (deltaphi < fDeltaPhiMinCut) || ( deltaphi > fDeltaPhiMaxCut) ) continue ;
+//          
+//            //Select good pair (aperture and invariant mass)
+//          if(GetNeutralMesonSelection()->SelectPair(gammai, gammaj)){
+//            
+//            if(GetDebug() > 2 ) printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillAOD() - Neutral Hadron Correlation: AOD Selected gamma pair: pt %2.2f, phi %2.2f, eta %2.2f, M %2.3f\n",
+//                                       (gammai+gammaj).Pt(),(gammai+gammaj).Phi(),(gammai+gammaj).Eta(), (gammai+gammaj).M());
+//            
+//            TLorentzVector pi0mom = gammai+gammaj;
+//            AliAODPWG4Particle pi0 = AliAODPWG4Particle(pi0mom);
+//              //pi0.SetLabel(calo->GetLabel());
+//            pi0.SetPdg(AliCaloPID::kPi0);
+//            pi0.SetDetector(detector);       
+//            if(IsDataMC()){
+//                //Check origin of the candidates
+//              
+//              Int_t label1 = calo->GetLabel();
+//              Int_t label2 = calo2->GetLabel();
+//              Int_t tag1 = GetMCAnalysisUtils()->CheckOrigin(label1, GetReader(), inputi);
+//              Int_t tag2 = GetMCAnalysisUtils()->CheckOrigin(label2, GetReader(), inputj);
+//              
+//              if(GetDebug() > 0) 
+//                printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillAOD() - Origin of: photon1 %d; photon2 %d \n",tag1, tag2);
+//              if(GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCPi0Decay) && GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCPi0Decay)){
+//                
+//                  //Check if pi0 mother is the same
+//                if(GetReader()->ReadStack()){ 
+//                  TParticle * mother1 = GetMCStack()->Particle(label1);//photon in kine tree
+//                  label1 = mother1->GetFirstMother();
+//                    //mother1 = GetMCStack()->Particle(label1);//pi0
+//                  
+//                  TParticle * mother2 = GetMCStack()->Particle(label2);//photon in kine tree
+//                  label2 = mother2->GetFirstMother();
+//                    //mother2 = GetMCStack()->Particle(label2);//pi0
+//                }
+//                else if(GetReader()->ReadAODMCParticles()){
+//                  AliAODMCParticle * mother1 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles(inputi))->At(label1);//photon in kine tree
+//                  label1 = mother1->GetMother();
+//                    //mother1 = GetMCStack()->Particle(label1);//pi0
+//                  AliAODMCParticle * mother2 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles(inputj))->At(label2);//photon in kine tree
+//                  label2 = mother2->GetMother();
+//                    //mother2 = GetMCStack()->Particle(label2);//pi0
+//                }
+//                
+//                  //printf("mother1 %d, mother2 %d\n",label1,label2);
+//                if(label1 == label2)
+//                  GetMCAnalysisUtils()->SetTagBit(tag,AliMCAnalysisUtils::kMCPi0);
+//              }
+//            }//Work with mc information also   
+//            pi0.SetTag(tag);
+//              //Set the indeces of the original caloclusters  
+//            pi0.SetCaloLabel(calo->GetID(), calo2->GetID());
+//            AddAODParticle(pi0);
+//            
+//            
+//          }//Pair selected
+//        }//if pair of gammas
+//      }//2nd loop
+//    }// if pdg = 22
+//  }//1st loop
+//  
+//  if(GetDebug() > 1) 
+//    printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillAOD() - End, %d pi0's found \n",GetOutputAODBranch()->GetEntriesFast());
+//}
+
+//____________________________________________________________________________
+void  AliAnaParticleHadronCorrelation::MakeNeutralCorrelation(AliAODPWG4ParticleCorrelation * const aodParticle, TObjArray* pi0list, const Bool_t bFillHisto)  
 {  
   // Neutral Pion Correlation Analysis
-  if(GetDebug() > 1) printf("Make trigger particle - neutral hadron correlation \n");
+  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 ; 
+  if (GetMixedEvent()) {
+    evtIndex11 = GetMixedEvent()->EventIndexForCaloCluster(aodParticle->GetCaloLabel(0)) ;
+    evtIndex12 = GetMixedEvent()->EventIndexForCaloCluster(aodParticle->GetCaloLabel(1)) ;    
+  }  
   
-  Double_t pt = -100.;
+  Double_t pt   = -100.;
+  Double_t px   = -100.;
+  Double_t py   = -100.;
   Double_t rat = -100.; 
-  Double_t phi = -100. ;
-  Double_t eta = -100. ;
+  Double_t phi = -100.;
+  Double_t eta = -100.;
+  Double_t xE  = -100.; 
+  Double_t cosi  = -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();
   
-  TLorentzVector gammai;
-  TLorentzVector gammaj;
   
-  Double_t vertex[] = {0,0,0};
-
-  if(!GetReader()->GetDataType()== AliCaloTrackReader::kMC) GetReader()->GetVertex(vertex);
+  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. ;
   
-  //Cluster loop, select pairs with good pt, phi and fill AODs or histograms
-  for(Int_t iclus = 0;iclus < pl->GetEntries() ; iclus ++ ){
-    AliAODCaloCluster * calo = (AliAODCaloCluster *) (pl->At(iclus)) ;
-    if(!bFillHisto){
-
-      //Cluster selection, not charged, with photon or pi0 id and in fidutial cut
-      Int_t pdg=0;
-      if(!SelectCluster(calo, vertex, gammai, pdg)) continue ;
+  Double_t ratDecay1  = -100.;  
+  Double_t ratDecay2  = -100.; 
+  Float_t deltaphi = -100. ;
+  Float_t deltaphiDecay1 = -100. ;
+  Float_t deltaphiDecay2 = -100. ;
+  TObjArray * clusters = 0x0 ;  
+  TLorentzVector photonMom ;   
+  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 = GetAODEMCAL() ;
+      else                                    clusters = GetAODPHOS() ;
+      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("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
+  
+  TObjArray * refpi0    =0x0;
+  Int_t nrefs = 0;
+  
+  //Loop on stored AOD pi0
+  Int_t naod = pi0list->GetEntriesFast();
+  if(GetDebug() > 0) printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillHistograms() -  aod branch entries %d\n", naod);
+  for(Int_t iaod = 0; iaod < naod ; iaod++){
+    AliAODPWG4Particle* pi0 =  (AliAODPWG4Particle*) (pi0list->At(iaod));
+    
+    Int_t evtIndex2 = 0 ; 
+    Int_t evtIndex3 = 0 ; 
+    if (GetMixedEvent()) {
+      evtIndex2 = GetMixedEvent()->EventIndexForCaloCluster(pi0->GetCaloLabel(0)) ;
+      evtIndex3 = GetMixedEvent()->EventIndexForCaloCluster(pi0->GetCaloLabel(1)) ;
       
-    if(GetDebug() > 2)
-      printf("neutral cluster: pt %f, phi %f, phi trigger %f. Cuts:  delta phi min %2.2f,  max%2.2f, pT min %2.2f \n",
-            gammai.Pt(),gammai.Phi(),phiTrig,fDeltaPhiMinCut, fDeltaPhiMaxCut, GetMinPt());
+      if (evtIndex11 == evtIndex2 || evtIndex12 == evtIndex2 || evtIndex11 == evtIndex3 || evtIndex12 == evtIndex3) // trigger and pi0 are not from different events
+        continue ; 
+    }      
     
-      //2 gamma overlapped, found with PID
-      if(pdg == AliCaloPID::kPi0){ 
-
-       //Select only hadrons in pt range
-       if(gammai.Pt() < GetMinPt() || gammai.Pt() > GetMaxPt()) continue ;
-
-       aodParticle->AddCluster(calo);
-       if(GetDebug() > 2) printf("Correlated with selected pi0 (pid): pt %f, phi %f",gammai.Pt(),gammai.Phi());
-
-      }// pdg = 111
+    //Int_t pdg = pi0->GetPdg();
+    //if(pdg != AliCaloPID::kPi0) continue;  
+    
+    pt  = pi0->Pt();
+    px  = pi0->Px();
+    py  = pi0->Py();    
+    if(pt < GetMinPt() || pt > GetMaxPt()) continue ;
+    //jumped out this event if near side associated partile pt larger than trigger
+    if(pt > ptTrig && TMath::Abs(phi-phiTrig)<TMath::PiOver2())  break ;
 
-      //Make invariant mass analysis
-      else if(pdg == AliCaloPID::kPhoton){     
-       //Search the photon companion in case it comes from  a Pi0 decay
-       //Apply several cuts to select the good pair;
-       for(Int_t jclus = iclus+1; jclus < pl->GetEntries() ; jclus ++ ){
-         AliAODCaloCluster * calo2 = (AliAODCaloCluster *) (pl->At(jclus)) ;
-         
-         //Cluster selection, not charged with photon or pi0 id and in fidutial cut
-         Int_t pdgj=0;
-         if(!SelectCluster(calo2,vertex, gammaj, pdgj)) continue ;
-
-         if(pdgj == AliCaloPID::kPhoton ){
-           
-           if((gammai+gammaj).Pt() < GetMinPt() || (gammai+gammaj).Pt() > GetMaxPt()) continue ;
-           
-           //Select good pair (aperture and invariant mass)
-           if(GetNeutralMesonSelection()->SelectPair(gammai, gammaj)){
-
-             if(GetDebug() > 2 ) printf("Neutral Hadron Correlation: AOD Selected gamma pair: pt %2.2f, phi %2.2f, eta %2.2f, M %2.3f",
-                                        (gammai+gammaj).Pt(),(gammai+gammaj).Phi(),(gammai+gammaj).Eta(), (gammai+gammaj).M());
-             Int_t labels[]={calo->GetLabel(0),calo2->GetLabel(0)};
-             Float_t pid[]={0,0,0,0,0,0,1,0,0,0,0,0,0};//Pi0 weight 1
-             Float_t pos[]={(gammai+gammaj).X(), (gammai+gammaj).Y(), (gammai+gammaj).Z()};
-
-             AliAODCaloCluster *caloCluster =  new AliAODCaloCluster(0,2,labels,(gammai+gammaj).E(), pos, pid,calo->GetType(),0);
-             aodParticle->AddCluster(caloCluster);
-           }//Pair selected
-         }//if pair of gammas
-       }//2nd loop
-      }// if pdg = 22
-    }// Fill AODs
-    else{ //Fill histograms
-      
-      calo->GetMomentum(gammai,vertex);//Assume that come from vertex in straight line
-      pt  = gammai.Pt();
+    //Selection within angular range
+    phi = pi0->Phi();
+    //Float_t deltaphi = TMath::Abs(phiTrig-phi);
+    //if( (deltaphi < fDeltaPhiMinCut) || ( deltaphi > fDeltaPhiMaxCut) ) continue ;
+    
+    if(bFillHisto){
       
-      if(pt < GetMinPt() || pt > GetMaxPt()) continue ;
+      deltaphi = phiTrig-phi;
+      if(deltaphi<-TMath::PiOver2()) deltaphi+=TMath::TwoPi();
+      if(deltaphi>3*TMath::PiOver2()) deltaphi-=TMath::TwoPi();
       
       rat = pt/ptTrig ;
-      phi = gammai.Phi() ;
-      eta = gammai.Eta() ;
-
-      if(GetDebug() > 2 ) printf("Neutral Hadron Correlation: Histograms selected gamma pair: pt %2.2f, phi %2.2f, eta %2.2f",pt,phi,eta);
+      phi = pi0->Phi() ;
+      eta = pi0->Eta() ;
+      xE   = -(px*pxTrig+py*pyTrig)/(ptTrig*ptTrig);
+      if(xE <0.)xE =-xE;
+      cosi = TMath::Log(1/xE);
       
-      fhEtaNeutral->Fill(ptTrig,eta);
-      fhPhiNeutral->Fill(ptTrig,phi);
+      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("deltaPhoton1 = %f, deltaPhoton2 = %f \n", deltaphiDecay1, deltaphiDecay2);
+          if( (deltaphiDecay1 > fDeltaPhiMinCut) && ( deltaphiDecay1 < fDeltaPhiMaxCut) )
+            fhPtImbalanceDecayNeutral->Fill(ptDecay1,ratDecay1); 
+          if( (deltaphiDecay2 > fDeltaPhiMinCut) && ( deltaphiDecay2 < fDeltaPhiMaxCut) )
+            fhPtImbalanceDecayNeutral->Fill(ptDecay2,ratDecay2);
+          if(GetDebug() > 1)printf("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,phiTrig-phi);
-      fhDeltaPhiNeutralPt->Fill(pt,phiTrig-phi);
-      //Selection within angular range
-      if(((phiTrig-phi)> fDeltaPhiMinCut) && ((phiTrig-phi)<fDeltaPhiMaxCut) ){
-       if(GetDebug() > 2 ) printf("Selected neutral for momentum imbalance: pt %2.2f, phi %2.2f, eta %2.2f ",pt,phi,eta);
-       fhPtImbalanceNeutral->Fill(ptTrig,rat);
-      }    
-    }//Fill histograms
-  }//1st loop
-  
+      fhDeltaPhiNeutral->Fill(ptTrig,deltaphi);
+      fhDeltaPhiDeltaEtaNeutral->Fill(deltaphi,etaTrig-eta);
+      
+      //delta phi cut for correlation
+      if( (deltaphi > fDeltaPhiMinCut) && ( deltaphi < fDeltaPhiMaxCut) ) {
+        fhDeltaPhiNeutralPt->Fill(pt,deltaphi);
+        fhPtImbalanceNeutral->Fill(ptTrig,rat); 
+        fhPtHbpNeutral->Fill(ptTrig,cosi); 
+      }
+      else {
+        fhDeltaPhiUeNeutralPt->Fill(pt,deltaphi);
+        fhPtImbalanceUeNeutral->Fill(ptTrig,rat);
+        fhPtHbpUeNeutral->Fill(ptTrig,cosi); 
+      }
+      //several UE calculation 
+      if(fMakeSeveralUE){
+        if((deltaphi<-fUeDeltaPhiMinCut) && (deltaphi >-fUeDeltaPhiMaxCut)){  
+          fhDeltaPhiUeLeftNeutral->Fill(pt,deltaphi);
+          fhPtImbalanceUeLeftNeutral->Fill(ptTrig,rat);
+          fhPtHbpUeLeftNeutral->Fill(ptTrig,cosi);
+        }
+        if((deltaphi>fUeDeltaPhiMinCut) && (deltaphi <fUeDeltaPhiMaxCut)){  
+          fhDeltaPhiUeRightNeutral->Fill(pt,deltaphi);
+          fhPtImbalanceUeRightNeutral->Fill(ptTrig,rat);
+          fhPtHbpUeRightNeutral->Fill(ptTrig,cosi);
+        }
+      } //several UE calculation
+         }
+    else{
+      nrefs++;
+      if(nrefs==1){
+        refpi0 = new TObjArray(0);
+        refpi0->SetName(GetAODObjArrayName()+"Pi0s");
+        refpi0->SetOwner(kFALSE);
+      }
+      refpi0->Add(pi0);
+    }//put references in trigger AOD 
+      
+     //if(GetDebug() > 2 ) printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelation() - Selected neutral for momentum imbalance: pt %2.2f, phi %2.2f, eta %2.2f \n",pt,phi,eta);
+      
+    }//loop
 }
+  
 
 //____________________________________________________________________________
-Bool_t  AliAnaParticleHadronCorrelation::SelectCluster(AliAODCaloCluster * calo, Double_t *vertex, TLorentzVector & mom, Int_t & pdg) const {
-   //Select cluster depending on its pid and acceptance selections
-   
-   //Skip matched clusters with tracks
-  if(calo->GetNTracksMatched() > 0) {
-  return kFALSE;} 
-   
-   //Check PID
-   calo->GetMomentum(mom,vertex);//Assume that come from vertex in straight line
-   pdg = AliCaloPID::kPhoton;   
-   if(IsCaloPIDOn()){
-     //Get most probable PID, 2 options check PID weights (in MC this option is mandatory)
-     //or redo PID, recommended option for EMCal.
-     TString detector = "";
-     if(calo->IsPHOSCluster()) detector= "PHOS";
-     else if(calo->IsEMCALCluster()) detector= "EMCAL";
-
-     if(!IsCaloPIDRecalculationOn() || GetReader()->GetDataType() == AliCaloTrackReader::kMC )
-       pdg = GetCaloPID()->GetPdg(detector,calo->PID(),mom.E());//PID with weights
-     else
-       pdg = GetCaloPID()->GetPdg(detector,mom,calo);//PID recalculated
-     
-     if(GetDebug() > 1) printf("PDG of identified particle %d\n",pdg);
-     
-     //If it does not pass pid, skip
-     if(pdg != AliCaloPID::kPhoton || pdg != AliCaloPID::kPi0) {
-       return kFALSE ;
-     }
-   }
-   
-   //Check acceptance selection
-   if(IsFidutialCutOn()){
-     Bool_t in = kFALSE;
-     if(calo->IsPHOSCluster())
-       in = GetFidutialCut()->IsInFidutialCut(mom,"PHOS") ;
-     else if(calo->IsEMCALCluster())
-       in = GetFidutialCut()->IsInFidutialCut(mom,"EMCAL") ;
-     if(! in ) { return kFALSE ;}
-   }
-   
-   if(GetDebug() > 1) printf("Correlation photon selection cuts passed: pT %3.2f, pdg %d\n",mom.Pt(), pdg);
-   
-   return kTRUE;
-   
- }
+//Bool_t  AliAnaParticleHadronCorrelation::SelectCluster(AliVCluster * calo, Double_t *vertex, TLorentzVector & mom, Int_t & pdg) {
+//  //Select cluster depending on its pid and acceptance selections
+//  
+//  //Skip matched clusters with tracks
+//  if(IsTrackMatched(calo)) return kFALSE;
+//  
+//  TString detector = "";
+//  if     (calo->IsPHOS())  detector= "PHOS";
+//  else if(calo->IsEMCAL()) detector= "EMCAL";
+//             
+//  //Check PID
+//  calo->GetMomentum(mom,vertex);//Assume that come from vertex in straight line
+//  pdg = AliCaloPID::kPhoton;   
+//  if(IsCaloPIDOn()){
+//    //Get most probable PID, 2 options check PID weights (in MC this option is mandatory)
+//    //or redo PID, recommended option for EMCal.
+//    
+//    if(!IsCaloPIDRecalculationOn() || GetReader()->GetDataType() == AliCaloTrackReader::kMC )
+//      pdg = GetCaloPID()->GetPdg(detector,calo->GetPID(),mom.E());//PID with weights
+//    else
+//      pdg = GetCaloPID()->GetPdg(detector,mom,calo);//PID recalculated
+//    
+//    if(GetDebug() > 5) printf("AliAnaParticleHadronCorrelation::SelectCluster() - PDG of identified particle %d\n",pdg);
+//    
+//    //If it does not pass pid, skip
+//    if(pdg != AliCaloPID::kPhoton && pdg != AliCaloPID::kPi0) {
+//      return kFALSE ;
+//    }
+//  }//PID on
+//  
+//  //Check acceptance selection
+//  if(IsFiducialCutOn()){
+//    Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,detector) ;
+//    if(! in ) return kFALSE ;
+//  }
+//  
+//  if(GetDebug() > 5) printf("AliAnaParticleHadronCorrelation::SelectCluster() - Correlation photon selection cuts passed: pT %3.2f, pdg %d\n",mom.Pt(), pdg);
+//  
+//  return kTRUE;
+//  
+//}