]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGGA/CaloTrackCorrelations/AliAnaParticleHadronCorrelation.cxx
change name of variable, conflicts with enum with detector types
[u/mrichter/AliRoot.git] / PWGGA / CaloTrackCorrelations / AliAnaParticleHadronCorrelation.cxx
index 73200b1dc9e4288dc19d8899df620ecf92509312..6afc9bca363d010a920c7f3c11dbbed3414f0d3f 100755 (executable)
@@ -61,6 +61,7 @@ ClassImp(AliAnaParticleHadronCorrelation)
     fUeDeltaPhiMaxCut(0.),          fUeDeltaPhiMinCut(0.), 
     fPi0AODBranchName(""),          fNeutralCorr(0),       
     fPi0Trigger(0),                 fDecayTrigger(0),
+    fNDecayBits(0),                 fDecayBits(),
     fMakeAbsoluteLeading(0),        fMakeNearSideLeading(0),      
     fLeadingTriggerIndex(-1),       fHMPIDCorrelation(0),  fFillBradHisto(0),
     fNAssocPtBins(0),               fAssocPtBinLimit(),
@@ -68,12 +69,13 @@ ClassImp(AliAnaParticleHadronCorrelation)
     fListMixTrackEvents(),          fListMixCaloEvents(),
     fUseMixStoredInReader(0),       fFillNeutralEventMixPool(0),
     fM02MaxCut(0),                  fM02MinCut(0),  
-    fFillPileUpHistograms(0),       fFillHighMultHistograms(0),
-    fSelectLeadingHadronAngle(0),
+    fSelectLeadingHadronAngle(0),   fFillLeadHadOppositeHisto(0),
     fMinLeadHadPhi(0),              fMaxLeadHadPhi(0),
     fMinLeadHadPt(0),               fMaxLeadHadPt(0),
     fFillEtaGapsHisto(1),           fFillMomImbalancePtAssocBinsHisto(0),
     fMCGenTypeMin(0),               fMCGenTypeMax(0),
+    fTrackVector(),                 fMomentum(),
+    fDecayMom1(),                   fDecayMom2(),
     //Histograms
     fhPtTriggerInput(0),            fhPtTriggerSSCut(0),
     fhPtTriggerIsoCut(0),           fhPtTriggerFidCut(0),
@@ -81,11 +83,13 @@ ClassImp(AliAnaParticleHadronCorrelation)
     fhPtTriggerVzBin(0),            fhPtTriggerBin(0),                 
     fhPhiTrigger(0),                fhEtaTrigger(0),   
     fhPtTriggerMC(),
+    fhPtDecayTrigger(),             fhPtDecayTriggerMC(),
     fhPtTriggerCentrality(0),       fhPtTriggerEventPlane(0), 
     fhTriggerEventPlaneCentrality(0),
     fhPtTriggerMixed(0),            fhPtTriggerMixedVzBin(0), fhPtTriggerMixedBin(0),              
     fhPhiTriggerMixed(0),           fhEtaTriggerMixed(0),
     fhPtLeadingOppositeHadron(0),   fhPtDiffPhiLeadingOppositeHadron(0), fhPtDiffEtaLeadingOppositeHadron(0),
+    fhPtNoLeadingOppositeHadron(0), fhEtaPhiNoLeadingOppositeHadron(0),
     fhDeltaPhiDeltaEtaCharged(0),
     fhPhiCharged(0),                fhEtaCharged(0), 
     fhDeltaPhiCharged(0),           fhDeltaEtaCharged(0), 
@@ -148,9 +152,10 @@ ClassImp(AliAnaParticleHadronCorrelation)
     fhDeltaPhiUeLeftNeutral(0),     fhXEUeLeftNeutral(0),
     fhPtHbpXEUeLeftNeutral(0),      fhZTUeLeftNeutral(0),
     fhPtHbpZTUeLeftNeutral(0),      fhPtPi0DecayRatio(0),
-    fhDeltaPhiDecayCharged(0),      fhXEDecayCharged(0),           fhZTDecayCharged(0),
-    fhDeltaPhiDecayNeutral(0),      fhXEDecayNeutral(0),           fhZTDecayNeutral(0),
-    fhDeltaPhiDecayChargedAssocPtBin(0), 
+    fhDeltaPhiPi0DecayCharged(0),   fhXEPi0DecayCharged(0),        fhZTPi0DecayCharged(0),
+    fhDeltaPhiPi0DecayNeutral(0),   fhXEPi0DecayNeutral(0),        fhZTPi0DecayNeutral(0),
+    fhDeltaPhiDecayCharged(),       fhXEDecayCharged(),            fhZTDecayCharged(),
+    fhDeltaPhiDecayChargedAssocPtBin(),
     fhMCPtTrigger(),                fhMCPhiTrigger(),              fhMCEtaTrigger(),
     fhMCPtTriggerNotLeading(),      fhMCPhiTriggerNotLeading(),    fhMCEtaTriggerNotLeading(),
     fhMCEtaCharged(),               fhMCPhiCharged(),
@@ -181,13 +186,16 @@ ClassImp(AliAnaParticleHadronCorrelation)
   //Initialize parameters
   InitParameters();
   
-  for(Int_t i = 0; i < 7; i++)
+  for(Int_t i = 0; i < fgkNmcTypes; i++)
   { 
     fhPtTriggerMC[i] = 0;
     fhXEChargedMC[i] = 0;
     fhDeltaPhiChargedMC[i] = 0;
+    for(Int_t ib = 0; ib < 4; ib++)  fhPtDecayTriggerMC[ib][i] = 0;
   }
-  
+
+  for(Int_t ib = 0; ib < 4; ib++)  fhPtDecayTrigger[ib] = 0;
+
   for(Int_t i = 0; i < 7; i++)
   {
     fhPtTriggerPileUp             [i] = 0 ;
@@ -250,7 +258,7 @@ AliAnaParticleHadronCorrelation::~AliAnaParticleHadronCorrelation()
 void AliAnaParticleHadronCorrelation::FillChargedAngularCorrelationHistograms(Float_t ptAssoc,  Float_t ptTrig,      Int_t   bin,
                                                                               Float_t phiAssoc, Float_t phiTrig,     Float_t deltaPhi,
                                                                               Float_t etaAssoc, Float_t etaTrig,  
-                                                                              Bool_t  decay,    Float_t hmpidSignal, Int_t  outTOF,
+                                                                              Int_t   decayTag, Float_t hmpidSignal, Int_t  outTOF,
                                                                               Int_t   cen,      Int_t   mcTag)
 {
   // Fill angular correlation related histograms
@@ -274,7 +282,7 @@ void AliAnaParticleHadronCorrelation::FillChargedAngularCorrelationHistograms(Fl
   
   // Pile up studies
   
-  if(fFillPileUpHistograms)
+  if(IsPileUpAnalysisOn())
   {
     if     (outTOF==1)
     {
@@ -318,9 +326,17 @@ void AliAnaParticleHadronCorrelation::FillChargedAngularCorrelationHistograms(Fl
   {
     Int_t mcIndex = GetMCTagHistogramIndex(mcTag);
     fhDeltaPhiChargedMC[mcIndex]->Fill(ptTrig , deltaPhi);
-  }  
+    if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCDecayPairLost) && mcIndex==2 )
+      fhDeltaPhiChargedMC[7]->Fill(ptTrig , deltaPhi);
+  }
   
-  if(fDecayTrigger && decay) fhDeltaPhiDecayCharged   ->Fill(ptTrig  , deltaPhi);      
+  if(fDecayTrigger && decayTag > 0)
+  {
+    for(Int_t ibit = 0; ibit<fNDecayBits; ibit++)
+    {
+      if(GetNeutralMesonSelection()->CheckDecayBit(decayTag,fDecayBits[ibit])) fhDeltaPhiDecayCharged[ibit]->Fill(ptTrig,deltaPhi);
+    }
+  }
   
   Double_t  dphiBrad = -100;
   if(fFillBradHisto)
@@ -354,7 +370,8 @@ void AliAnaParticleHadronCorrelation::FillChargedAngularCorrelationHistograms(Fl
     if (fFillBradHisto)
       fhDeltaPhiBradAssocPtBin        [bin]->Fill(ptTrig, dphiBrad);
     
-    if(fDecayTrigger && decay)
+    if(fDecayTrigger && decayTag > 0 && fNDecayBits > 0 &&
+       GetNeutralMesonSelection()->CheckDecayBit(decayTag,fDecayBits[0]))
       fhDeltaPhiDecayChargedAssocPtBin[bin]->Fill(ptTrig, deltaPhi);      
     
     if(fHMPIDCorrelation)
@@ -374,7 +391,7 @@ void AliAnaParticleHadronCorrelation::FillChargedAngularCorrelationHistograms(Fl
   }
   
   //fill different multiplicity/centrality histogram
-  if(fFillHighMultHistograms && cen >= 0 && cen < GetNCentrBin())
+  if(IsHighMultiplicityAnalysisOn() && cen >= 0 && cen < GetNCentrBin())
   {
     fhDeltaPhiChargedMult[cen]->Fill(ptTrig,deltaPhi);
     fhDeltaEtaChargedMult[cen]->Fill(ptTrig,deltaEta);
@@ -384,7 +401,7 @@ void AliAnaParticleHadronCorrelation::FillChargedAngularCorrelationHistograms(Fl
 //___________________________________________________________________________________________________________________________________
 Bool_t AliAnaParticleHadronCorrelation::FillChargedMCCorrelationHistograms(Float_t mcAssocPt, Float_t mcAssocPhi, Float_t mcAssocEta,
                                                                            Float_t mcTrigPt,  Float_t mcTrigPhi,  Float_t mcTrigEta,
-                                                                           Int_t histoIndex)
+                                                                           Int_t histoIndex, Bool_t lostDecayPair)
 {
   // Fill MC histograms independently of AOD or ESD
   
@@ -397,26 +414,29 @@ Bool_t AliAnaParticleHadronCorrelation::FillChargedMCCorrelationHistograms(Float
   if( fMakeAbsoluteLeading && mcAssocPt > mcTrigPt )     lead = kFALSE; // skip event
   
   // Skip this event if near side associated particle pt larger than trigger
-  if( fMakeNearSideLeading && mcAssocPt > mcTrigPt &&
-     TMath::Abs(mcAssocPhi-mcTrigPhi)<TMath::PiOver2() ) lead = kFALSE; // skip event
+  if( mcAssocPhi < 0 ) mcAssocPhi+=TMath::TwoPi();
+  
+  Float_t mcdeltaPhi= mcTrigPhi-mcAssocPhi;
+  if(mcdeltaPhi <= -TMath::PiOver2()) mcdeltaPhi+=TMath::TwoPi();
+  if(mcdeltaPhi > 3*TMath::PiOver2()) mcdeltaPhi-=TMath::TwoPi();
+
+  if( fMakeNearSideLeading)
+  {
+    if( mcAssocPt > mcTrigPt && mcdeltaPhi < TMath::PiOver2() ) lead = kFALSE; // skip event
+  }
   
   //
   // Select only hadrons in pt range
   if ( mcAssocPt < fMinAssocPt || mcAssocPt > fMaxAssocPt ) return lead ; // exclude but continue
   if ( mcAssocPt < GetReader()->GetCTSPtMin())              return lead ;
 
-  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 && 
+     mcdeltaPhi                       < 1e-6 &&
      TMath::Abs(mcAssocEta-mcTrigEta) < 1e-6)            return lead ; // exclude but continue
   
-  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);
@@ -453,6 +473,28 @@ Bool_t AliAnaParticleHadronCorrelation::FillChargedMCCorrelationHistograms(Float
     fhMCPtTrigPout       [histoIndex]->Fill(mcTrigPt, mcpout) ;
   }
 
+  if(histoIndex==2 && lostDecayPair && 7 >= fMCGenTypeMin && 7 <= fMCGenTypeMax )
+  {
+    fhMCEtaCharged     [7]->Fill(mcAssocPt, mcAssocEta);
+    fhMCPhiCharged     [7]->Fill(mcAssocPt, mcAssocPhi);
+    fhMCDeltaEtaCharged[7]->Fill(mcTrigPt , mcTrigEta-mcAssocEta);
+    fhMCDeltaPhiCharged[7]->Fill(mcTrigPt , mcdeltaPhi);
+    fhMCPtAssocDeltaPhi[7]->Fill(mcAssocPt, mcdeltaPhi);
+    
+    fhMCDeltaPhiDeltaEtaCharged[7]->Fill(mcdeltaPhi,mcTrigEta-mcAssocEta);
+    
+    //delta phi cut for correlation
+    if( (mcdeltaPhi > fDeltaPhiMinCut) && (mcdeltaPhi < fDeltaPhiMaxCut) )
+    {
+      fhMCDeltaPhiChargedPt[7]->Fill(mcAssocPt,mcdeltaPhi);
+      fhMCPtXECharged      [7]->Fill(mcTrigPt, mcxE);
+      fhMCPtHbpXECharged   [7]->Fill(mcTrigPt, mchbpXE);
+      fhMCPtZTCharged      [7]->Fill(mcTrigPt, mczT);
+      fhMCPtHbpZTCharged   [7]->Fill(mcTrigPt, mchbpZT);
+      fhMCPtTrigPout       [7]->Fill(mcTrigPt, mcpout) ;
+    }
+  }
+  
   // Underlying event
   
   // Right
@@ -474,6 +516,17 @@ Bool_t AliAnaParticleHadronCorrelation::FillChargedMCCorrelationHistograms(Float
     if(mcUezT > 0) fhMCPtHbpZTUeCharged[histoIndex]->Fill(mcTrigPt,TMath::Log(1/mcUezT));
     
     fhMCUePart[histoIndex]->Fill(mcTrigPt);
+    
+    if(histoIndex==2 && lostDecayPair && 7 >= fMCGenTypeMin && 7 <= fMCGenTypeMax )
+    {
+      fhMCPtXEUeCharged[7]->Fill(mcTrigPt,mcUexE);
+      if(mcUexE > 0) fhMCPtHbpXEUeCharged[histoIndex]->Fill(mcTrigPt,TMath::Log(1/mcUexE));
+      
+      fhMCPtZTUeCharged[7]->Fill(mcTrigPt,mcUezT);
+      if(mcUezT > 0) fhMCPtHbpZTUeCharged[histoIndex]->Fill(mcTrigPt,TMath::Log(1/mcUezT));
+      
+      fhMCUePart[7]->Fill(mcTrigPt);
+    }
   }
 
   if(fMakeSeveralUE)
@@ -494,6 +547,15 @@ Bool_t AliAnaParticleHadronCorrelation::FillChargedMCCorrelationHistograms(Float
       
       fhMCPtZTUeLeftCharged[histoIndex]->Fill(mcTrigPt,mcUezT);
       if(mcUezT > 0) fhMCPtHbpZTUeLeftCharged[histoIndex]->Fill(mcTrigPt,TMath::Log(1/mcUezT));
+      
+      if(histoIndex==2 && lostDecayPair && 7 >= fMCGenTypeMin && 7 <= fMCGenTypeMax )
+      {
+        fhMCPtXEUeLeftCharged[7]->Fill(mcTrigPt,mcUexE);
+        if(mcUexE > 0) fhMCPtHbpXEUeLeftCharged[7]->Fill(mcTrigPt,TMath::Log(1/mcUexE));
+        
+        fhMCPtZTUeLeftCharged[7]->Fill(mcTrigPt,mcUezT);
+        if(mcUezT > 0) fhMCPtHbpZTUeLeftCharged[7]->Fill(mcTrigPt,TMath::Log(1/mcUezT));
+      }
     }
   }
   
@@ -504,7 +566,7 @@ Bool_t AliAnaParticleHadronCorrelation::FillChargedMCCorrelationHistograms(Float
 void AliAnaParticleHadronCorrelation::FillChargedMomentumImbalanceHistograms(Float_t ptTrig,   Float_t ptAssoc, 
                                                                              Float_t deltaPhi,
                                                                              Int_t   cen,      Int_t   charge,
-                                                                             Int_t   bin,      Bool_t  decay,
+                                                                             Int_t   bin,      Int_t   decayTag,
                                                                              Int_t   outTOF,   Int_t   mcTag)
 
 {
@@ -541,10 +603,12 @@ void AliAnaParticleHadronCorrelation::FillChargedMomentumImbalanceHistograms(Flo
   {
     Int_t mcIndex = GetMCTagHistogramIndex(mcTag);
     fhXEChargedMC[mcIndex]->Fill(ptTrig , xE);
+    if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCDecayPairLost) && mcIndex==2 )
+      fhXEChargedMC[7]->Fill(ptTrig , xE);
   }
   
   // Pile up studies
-  if(fFillPileUpHistograms) 
+  if(IsPileUpAnalysisOn())
   {
     if     (outTOF==1)
     {
@@ -576,11 +640,17 @@ void AliAnaParticleHadronCorrelation::FillChargedMomentumImbalanceHistograms(Flo
     if(GetReader()->IsPileUpFromNotSPDAndNotEMCal())  { fhXEChargedPileUp[6]->Fill(ptTrig,xE); fhZTChargedPileUp[6]->Fill(ptTrig,zT); fhPtTrigChargedPileUp[6]->Fill(ptTrig,ptAssoc); }
   }
   
-  if(fDecayTrigger && decay)
-  {          
-    fhXEDecayCharged->Fill(ptTrig,xE); 
-    fhZTDecayCharged->Fill(ptTrig,zT);
-  } // photon decay pi0/eta trigger        
+  if(fDecayTrigger && decayTag > 0)
+  {
+    for(Int_t ibit = 0; ibit<fNDecayBits; ibit++)
+    {
+      if(GetNeutralMesonSelection()->CheckDecayBit(decayTag,fDecayBits[ibit]))
+      {
+        fhXEDecayCharged[ibit]->Fill(ptTrig,xE);
+        fhZTDecayCharged[ibit]->Fill(ptTrig,zT);
+      }
+    }
+  } // photon decay pi0/eta trigger
   
   if(bin >= 0 && fFillMomImbalancePtAssocBinsHisto)//away side
   {
@@ -607,7 +677,7 @@ void AliAnaParticleHadronCorrelation::FillChargedMomentumImbalanceHistograms(Flo
   }
   
   //fill different multiplicity/centrality histogram
-  if(fFillHighMultHistograms && cen >= 0 && cen < GetNCentrBin())
+  if(IsHighMultiplicityAnalysisOn() && cen >= 0 && cen < GetNCentrBin())
   {
     fhXEMult[cen]->Fill(ptTrig,xE);
     fhZTMult[cen]->Fill(ptTrig,zT);
@@ -640,7 +710,7 @@ void AliAnaParticleHadronCorrelation::FillChargedUnderlyingEventHistograms(Float
   
   // Pile up studies
   
-  if(fFillPileUpHistograms)
+  if(IsPileUpAnalysisOn())
   {
     if     (outTOF==1)
     {
@@ -670,7 +740,7 @@ void AliAnaParticleHadronCorrelation::FillChargedUnderlyingEventHistograms(Float
   }
   
   //fill different multiplicity/centrality histogram
-  if(fFillHighMultHistograms && cen >= 0 && cen < GetNCentrBin())
+  if(IsHighMultiplicityAnalysisOn() && cen >= 0 && cen < GetNCentrBin())
   {
     fhXEUeMult[cen]->Fill(ptTrig,uexE);
     fhZTUeMult[cen]->Fill(ptTrig,uezT);
@@ -758,26 +828,24 @@ void AliAnaParticleHadronCorrelation::FillChargedUnderlyingEventSidesHistograms(
   }
 } 
 
-//______________________________________________________________________________________________________________________________
-void AliAnaParticleHadronCorrelation::FillDecayPhotonCorrelationHistograms(Float_t ptAssoc,     Float_t phiAssoc, 
-                                                                           TLorentzVector mom1, TLorentzVector mom2,
-                                                                           Bool_t bChargedOrNeutral)
+//_____________________________________________________________________________________________________________________________________
+void AliAnaParticleHadronCorrelation::FillDecayPhotonCorrelationHistograms(Float_t ptAssoc, Float_t phiAssoc, 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 ptDecay1 = fDecayMom1.Pt();
+  Float_t ptDecay2 = fDecayMom2.Pt();
   
   Float_t zTDecay1 = -100, zTDecay2 = -100;
   if(ptDecay1 > 0) zTDecay1 = ptAssoc/ptDecay1 ;
   if(ptDecay2 > 0) zTDecay2 = ptAssoc/ptDecay2 ;
   
-  Float_t deltaPhiDecay1 = mom1.Phi()-phiAssoc;
+  Float_t deltaPhiDecay1 = fDecayMom1.Phi()-phiAssoc;
   if(deltaPhiDecay1< -TMath::PiOver2()) deltaPhiDecay1+=TMath::TwoPi();
   if(deltaPhiDecay1>3*TMath::PiOver2()) deltaPhiDecay1-=TMath::TwoPi();
   
-  Float_t deltaPhiDecay2 = mom2.Phi()-phiAssoc;
+  Float_t deltaPhiDecay2 = fDecayMom2.Phi()-phiAssoc;
   if(deltaPhiDecay2< -TMath::PiOver2()) deltaPhiDecay2+=TMath::TwoPi();
   if(deltaPhiDecay2>3*TMath::PiOver2()) deltaPhiDecay2-=TMath::TwoPi();
   
@@ -786,38 +854,38 @@ void AliAnaParticleHadronCorrelation::FillDecayPhotonCorrelationHistograms(Float
   
   if(bChargedOrNeutral) // correlate with charges
   {
-    fhDeltaPhiDecayCharged->Fill(ptDecay1, deltaPhiDecay1);
-    fhDeltaPhiDecayCharged->Fill(ptDecay2, deltaPhiDecay2);
+    fhDeltaPhiPi0DecayCharged->Fill(ptDecay1, deltaPhiDecay1);
+    fhDeltaPhiPi0DecayCharged->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); 
+      fhZTPi0DecayCharged->Fill(ptDecay1,zTDecay1); 
+      fhXEPi0DecayCharged->Fill(ptDecay1,xEDecay1); 
     }
     if( (deltaPhiDecay2 > fDeltaPhiMinCut) && ( deltaPhiDecay2 < fDeltaPhiMaxCut) )
     {
-      fhZTDecayCharged->Fill(ptDecay2,zTDecay2);
-      fhXEDecayCharged->Fill(ptDecay2,xEDecay2);
+      fhZTPi0DecayCharged->Fill(ptDecay2,zTDecay2);
+      fhXEPi0DecayCharged->Fill(ptDecay2,xEDecay2);
     }
   }
   else // correlate with neutrals
   {
-    fhDeltaPhiDecayNeutral->Fill(ptDecay1, deltaPhiDecay1);
-    fhDeltaPhiDecayNeutral->Fill(ptDecay2, deltaPhiDecay2);
+    fhDeltaPhiPi0DecayNeutral->Fill(ptDecay1, deltaPhiDecay1);
+    fhDeltaPhiPi0DecayNeutral->Fill(ptDecay2, deltaPhiDecay2);
     
     if(GetDebug() > 1) printf("AliAnaParticleHadronCorrelation::FillDecayPhotonHistograms(Neutral corr) - deltaPhoton1 = %f, deltaPhoton2 = %f \n", deltaPhiDecay1, deltaPhiDecay2);
     
     if( (deltaPhiDecay1 > fDeltaPhiMinCut) && ( deltaPhiDecay1 < fDeltaPhiMaxCut) )
     {
-      fhZTDecayNeutral->Fill(ptDecay1,zTDecay1);
-      fhXEDecayNeutral->Fill(ptDecay1,xEDecay1);
+      fhZTPi0DecayNeutral->Fill(ptDecay1,zTDecay1);
+      fhXEPi0DecayNeutral->Fill(ptDecay1,xEDecay1);
     }
     if( (deltaPhiDecay2 > fDeltaPhiMinCut) && ( deltaPhiDecay2 < fDeltaPhiMaxCut) )
     {
-      fhZTDecayNeutral->Fill(ptDecay2,zTDecay2);
-      fhXEDecayNeutral->Fill(ptDecay2,xEDecay2);
+      fhZTPi0DecayNeutral->Fill(ptDecay2,zTDecay2);
+      fhXEPi0DecayNeutral->Fill(ptDecay2,xEDecay2);
     }    
   }
 }
@@ -903,20 +971,18 @@ void AliAnaParticleHadronCorrelation::FillChargedEventMixPool()
   
   TList * pool = fListMixTrackEvents[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();
+    fTrackVector.SetXYZ(track->Px(),track->Py(),track->Pz());
+    Float_t pt   = fTrackVector.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");
+    AliAODPWG4Particle * mixedTrack = new AliAODPWG4Particle(track->Px(),track->Py(),track->Pz(),0);
+    mixedTrack->SetDetectorTag(kCTS);
     mixedTrack->SetChargedBit(track->Charge()>0);
     mixEventTracks->Add(mixedTrack);
   }
@@ -982,8 +1048,6 @@ void AliAnaParticleHadronCorrelation::FillNeutralEventMixPool()
   
   TList * poolCalo = fListMixCaloEvents[eventBin];
   
-  TLorentzVector mom;
-
   for(Int_t ipr = 0;ipr <  pl->GetEntriesFast() ; ipr ++ )
   {
     AliVCluster * calo = (AliVCluster *) (pl->At(ipr)) ;
@@ -994,20 +1058,20 @@ void AliAnaParticleHadronCorrelation::FillNeutralEventMixPool()
     //Cluster momentum calculation
     if(GetReader()->GetDataType() != AliCaloTrackReader::kMC)
     {
-      calo->GetMomentum(mom,GetVertex(0)) ;
+      calo->GetMomentum(fMomentum,GetVertex(0)) ;
     }//Assume that come from vertex in straight line
     else
     {
       Double_t vertex[]={0,0,0};
-      calo->GetMomentum(mom,vertex) ;
+      calo->GetMomentum(fMomentum,vertex) ;
     }
     
-    Float_t pt = mom.Pt();
+    Float_t pt = fMomentum.Pt();
     //Select only clusters in pt range
     if(pt < fMinAssocPt || pt > fMaxAssocPt) continue ;
     
-    AliAODPWG4Particle * mixedCalo = new AliAODPWG4Particle(mom);
-    mixedCalo->SetDetector("EMCAL");
+    AliAODPWG4Particle * mixedCalo = new AliAODPWG4Particle(fMomentum);
+    mixedCalo->SetDetectorTag(kEMCAL);
     mixEventCalo->Add(mixedCalo);
   }
   
@@ -1036,57 +1100,84 @@ Bool_t AliAnaParticleHadronCorrelation::FindLeadingOppositeHadronInWindow(AliAOD
   // Select events where the leading charged particle in the opposite hemisphere
   // to the trigger particle is in a window centered at 180 from the trigger
   
-  Float_t phiTrig    = particle->Phi();
   Float_t etaTrig    = particle->Eta();
   Float_t ptTrig     = particle->Pt();
-  Float_t ptLeadHad  = -100 ;
+  Float_t phiTrig    = particle->Phi();
+  if(phiTrig < 0 ) phiTrig+= TMath::TwoPi();
+
+  Float_t ptLeadHad  = 0 ;
+  Float_t dphiLeadHad= -100 ;
   Float_t phiLeadHad = -100 ;
   Float_t etaLeadHad = -100 ;
-  TVector3 p3;
-  
+  Int_t   nTrack     = 0;
+
   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]);
+    fTrackVector.SetXYZ(track->Px(),track->Py(),track->Pz());
     
-    Float_t pt   = p3.Pt();
-    Float_t phi  = p3.Phi() ;
+    Float_t pt   = fTrackVector.Pt();
+    Float_t phi  = fTrackVector.Phi() ;
     if(phi < 0 ) phi+= TMath::TwoPi();
     
-    if(pt > ptLeadHad && TMath::Abs(phi-phiTrig) > TMath::PiOver2()) // in opposite hemisphere
+    Float_t deltaPhi = phiTrig-phi;
+    //
+    // Calculate deltaPhi shift so that for the particles on the opposite side
+    // it is defined between 90 and 270 degrees
+    // Shift [-360,-90]  to [0, 270]
+    // and [270,360] to [-90,0]
+    if(deltaPhi <= -TMath::PiOver2()) deltaPhi+=TMath::TwoPi();
+    if(deltaPhi > 3*TMath::PiOver2()) deltaPhi-=TMath::TwoPi();
+
+    if(pt > ptLeadHad && deltaPhi > TMath::PiOver2()) // in opposite hemisphere
     {
       ptLeadHad  = pt ;
       phiLeadHad = phi;
-      etaLeadHad = p3.Eta();
+      dphiLeadHad= deltaPhi;
+      etaLeadHad = fTrackVector.Eta();
+      nTrack++;
     }
   }// track loop
   
-  fhPtLeadingOppositeHadron       ->Fill(ptTrig, ptLeadHad);
-  fhPtDiffPhiLeadingOppositeHadron->Fill(ptTrig,phiLeadHad-phiTrig);
-  fhPtDiffEtaLeadingOppositeHadron->Fill(ptTrig,etaLeadHad-etaTrig);
+  if(fFillLeadHadOppositeHisto)
+  {
+    if(nTrack == 0)
+    {
+      fhPtNoLeadingOppositeHadron    ->Fill(ptTrig);
+      fhEtaPhiNoLeadingOppositeHadron->Fill(etaTrig,phiTrig);
+    }
+    else
+    {
+      fhPtLeadingOppositeHadron       ->Fill(ptTrig,  ptLeadHad);
+      fhPtDiffPhiLeadingOppositeHadron->Fill(ptTrig,dphiLeadHad);
+      fhPtDiffEtaLeadingOppositeHadron->Fill(ptTrig, etaLeadHad-etaTrig);
+    }
+  }
   
   if(GetDebug() > 1 )
   {
-    printf("AliAnaParticleHadronCorrelation::FindLeadingOppositeHadronInWindow() pT %2.2f, phi %2.2f, eta %2.2f\n",
-           ptLeadHad,phiLeadHad*TMath::RadToDeg(),etaLeadHad);
+    printf("AliAnaParticleHadronCorrelation::FindLeadingOppositeHadronInWindow() pT %2.2f, phi %2.2f, eta %2.2f, nTracks away %d, total tracks %d\n",
+           ptLeadHad,phiLeadHad*TMath::RadToDeg(),etaLeadHad,nTrack, GetTrackMultiplicity());
     
     printf("\t  pT trig %2.2f, Dphi (trigger-hadron) %2.2f, Deta (trigger-hadron) %2.2f\n",
-           ptTrig, (phiLeadHad-phiTrig)*TMath::RadToDeg(), etaLeadHad-etaTrig);
+           ptTrig, dphiLeadHad*TMath::RadToDeg(), etaLeadHad-etaTrig);
     printf("\t cuts pT: min %2.2f, max %2.2f; DPhi: min %2.2f, max %2.2f\n",fMinLeadHadPt,fMaxLeadHadPt,fMinLeadHadPhi*TMath::RadToDeg(),fMaxLeadHadPhi*TMath::RadToDeg());
   }
   
-  if( ptLeadHad < fMinLeadHadPt ||
-      ptLeadHad > fMaxLeadHadPt ) return kFALSE;
+  // reject the trigger if the leading hadron is not in the requested pt or phi window and
+  
+  if( nTrack == 0 ) return kFALSE; // No track found in opposite hemisphere
+  
+  if( ptLeadHad < fMinLeadHadPt || ptLeadHad > fMaxLeadHadPt ) return kFALSE;
   
   //printf("Accept leading hadron pT \n");
   
-  if( TMath::Abs(phiLeadHad-phiTrig) < fMinLeadHadPhi ||
-      TMath::Abs(phiLeadHad-phiTrig) > fMaxLeadHadPhi ) return kFALSE;
+  if( dphiLeadHad < fMinLeadHadPhi || dphiLeadHad > fMaxLeadHadPhi ) return kFALSE;
   
   //printf("Accept leading hadron phi \n");
   
+  
   return kTRUE ;
 }
 
@@ -1161,7 +1252,7 @@ TList *  AliAnaParticleHadronCorrelation::GetCreateOutputObjects()
   
   Int_t nMixBins = GetNCentrBin()*GetNZvertBin()*GetNRPBin();
   
-  TString nameMC[]     = {"Photon","Pi0","Pi0Decay","EtaDecay","OtherDecay","Electron","Hadron"};
+  TString nameMC[]     = {"Photon","Pi0","Pi0Decay","EtaDecay","OtherDecay","Electron","Hadron","Pi0DecayLostPair"};
   TString pileUpName[] = {"SPD","EMCAL","SPDOrEMCAL","SPDAndEMCAL","SPDAndNotEMCAL","EMCALAndNotSPD","NotSPDAndNotEMCAL"} ;
   
   // For vz dependent histograms, if option ON
@@ -1202,7 +1293,7 @@ TList *  AliAnaParticleHadronCorrelation::GetCreateOutputObjects()
   
   if(IsDataMC())
   {
-    for(Int_t i=0; i < 7; i++)
+    for(Int_t i=0; i < fgkNmcTypes; i++)
     {
       fhPtTriggerMC[i]  = new TH1F(Form("hPtTrigger_MC%s",nameMC[i].Data()),
                                    Form("#it{p}_{T} distribution of trigger particles, trigger origin is %s",nameMC[i].Data()),
@@ -1212,6 +1303,30 @@ TList *  AliAnaParticleHadronCorrelation::GetCreateOutputObjects()
     }
   }
   
+  if(fDecayTrigger)
+  {
+    for(Int_t ibit = 0; ibit < fNDecayBits; ibit++)
+    {
+      fhPtDecayTrigger[ibit]  = new TH1F(Form("hPtDecayTrigger_bit%d",fDecayBits[ibit]),
+                                         Form("#it{p}_{T} distribution of trigger particles, decay Bit %d",fDecayBits[ibit]),
+                                         nptbins,ptmin,ptmax);
+      fhPtDecayTrigger[ibit]->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
+      outputContainer->Add(fhPtDecayTrigger[ibit]);
+      
+      if(IsDataMC())
+      {
+        for(Int_t i=0; i < fgkNmcTypes; i++)
+        {
+          fhPtDecayTriggerMC[ibit][i]  = new TH1F(Form("hPtDecayTrigger_bit%d_MC%s",fDecayBits[ibit], nameMC[i].Data()),
+                                            Form("#it{p}_{T} distribution of trigger particles, decay Bit %d, trigger origin is %s",fDecayBits[ibit], nameMC[i].Data()),
+                                            nptbins,ptmin,ptmax);
+          fhPtDecayTriggerMC[ibit][i]->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
+          outputContainer->Add(fhPtDecayTriggerMC[ibit][i]);
+        }
+      }
+    }
+  }
+  
   if(fCorrelVzBin)
   {
     fhPtTriggerVzBin  = new TH2F("hPtTriggerVzBin","#it{p}_{T} distribution of trigger particles vs vz bin", nptbins,ptmin,ptmax,GetNZvertBin(),0,GetNZvertBin());
@@ -1233,7 +1348,7 @@ TList *  AliAnaParticleHadronCorrelation::GetCreateOutputObjects()
   fhEtaTrigger->SetYTitle("#eta ");
   outputContainer->Add(fhEtaTrigger);
   
-  if(fFillHighMultHistograms)
+  if(IsHighMultiplicityAnalysisOn())
   {
     fhPtTriggerCentrality   = new TH2F("hPtTriggerCentrality","Trigger particle #it{p}_{T} vs centrality",nptbins,ptmin,ptmax,100,0.,100) ;
     fhPtTriggerCentrality->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
@@ -1252,16 +1367,28 @@ TList *  AliAnaParticleHadronCorrelation::GetCreateOutputObjects()
   }
 
   // Leading hadron in oposite side
-  if(fSelectLeadingHadronAngle)
+  if(fFillLeadHadOppositeHisto)
   {
     fhPtLeadingOppositeHadron  = new TH2F("hPtTriggerPtLeadingOppositeHadron","Leading hadron opposite to trigger vs trigger #it{p}_{T}",
                                           nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
     fhPtLeadingOppositeHadron->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
     fhPtLeadingOppositeHadron->SetYTitle("#it{p}_{T}^{lead hadron} (GeV/#it{c})");
     outputContainer->Add(fhPtLeadingOppositeHadron);
+
+    fhPtNoLeadingOppositeHadron  = new TH1F("hPtTriggerNoLeadingOppositeHadron","No Leading hadron opposite to trigger #it{p}_{T}",
+                                          nptbins,ptmin,ptmax);
+    fhPtNoLeadingOppositeHadron->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
+    outputContainer->Add(fhPtNoLeadingOppositeHadron);
+
+    fhEtaPhiNoLeadingOppositeHadron  = new TH2F("hEtaPhiTriggerNoLeadingOppositeHadron","No Leading hadron opposite to trigger #eta:#phi",
+                                            netabins,etamin,etamax,nphibins,phimin,phimax);
+    fhEtaPhiNoLeadingOppositeHadron->SetXTitle("#eta");
+    fhEtaPhiNoLeadingOppositeHadron->SetYTitle("#phi");
+    outputContainer->Add(fhEtaPhiNoLeadingOppositeHadron);
+
     
     fhPtDiffPhiLeadingOppositeHadron  = new TH2F("hPtTriggerDiffPhiTriggerLeadingOppositeHadron","#phi_{trigger}-#phi_{leading opposite hadron} vs #it{p}_{T}^{trig}",
-                                                 nptbins,ptmin,ptmax,ndeltaphibins ,deltaphimin,deltaphimax);
+                                                 nptbins,ptmin,ptmax,ndeltaphibins,deltaphimin,deltaphimax);
     fhPtDiffPhiLeadingOppositeHadron->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
     fhPtDiffPhiLeadingOppositeHadron->SetYTitle("#phi_{trigger}-#phi_{leading opposite hadron} (rad)");
     outputContainer->Add(fhPtDiffPhiLeadingOppositeHadron);
@@ -1417,7 +1544,7 @@ TList *  AliAnaParticleHadronCorrelation::GetCreateOutputObjects()
   
   if(IsDataMC())
   {
-    for(Int_t i=0; i < 7; i++)
+    for(Int_t i=0; i < fgkNmcTypes; i++)
     {
       
       fhDeltaPhiChargedMC[i]  = new TH2F(Form("hDeltaPhiCharged_MC%s",nameMC[i].Data()),
@@ -1588,7 +1715,7 @@ TList *  AliAnaParticleHadronCorrelation::GetCreateOutputObjects()
     outputContainer->Add(fhPtHbpZTUeLeftCharged) ;
   }
   
-  if(fFillPileUpHistograms)
+  if(IsPileUpAnalysisOn())
   {
     fhDeltaPhiChargedOtherBC  = new TH2F
     ("hDeltaPhiChargedOtherBC","#phi_{trigger} - #phi_{h^{#pm}} vs #it{p}_{T trigger}, track BC!=0",
@@ -1817,7 +1944,7 @@ TList *  AliAnaParticleHadronCorrelation::GetCreateOutputObjects()
     }
   }
 
-  if(fFillHighMultHistograms)
+  if(IsHighMultiplicityAnalysisOn())
   {
     Int_t nMultiBins = GetNCentrBin();
     fhDeltaPhiChargedMult = new TH2F*[nMultiBins] ;
@@ -1902,12 +2029,11 @@ TList *  AliAnaParticleHadronCorrelation::GetCreateOutputObjects()
   if(fFillBradHisto)
     fhDeltaPhiBradAssocPtBin = new TH2F*[fNAssocPtBins*nz];
   
-  if(fPi0Trigger || fDecayTrigger)
-  {
-    fhDeltaPhiAssocPtBin       = new TH2F*[fNAssocPtBins*nz];
-    fhDeltaPhiAssocPtBinDEta08 = new TH2F*[fNAssocPtBins*nz];
-    fhDeltaPhiDecayChargedAssocPtBin = new TH2F*[fNAssocPtBins*nz];
-  }
+
+  fhDeltaPhiAssocPtBin       = new TH2F*[fNAssocPtBins*nz];
+  if(fFillEtaGapsHisto)fhDeltaPhiAssocPtBinDEta08       = new TH2F*[fNAssocPtBins*nz];
+  if(fDecayTrigger)    fhDeltaPhiDecayChargedAssocPtBin = new TH2F*[fNAssocPtBins*nz];
+  
   
   if(fHMPIDCorrelation)
   {
@@ -1962,10 +2088,10 @@ TList *  AliAnaParticleHadronCorrelation::GetCreateOutputObjects()
         outputContainer->Add(fhDeltaPhiAssocPtBinDEta0[bin]) ;
       }
       
-      if(fPi0Trigger || fDecayTrigger)
+      if(fDecayTrigger)
       {
-        fhDeltaPhiDecayChargedAssocPtBin[bin] = new TH2F(Form("hDeltaPhiPtDecayChargedAssocPt%2.1f_%2.1f%s", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],sz.Data()),
-                                                         Form("#Delta #phi vs #it{p}_{T trigger} tagged as decay for associated #it{p}_{T} bin [%2.1f,%2.1f]%s", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],tz.Data()),
+        fhDeltaPhiDecayChargedAssocPtBin[bin] = new TH2F(Form("hDeltaPhiPtDecayChargedAssocPt%2.1f_%2.1f%s_bit%d", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],sz.Data(),fDecayBits[0]),
+                                                         Form("#Delta #phi vs #it{p}_{T trigger} tagged as decay for associated #it{p}_{T} bin [%2.1f,%2.1f]%s, Bit %d", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],tz.Data(),fDecayBits[0]),
                                                          nptbins, ptmin, ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
         fhDeltaPhiDecayChargedAssocPtBin[bin]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
         fhDeltaPhiDecayChargedAssocPtBin[bin]->SetYTitle("#Delta #phi (rad)");
@@ -2049,39 +2175,68 @@ TList *  AliAnaParticleHadronCorrelation::GetCreateOutputObjects()
     }
   }
 
-  if(fPi0Trigger || fDecayTrigger)
+  if(fPi0Trigger)
   {
-    if(fPi0Trigger)
-    {
-      fhPtPi0DecayRatio  = new TH2F
-      ("hPtPi0DecayRatio","#it{p}_{T} of #pi^{0} and the ratio of pt for two decay",
-       nptbins,ptmin,ptmax, 100,0.,2.);
-      fhPtPi0DecayRatio->SetXTitle("#it{p}_{T}^{#pi^{0}} (GeV/#it{c})");
-      fhPtPi0DecayRatio->SetYTitle("#it{p}_{T}^{Decay}/#it{p}_{T}^{#pi^{0}}");
-      outputContainer->Add(fhPtPi0DecayRatio) ;
-    }
-    
-    fhDeltaPhiDecayCharged  = new TH2F
-    ("hDeltaPhiDecayCharged","#phi_{Decay} - #phi_{h^{#pm}} vs #it{p}_{T Decay}",
+    fhPtPi0DecayRatio  = new TH2F
+    ("hPtPi0DecayRatio","#it{p}_{T} of #pi^{0} and the ratio of pt for two decay",
+     nptbins,ptmin,ptmax, 100,0.,2.);
+    fhPtPi0DecayRatio->SetXTitle("#it{p}_{T}^{#pi^{0}} (GeV/#it{c})");
+    fhPtPi0DecayRatio->SetYTitle("#it{p}_{T}^{Decay}/#it{p}_{T}^{#pi^{0}}");
+    outputContainer->Add(fhPtPi0DecayRatio) ;
+    
+    fhDeltaPhiPi0DecayCharged  = new TH2F
+    ("hDeltaPhiPi0DecayCharged","#phi_{Decay} - #phi_{h^{#pm}} vs #it{p}_{T Decay}",
      nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
-    fhDeltaPhiDecayCharged->SetYTitle("#Delta #phi (rad)");
-    fhDeltaPhiDecayCharged->SetXTitle("#it{p}_{T Decay} (GeV/#it{c})");
+    fhDeltaPhiPi0DecayCharged->SetYTitle("#Delta #phi (rad)");
+    fhDeltaPhiPi0DecayCharged->SetXTitle("#it{p}_{T Decay} (GeV/#it{c})");
     
-    fhXEDecayCharged  =
-    new TH2F("hXEDecayCharged","#it{x}_{#it{E}}  Decay",
+    fhXEPi0DecayCharged  =
+    new TH2F("hXEPi0DecayCharged","#it{x}_{#it{E}}  Decay",
              nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
-    fhXEDecayCharged->SetYTitle("#it{x}_{#it{E}}");
-    fhXEDecayCharged->SetXTitle("#it{p}_{T decay} (GeV/#it{c})");
+    fhXEPi0DecayCharged->SetYTitle("#it{x}_{#it{E}}");
+    fhXEPi0DecayCharged->SetXTitle("#it{p}_{T decay} (GeV/#it{c})");
     
-    fhZTDecayCharged  =
-    new TH2F("hZTDecayCharged","#it{z}_{trigger h^{#pm}} = #it{p}_{T h^{#pm}} / #it{p}_{T Decay}",
+    fhZTPi0DecayCharged  =
+    new TH2F("hZTPi0DecayCharged","#it{z}_{trigger h^{#pm}} = #it{p}_{T h^{#pm}} / #it{p}_{T Decay}",
              nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
-    fhZTDecayCharged->SetYTitle("#it{z}_{decay h^{#pm}}");
-    fhZTDecayCharged->SetXTitle("#it{p}_{T decay} (GeV/#it{c})");
+    fhZTPi0DecayCharged->SetYTitle("#it{z}_{decay h^{#pm}}");
+    fhZTPi0DecayCharged->SetXTitle("#it{p}_{T decay} (GeV/#it{c})");
     
-    outputContainer->Add(fhDeltaPhiDecayCharged) ;
-    outputContainer->Add(fhXEDecayCharged) ;
-    outputContainer->Add(fhZTDecayCharged) ;
+    outputContainer->Add(fhDeltaPhiPi0DecayCharged) ;
+    outputContainer->Add(fhXEPi0DecayCharged) ;
+    outputContainer->Add(fhZTPi0DecayCharged) ;
+  }
+  
+  if(fDecayTrigger)
+  {
+    for(Int_t ibit = 0; ibit< fNDecayBits; ibit++)
+    {
+      fhDeltaPhiDecayCharged[ibit]  = new TH2F
+      (Form("hDeltaPhiDecayCharged_bit%d",fDecayBits[ibit]),
+       Form("#phi_{Decay} - #phi_{h^{#pm}} vs #it{p}_{T Decay}, Bit %d",fDecayBits[ibit]),
+       nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
+      fhDeltaPhiDecayCharged[ibit]->SetYTitle("#Delta #phi (rad)");
+      fhDeltaPhiDecayCharged[ibit]->SetXTitle("#it{p}_{T Decay} (GeV/#it{c})");
+      
+      fhXEDecayCharged[ibit]  =
+      new TH2F(Form("hXEDecayCharged_bit%d",fDecayBits[ibit]),
+               Form("#it{x}_{#it{E}}  Decay, Bit %d",fDecayBits[ibit]),
+               nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
+      fhXEDecayCharged[ibit]->SetYTitle("#it{x}_{#it{E}}");
+      fhXEDecayCharged[ibit]->SetXTitle("#it{p}_{T decay} (GeV/#it{c})");
+      
+      fhZTDecayCharged[ibit]  =
+      new TH2F(Form("hZTDecayCharged_bit%d",fDecayBits[ibit]),
+               Form("#it{z}_{trigger h^{#pm}} = #it{p}_{T h^{#pm}} / #it{p}_{T Decay}, Bit %d",fDecayBits[ibit]),
+               nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
+      fhZTDecayCharged[ibit]->SetYTitle("#it{z}_{decay h^{#pm}}");
+      fhZTDecayCharged[ibit]->SetXTitle("#it{p}_{T decay} (GeV/#it{c})");
+      
+      outputContainer->Add(fhDeltaPhiDecayCharged[ibit]) ;
+      outputContainer->Add(fhXEDecayCharged[ibit]) ;
+      outputContainer->Add(fhZTDecayCharged[ibit]) ;
+    }
   }
   
   //Correlation with neutral hadrons
@@ -2232,30 +2387,29 @@ TList *  AliAnaParticleHadronCorrelation::GetCreateOutputObjects()
       outputContainer->Add(fhPtHbpZTUeLeftNeutral) ;
     }
     
-    if(fPi0Trigger || fDecayTrigger)
+    if(fPi0Trigger)
     {
-      fhDeltaPhiDecayNeutral  = new TH2F
-      ("hDeltaPhiDecayNeutral","#phi_{Decay} - #phi_{h^{0}} vs #it{p}_{T Decay}",
+      fhDeltaPhiPi0DecayNeutral  = new TH2F
+      ("hDeltaPhiPi0DecayNeutral","#phi_{Decay} - #phi_{h^{0}} vs #it{p}_{T Decay}",
        nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
-      fhDeltaPhiDecayNeutral->SetYTitle("#Delta #phi (rad)");
-      fhDeltaPhiDecayNeutral->SetXTitle("#it{p}_{T Decay} (GeV/#it{c})");
+      fhDeltaPhiPi0DecayNeutral->SetYTitle("#Delta #phi (rad)");
+      fhDeltaPhiPi0DecayNeutral->SetXTitle("#it{p}_{T Decay} (GeV/#it{c})");
       
-      fhXEDecayNeutral  =
-      new TH2F("hXEDecayNeutral","#it{x}_{#it{E}} for decay trigger",
+      fhXEPi0DecayNeutral  =
+      new TH2F("hXEPi0DecayNeutral","#it{x}_{#it{E}} for decay trigger",
                nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
-      fhXEDecayNeutral->SetYTitle("#it{x}_{#it{E}}");
-      fhXEDecayNeutral->SetXTitle("#it{p}_{T decay}");
+      fhXEPi0DecayNeutral->SetYTitle("#it{x}_{#it{E}}");
+      fhXEPi0DecayNeutral->SetXTitle("#it{p}_{T decay}");
       
-      fhZTDecayNeutral  =
-      new TH2F("hZTDecayNeutral","#it{z}_{trigger h^{0}} = #it{p}_{T h^{0}} / #it{p}_{T Decay}",
+      fhZTPi0DecayNeutral  =
+      new TH2F("hZTPi0DecayNeutral","#it{z}_{trigger h^{0}} = #it{p}_{T h^{0}} / #it{p}_{T Decay}",
                nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
-      fhZTDecayNeutral->SetYTitle("#it{z}_{h^{0}}");
-      fhZTDecayNeutral->SetXTitle("#it{p}_{T decay}");
-      
-      outputContainer->Add(fhDeltaPhiDecayNeutral) ;
-      outputContainer->Add(fhXEDecayNeutral) ;
-      outputContainer->Add(fhZTDecayNeutral) ;
+      fhZTPi0DecayNeutral->SetYTitle("#it{z}_{h^{0}}");
+      fhZTPi0DecayNeutral->SetXTitle("#it{p}_{T decay}");
       
+      outputContainer->Add(fhDeltaPhiPi0DecayNeutral) ;
+      outputContainer->Add(fhXEPi0DecayNeutral) ;
+      outputContainer->Add(fhZTPi0DecayNeutral) ;
     }
   }//Correlation with neutral hadrons
   
@@ -2409,29 +2563,29 @@ TList *  AliAnaParticleHadronCorrelation::GetCreateOutputObjects()
       fhMCUePart[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
       
       fhMCPtXEUeCharged[i]  =
-      new TH2F(Form("hMCPtXEUeCharged%s",right.Data()),
-               Form("MC %s: #it{x}_{#it{E}} with charged hadrons, Underlying Event",nameMC[i].Data()),
+      new TH2F(Form("hMCPtXEUeCharged%s_%s",right.Data(),nameMC[i].Data()),
+               Form("MC %s: #it{x}_{#it{E}} with charged hadrons, Underlying Event %s",nameMC[i].Data(),right.Data()),
                nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
       fhMCPtXEUeCharged[i]->SetYTitle("#it{x}_{#it{E}}");
       fhMCPtXEUeCharged[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
       
       fhMCPtHbpXEUeCharged[i] =
-      new TH2F(Form("hMCPtHbpXEUeCharged%s",right.Data()),
-               Form("MC %s: #xi = ln(1/#it{x}_{#it{E}}) with charged hadrons, Underlying Event",nameMC[i].Data()),
+      new TH2F(Form("hMCPtHbpXEUeCharged%s_%s",right.Data(),nameMC[i].Data()),
+               Form("MC %s: #xi = ln(1/#it{x}_{#it{E}}) with charged hadrons, Underlying Event %s",nameMC[i].Data(),right.Data()),
                nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
       fhMCPtHbpXEUeCharged[i]->SetYTitle("ln(1/#it{x}_{#it{E}})");
       fhMCPtHbpXEUeCharged[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
       
       fhMCPtZTUeCharged[i] =
-      new TH2F(Form("hMCPtZTUeCharged%s",right.Data()),
-               Form("MC %s: #it{z}_{T} with charged hadrons, Underlying Event",nameMC[i].Data()),
+      new TH2F(Form("hMCPtZTUeCharged%s_%s",right.Data(),nameMC[i].Data()),
+               Form("MC %s: #it{z}_{T} with charged hadrons, Underlying Event %s",nameMC[i].Data(),right.Data()),
                nptbins,ptmin,ptmax,nxeztbins,xeztmin,xeztmax);
       fhMCPtZTUeCharged[i]->SetYTitle("#it{z}_{T}");
       fhMCPtZTUeCharged[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
       
       fhMCPtHbpZTUeCharged[i] =
-      new TH2F(Form("hMCPtHbpZTUeCharged%s",right.Data()),
-               Form("MC %s: #xi = ln(1/#it{z}_{T}) with charged hadrons, Underlying Event",nameMC[i].Data()),
+      new TH2F(Form("hMCPtHbpZTUeCharged%s_%s",right.Data(),nameMC[i].Data()),
+               Form("MC %s: #xi = ln(1/#it{z}_{T}) with charged hadrons, Underlying Event %s",nameMC[i].Data(),right.Data()),
                nptbins,ptmin,ptmax,nhbpbins,hbpmin,hbpmax);
       fhMCPtHbpZTUeCharged[i]->SetYTitle("ln(1/#it{z}_{T})");
       fhMCPtHbpZTUeCharged[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
@@ -2700,35 +2854,30 @@ TList *  AliAnaParticleHadronCorrelation::GetCreateOutputObjects()
   
 }
 
-//_________________________________________________________________________________________________
-Bool_t AliAnaParticleHadronCorrelation::GetDecayPhotonMomentum(AliAODPWG4Particle* trigger,
-                                                               TLorentzVector & mom1,
-                                                               TLorentzVector & mom2)
+//_____________________________________________________________________________________________________________________
+Bool_t AliAnaParticleHadronCorrelation::GetDecayPhotonMomentum(Int_t indexPhoton1, Int_t indexPhoton2, Int_t idetector)
 {
   // 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);
-  
   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()  ;
+  if(idetector==kEMCAL) 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(photon->GetID()==indexPhoton2) photon->GetMomentum(mom1,GetVertex(0)) ;
+    if(photon->GetID()==indexPhoton1) photon->GetMomentum(fDecayMom1,GetVertex(0)) ;
+    if(photon->GetID()==indexPhoton2) photon->GetMomentum(fDecayMom2,GetVertex(0)) ;
     
-    if(GetDebug() > 1)printf("AliAnaParticleHadronCorrelation::GetDecayPhotonMomentum() - Photon1 = %f, Photon2 = %f \n", mom1.Pt(), mom2.Pt());
+    if(GetDebug() > 1)printf("AliAnaParticleHadronCorrelation::GetDecayPhotonMomentum() - Photon1 = %f, Photon2 = %f \n", fDecayMom1.Pt(), fDecayMom2.Pt());
     
   } //cluster loop        
   
@@ -2816,6 +2965,7 @@ void AliAnaParticleHadronCorrelation::InitParameters()
   fM02MaxCut   = -1 ;
   
   fSelectLeadingHadronAngle = kFALSE;
+  fFillLeadHadOppositeHisto = kFALSE;
   fMinLeadHadPhi = 150*TMath::DegToRad();
   fMaxLeadHadPhi = 210*TMath::DegToRad();
 
@@ -2825,6 +2975,11 @@ void AliAnaParticleHadronCorrelation::InitParameters()
   fMCGenTypeMin = 0;
   fMCGenTypeMax = 6;
   
+  fNDecayBits = 1;
+  fDecayBits[0] = AliNeutralMesonSelection::kPi0;
+  fDecayBits[1] = AliNeutralMesonSelection::kEta;
+  fDecayBits[2] = AliNeutralMesonSelection::kPi0Side;
+  fDecayBits[3] = AliNeutralMesonSelection::kEtaSide;
 }
 
 //_________________________________________________________________________
@@ -2869,7 +3024,6 @@ Bool_t AliAnaParticleHadronCorrelation::IsTriggerTheEventLeadingParticle()
   
   // Compare if it is the leading of all tracks
   
-  TVector3 p3;
   for(Int_t ipr = 0;ipr < GetCTSTracks()->GetEntriesFast() ; ipr ++ )
   {
     AliVTrack * track = (AliVTrack *) (GetCTSTracks()->At(ipr)) ;
@@ -2877,16 +3031,19 @@ Bool_t AliAnaParticleHadronCorrelation::IsTriggerTheEventLeadingParticle()
     if(track->GetID() == pLeading->GetTrackLabel(0) || track->GetID() == pLeading->GetTrackLabel(1) ||
        track->GetID() == pLeading->GetTrackLabel(2) || track->GetID() == pLeading->GetTrackLabel(3)   ) continue ;
     
-    Double_t mom[3] = {track->Px(),track->Py(),track->Pz()};
-    p3.SetXYZ(mom[0],mom[1],mom[2]);
-    Float_t pt   = p3.Pt();
-    Float_t phi  = p3.Phi() ;
+    fTrackVector.SetXYZ(track->Px(),track->Py(),track->Pz());
+    Float_t pt   = fTrackVector.Pt();
+    Float_t phi  = fTrackVector.Phi() ;
     if(phi < 0) phi+=TMath::TwoPi();
     
     //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;
+      Float_t deltaPhi = phiTrig-phi;
+      if(deltaPhi <= -TMath::PiOver2()) deltaPhi+=TMath::TwoPi();
+      if(deltaPhi > 3*TMath::PiOver2()) deltaPhi-=TMath::TwoPi();
+      
+      if(pt > ptTrig && deltaPhi < TMath::PiOver2())  return kFALSE;
     }
     //jump out this event if there is any other particle with pt larger than trigger
     else
@@ -2901,24 +3058,23 @@ Bool_t AliAnaParticleHadronCorrelation::IsTriggerTheEventLeadingParticle()
   {
     // Select the calorimeter cluster list
     TObjArray * nePl = 0x0;
-    if      (pLeading->GetDetector() == "PHOS" )
+    if      (pLeading->GetDetectorTag() == kPHOS )
       nePl = GetPHOSClusters();
     else
       nePl = GetEMCALClusters();
     
     if(!nePl) return kTRUE; // Do the selection just with the tracks if no calorimeter is available.
     
-    TLorentzVector lv;
     for(Int_t ipr = 0;ipr < nePl->GetEntriesFast() ; ipr ++ )
     {
       AliVCluster * cluster = (AliVCluster *) (nePl->At(ipr)) ;
       
       if(cluster->GetID() == pLeading->GetCaloLabel(0) || cluster->GetID() == pLeading->GetCaloLabel(1) ) continue ;
       
-      cluster->GetMomentum(lv,GetVertex(0));
+      cluster->GetMomentum(fMomentum,GetVertex(0));
       
-      Float_t pt   = lv.Pt();
-      Float_t phi  = lv.Phi() ;
+      Float_t pt   = fMomentum.Pt();
+      Float_t phi  = fMomentum.Phi() ;
       if(phi < 0) phi+=TMath::TwoPi();
       
       if(IsTrackMatched(cluster,GetReader()->GetInputEvent())) continue ; // avoid charged clusters, already covered by tracks, or cluster merging with track.
@@ -2927,7 +3083,11 @@ Bool_t AliAnaParticleHadronCorrelation::IsTriggerTheEventLeadingParticle()
       // not really needed for calorimeter, unless DCal is included
       if (fMakeNearSideLeading)
       {
-        if(pt > ptTrig && TMath::Abs(phi-phiTrig) < TMath::PiOver2()) return kFALSE ;
+        Float_t deltaPhi = phiTrig-phi;
+        if(deltaPhi <= -TMath::PiOver2()) deltaPhi+=TMath::TwoPi();
+        if(deltaPhi > 3*TMath::PiOver2()) deltaPhi-=TMath::TwoPi();
+        
+        if(pt > ptTrig && deltaPhi < TMath::PiOver2()) return kFALSE ;
       }
       //jump out this event if there is any other particle with pt larger than trigger
       else
@@ -3003,7 +3163,7 @@ void  AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms()
   
   Float_t cen         = GetEventCentrality();
   Float_t ep          = GetEventPlaneAngle();
-  if(fFillHighMultHistograms) fhTriggerEventPlaneCentrality->Fill(cen,ep);
+  if(IsHighMultiplicityAnalysisOn()) fhTriggerEventPlaneCentrality->Fill(cen,ep);
 
   Int_t   mixEventBin = GetEventMixBin();
   Int_t   vzbin       = GetEventVzBin();
@@ -3032,22 +3192,25 @@ void  AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms()
     //
     Int_t clID1  = particle->GetCaloLabel(0) ;
     Int_t clID2  = particle->GetCaloLabel(1) ; // for photon clusters should not be set.
-    if( GetDebug() > 1 ) printf("%s Trigger : id1 %d, id2 %d, min %f, max %f, det %s\n",
-           GetInputAODName().Data(),clID1,clID2,fM02MinCut,fM02MaxCut,(particle->GetDetector()).Data());
+    if( GetDebug() > 1 ) printf("%s Trigger : id1 %d, id2 %d, min %f, max %f, det %d\n",
+           GetInputAODName().Data(),clID1,clID2,fM02MinCut,fM02MaxCut,particle->GetDetectorTag());
     
     if(clID1 > 0 && clID2 < 0 && fM02MaxCut > 0 && fM02MinCut > 0)
     {
-      Int_t iclus = -1;
-      TObjArray* clusters = 0x0;
-      if     (particle->GetDetector() == "EMCAL") clusters = GetEMCALClusters();
-      else if(particle->GetDetector() == "PHOS" ) clusters = GetPHOSClusters();
+//      Int_t iclus = -1;
+//      TObjArray* clusters = 0x0;
+//      if     (particle->GetDetectorTag() == kEMCAL) clusters = GetEMCALClusters();
+//      else if(particle->GetDetectorTag() == kPHOS ) clusters = GetPHOSClusters();
+//      
+//      if(clusters)
+//      {
+//        AliVCluster *cluster = FindCluster(clusters,clID1,iclus);
+//        Float_t m02 = cluster->GetM02();
+//        if(m02 > fM02MaxCut || m02 < fM02MinCut) continue ;
+//      }
       
-      if(clusters)
-      {
-        AliVCluster *cluster = FindCluster(clusters,clID1,iclus);
-        Float_t m02 = cluster->GetM02();
-        if(m02 > fM02MaxCut || m02 < fM02MinCut) continue ;
-      }
+      Float_t m02 = particle->GetM02();
+      if(m02 > fM02MaxCut || m02 < fM02MinCut) continue ;
       
       fhPtTriggerSSCut->Fill(pt);
     }
@@ -3067,7 +3230,7 @@ void  AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms()
     //
     if(IsFiducialCutOn())
     {
-      Bool_t in = GetFiducialCut()->IsInFiducialCut(*particle->Momentum(),particle->GetDetector()) ;
+      Bool_t in = GetFiducialCut()->IsInFiducialCut(particle->Eta(),particle->Phi(),particle->GetDetectorTag()) ;
       if(! in ) continue ;
     }
     
@@ -3079,9 +3242,11 @@ void  AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms()
     // Find the leading hadron in the opposite hemisphere to the triggeer
     // and accept the trigger if leading is in defined window.
     Bool_t okLeadHad = kTRUE;
-    if(fSelectLeadingHadronAngle)
+    if(fSelectLeadingHadronAngle || fFillLeadHadOppositeHisto)
+    {
       okLeadHad = FindLeadingOppositeHadronInWindow(particle);
-    if(!okLeadHad) continue;
+      if(!okLeadHad && fSelectLeadingHadronAngle) continue;
+    }
     
     //
     // Charged particles correlation
@@ -3090,10 +3255,13 @@ void  AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms()
     
     // MC
     Int_t mcIndex = -1;
+    Int_t mcTag   = particle->GetTag();
+    Bool_t lostDecayPair = kFALSE;
     if(IsDataMC())
     {
-      mcIndex = GetMCTagHistogramIndex(particle->GetTag());
-      MakeMCChargedCorrelation(particle->GetLabel(), mcIndex);
+      mcIndex = GetMCTagHistogramIndex(mcTag);
+      lostDecayPair = GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCDecayPairLost);
+      MakeMCChargedCorrelation(particle->GetLabel(), mcIndex,lostDecayPair);
     }
     
     // Do own mixed event with charged,
@@ -3114,8 +3282,33 @@ void  AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms()
     // pT of the trigger, vs trigger origin if MC
     //
     fhPtTrigger->Fill(pt);
-    if(IsDataMC() && mcIndex >=0 && mcIndex < 7)
+    if(IsDataMC() && mcIndex >=0 && mcIndex < fgkNmcTypes)
+    {
       fhPtTriggerMC[mcIndex]->Fill(pt);
+      if( lostDecayPair && mcIndex==2 )
+        fhPtTriggerMC[7]->Fill(pt);
+    }
+    
+    if(fDecayTrigger)
+    {
+      Int_t decayTag = particle->DecayTag();
+      if(decayTag < 0) decayTag = 0;
+      
+      for(Int_t ibit = 0; ibit<fNDecayBits; ibit++)
+      {
+        if(GetNeutralMesonSelection()->CheckDecayBit(decayTag,fDecayBits[ibit]))
+        {
+          fhPtDecayTrigger[ibit]->Fill(pt);
+          
+          if(IsDataMC() && mcIndex >=0 && mcIndex < fgkNmcTypes)
+          {
+            fhPtDecayTriggerMC[ibit][mcIndex]->Fill(pt);
+            if(lostDecayPair && mcIndex==2 )
+              fhPtDecayTriggerMC[ibit][7]->Fill(pt);
+          }
+        }// check bit
+      }// bit loop
+    }
     
     //
     // Acceptance of the trigger
@@ -3134,7 +3327,7 @@ void  AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms()
     if(fCorrelVzBin)
       fhPtTriggerVzBin->Fill(pt,vzbin);
     
-    if(fFillHighMultHistograms)
+    if(IsHighMultiplicityAnalysisOn())
     {
       fhPtTriggerCentrality->Fill(pt,cen);
       fhPtTriggerEventPlane->Fill(pt,ep);
@@ -3143,7 +3336,7 @@ void  AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms()
     //----------------------------------
     // Trigger particle pT vs pile-up
     
-    if(fFillPileUpHistograms)
+    if(IsPileUpAnalysisOn())
     {
       Int_t vtxBC = GetReader()->GetVertexBC();
       if(vtxBC == 0 || vtxBC==AliVTrack::kTOFBCNA)     fhPtTriggerVtxBC0->Fill(pt);
@@ -3174,17 +3367,29 @@ void  AliAnaParticleHadronCorrelation::MakeChargedCorrelation(AliAODPWG4Particle
   Float_t phiTrig = aodParticle->Phi();
   Float_t etaTrig = aodParticle->Eta();
   Float_t ptTrig  = aodParticle->Pt();  
-  Bool_t   decay  = aodParticle->IsTagged();
   Int_t    mcTag  = aodParticle->GetTag();
   Double_t bz     = GetReader()->GetInputEvent()->GetMagneticField();
 
+  
+  Int_t   decayTag = 0;
+  if(fDecayTrigger)
+  {
+    //decay = aodParticle->IsTagged();
+    decayTag = aodParticle->DecayTag();
+    if(decayTag < 0) decayTag = 0;
+//    printf("Correlation: pT %2.2f, BTag %d, Tagged %d\n",ptTrig, decayTag, aodParticle->IsTagged());
+//    printf("\t check bit Pi0 %d, Eta %d,  Pi0Side %d, EtaSide %d\n",
+//           GetNeutralMesonSelection()->CheckDecayBit(decayTag,AliNeutralMesonSelection::kPi0),
+//           GetNeutralMesonSelection()->CheckDecayBit(decayTag,AliNeutralMesonSelection::kEta),
+//           GetNeutralMesonSelection()->CheckDecayBit(decayTag,AliNeutralMesonSelection::kPi0Side),
+//           GetNeutralMesonSelection()->CheckDecayBit(decayTag,AliNeutralMesonSelection::kEtaSide));
+  }
+
   Float_t pt       = -100. ;
   Float_t phi      = -100. ;
   Float_t eta      = -100. ;
   Float_t deltaPhi = -100. ;
   
-  TVector3 p3;  
-  TLorentzVector photonMom ;   
   TObjArray * reftracks = 0x0;
   Int_t nrefs           = 0;
   
@@ -3200,20 +3405,23 @@ void  AliAnaParticleHadronCorrelation::MakeChargedCorrelation(AliAODPWG4Particle
     evtIndex13 = GetMixedEvent()->EventIndex(aodParticle->GetTrackLabel(0)) ;
   }
   
+  // Track multiplicity or cent bin
+  Int_t cenbin = 0;
+  if(IsHighMultiplicityAnalysisOn()) cenbin = GetEventCentralityBin();
+
   //
-  // In case of pi0/eta trigger, we may want to check their decay 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 )
   {
-    decayFound = GetDecayPhotonMomentum(aodParticle,decayMom1, decayMom2);
+    decayFound = GetDecayPhotonMomentum(aodParticle->GetCaloLabel(0),aodParticle->GetCaloLabel(1),aodParticle->GetDetectorTag());
     if(decayFound)
     {
-      fhPtPi0DecayRatio->Fill(ptTrig, decayMom1.Pt()/ptTrig);
-      fhPtPi0DecayRatio->Fill(ptTrig, decayMom2.Pt()/ptTrig);
+      fhPtPi0DecayRatio->Fill(ptTrig, fDecayMom1.Pt()/ptTrig);
+      fhPtPi0DecayRatio->Fill(ptTrig, fDecayMom2.Pt()/ptTrig);
     }
   }
   
@@ -3225,11 +3433,10 @@ void  AliAnaParticleHadronCorrelation::MakeChargedCorrelation(AliAODPWG4Particle
   {
     AliVTrack * track = (AliVTrack *) (GetCTSTracks()->At(ipr)) ;
     
-    Double_t mom[3] = {track->Px(),track->Py(),track->Pz()};
-    p3.SetXYZ(mom[0],mom[1],mom[2]);
-    pt   = p3.Pt();
-    eta  = p3.Eta();
-    phi  = p3.Phi() ;
+    fTrackVector.SetXYZ(track->Px(),track->Py(),track->Pz());
+    pt   = fTrackVector.Pt();
+    eta  = fTrackVector.Eta();
+    phi  = fTrackVector.Phi() ;
     if(phi < 0) phi+=TMath::TwoPi();
     
     //Select only hadrons in pt range
@@ -3297,9 +3504,6 @@ void  AliAnaParticleHadronCorrelation::MakeChargedCorrelation(AliAODPWG4Particle
     if     (okTOF && trackBC!=0) outTOF = 1;
     else if(okTOF && trackBC==0) outTOF = 0;
     
-    // Track multiplicity or cent bin
-    Int_t cenbin = GetEventCentralityBin(); // combine with vz assoc bin???
-
     //----------------
     // Fill Histograms
 
@@ -3318,7 +3522,7 @@ void  AliAnaParticleHadronCorrelation::MakeChargedCorrelation(AliAODPWG4Particle
     if(deltaPhi > 3*TMath::PiOver2()) deltaPhi-=TMath::TwoPi();
 
     FillChargedAngularCorrelationHistograms(pt,  ptTrig,  bin, phi, phiTrig,  deltaPhi,
-                                            eta, etaTrig, decay, track->GetHMPIDsignal(),
+                                            eta, etaTrig, decayTag, track->GetHMPIDsignal(),
                                             outTOF, cenbin, mcTag);
     
     //
@@ -3330,7 +3534,7 @@ void  AliAnaParticleHadronCorrelation::MakeChargedCorrelation(AliAODPWG4Particle
     //
     if  ( (deltaPhi > fDeltaPhiMinCut)   && (deltaPhi < fDeltaPhiMaxCut)   )
       FillChargedMomentumImbalanceHistograms(ptTrig, pt, deltaPhi, cenbin, track->Charge(),
-                                             assocBin, decay, outTOF, mcTag);
+                                             assocBin, decayTag, outTOF, mcTag);
     
     //
     // Underlying event, right side, default case
@@ -3347,7 +3551,7 @@ void  AliAnaParticleHadronCorrelation::MakeChargedCorrelation(AliAODPWG4Particle
     
     //
     if(fPi0Trigger && decayFound)
-      FillDecayPhotonCorrelationHistograms(pt, phi, decayMom1,decayMom2, kTRUE) ;
+      FillDecayPhotonCorrelationHistograms(pt, phi, kTRUE) ;
     
     //
     // Add track reference to array
@@ -3470,11 +3674,13 @@ void AliAnaParticleHadronCorrelation::MakeChargedMixCorrelation(AliAODPWG4Partic
     //
     if( OnlyIsolated() )
     {
-      Int_t n=0; Int_t nfrac = 0; Bool_t isolated = kFALSE; Float_t coneptsum = 0;
+      Int_t   n=0, nfrac = 0;
+      Bool_t  isolated = kFALSE;
+      Float_t coneptsum = 0, coneptlead = 0;
       GetIsolationCut()->MakeIsolationCut(bgTracks,bgCalo,
                                           GetReader(), GetCaloPID(),
                                           kFALSE, aodParticle, "",
-                                          n,nfrac,coneptsum, isolated);
+                                          n,nfrac,coneptsum,coneptlead,isolated);
       
       //printf("AliAnaParticleHadronCorrelation::MakeChargedMixCorrelation() - Isolated? %d - cone %f, ptthres %f",
       //       isolated,GetIsolationCut()->GetConeSize(),GetIsolationCut()->GetPtThreshold());
@@ -3502,7 +3708,11 @@ void AliAnaParticleHadronCorrelation::MakeChargedMixCorrelation(AliAODPWG4Partic
         
         if (fMakeNearSideLeading)
         {
-          if(ptAssoc > ptTrig && TMath::Abs(phiAssoc-phiTrig) < TMath::PiOver2())  
+          deltaPhi = phiTrig-phiAssoc;
+          if(deltaPhi <= -TMath::PiOver2()) deltaPhi+=TMath::TwoPi();
+          if(deltaPhi > 3*TMath::PiOver2()) deltaPhi-=TMath::TwoPi();
+          
+          if(ptAssoc > ptTrig && deltaPhi < TMath::PiOver2())
           {
             leading = kFALSE;
             break;
@@ -3525,7 +3735,6 @@ void AliAnaParticleHadronCorrelation::MakeChargedMixCorrelation(AliAODPWG4Partic
       if(neutralMix && fCheckLeadingWithNeutralClusters && bgCalo)
       {
         Int_t nClusters=bgCalo->GetEntriesFast();
-        TLorentzVector mom ;
         for(Int_t jlead = 0;jlead <nClusters; jlead++ )
         {
           AliAODPWG4Particle *cluster= (AliAODPWG4Particle*) bgCalo->At(jlead) ;
@@ -3536,7 +3745,11 @@ void AliAnaParticleHadronCorrelation::MakeChargedMixCorrelation(AliAODPWG4Partic
           
           if (fMakeNearSideLeading)
           {
-            if(ptAssoc > ptTrig && TMath::Abs(phiAssoc-phiTrig) < TMath::PiOver2())
+            deltaPhi = phiTrig-phiAssoc;
+            if(deltaPhi <= -TMath::PiOver2()) deltaPhi+=TMath::TwoPi();
+            if(deltaPhi > 3*TMath::PiOver2()) deltaPhi-=TMath::TwoPi();
+            
+            if(ptAssoc > ptTrig && deltaPhi < TMath::PiOver2())
             {
               leading = kFALSE;
               break;
@@ -3702,15 +3915,12 @@ void AliAnaParticleHadronCorrelation::MakeNeutralCorrelation(AliAODPWG4ParticleC
   Float_t etaTrig = aodParticle->Eta();
   Float_t deltaPhi= -100. ;
   Float_t deltaEta= -100. ;
-
-  TLorentzVector photonMom ;
        
   // 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) decayFound = GetDecayPhotonMomentum(aodParticle,decayMom1, decayMom2);
+  if(fPi0Trigger) decayFound = GetDecayPhotonMomentum(aodParticle->GetCaloLabel(0),aodParticle->GetCaloLabel(1),aodParticle->GetDetectorTag());
   
   TObjArray * refpi0 = 0x0;
   Int_t nrefs        = 0;
@@ -3811,7 +4021,7 @@ void AliAnaParticleHadronCorrelation::MakeNeutralCorrelation(AliAODPWG4ParticleC
     // Decay photon correlations
     //
     if(fPi0Trigger && decayFound)
-      FillDecayPhotonCorrelationHistograms(pt, phi, decayMom1,decayMom2,kFALSE) ;
+      FillDecayPhotonCorrelationHistograms(pt, phi, kFALSE) ;
     
     if(fFillAODWithReferences)
     {
@@ -3837,8 +4047,8 @@ void AliAnaParticleHadronCorrelation::MakeNeutralCorrelation(AliAODPWG4ParticleC
   }
 }
   
-//____________________________________________________________________________________________
-void  AliAnaParticleHadronCorrelation::MakeMCChargedCorrelation(Int_t label, Int_t histoIndex)
+//__________________________________________________________________________________________________________________
+void  AliAnaParticleHadronCorrelation::MakeMCChargedCorrelation(Int_t label, Int_t histoIndex, Bool_t lostDecayPair)
 {
   // Charged Hadron Correlation Analysis with MC information
   
@@ -3906,7 +4116,6 @@ void  AliAnaParticleHadronCorrelation::MakeMCChargedCorrelation(Int_t label, Int
     for (iParticle = 0 ; iParticle <  nTracks ; iParticle++)
     {
       TParticle * particle = stack->Particle(iParticle);
-      TLorentzVector momentum;
       
       //keep only final state particles
       if( particle->GetStatusCode() != 1 ) continue ;
@@ -3916,11 +4125,11 @@ void  AliAnaParticleHadronCorrelation::MakeMCChargedCorrelation(Int_t label, Int
       Int_t charge = (Int_t) TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
       if(charge == 0) continue;
       
-      particle->Momentum(momentum);
+      particle->Momentum(fMomentum);
       
       //Particles in CTS acceptance, make sure to use the same selection as in the reader
-      Bool_t inCTS =  GetReader()->GetFiducialCut()->IsInFiducialCut(momentum,"CTS");
-      //printf("Accepted? %d; pt %2.2f, eta %2.2f, phi %2.2f\n",inCTS,momentum.Pt(),momentum.Eta(),momentum.Phi()*TMath::RadToDeg());
+      Bool_t inCTS =  GetReader()->GetFiducialCut()->IsInFiducialCut(fMomentum.Eta(),fMomentum.Phi(),kCTS);
+      //printf("Accepted? %d; pt %2.2f, eta %2.2f, phi %2.2f\n",inCTS,fMomentum.Pt(),fMomentum.Eta(),fMomentum.Phi()*TMath::RadToDeg());
       if( !inCTS ) continue;
       
       // Remove conversions
@@ -3931,7 +4140,7 @@ void  AliAnaParticleHadronCorrelation::MakeMCChargedCorrelation(Int_t label, Int
       Float_t phi = particle->Phi();
       if(phi < 0) phi+=TMath::TwoPi();
       
-      Bool_t lead = FillChargedMCCorrelationHistograms(particle->Pt(),phi,particle->Eta(),ptprim,phiprim,etaprim,histoIndex);
+      Bool_t lead = FillChargedMCCorrelationHistograms(particle->Pt(),phi,particle->Eta(),ptprim,phiprim,etaprim,histoIndex,lostDecayPair);
       if(!lead) leadTrig = kFALSE;
       //if ( !lead && (fMakeAbsoluteLeading || fMakeNearSideLeading) ) return;
       
@@ -3978,11 +4187,11 @@ void  AliAnaParticleHadronCorrelation::MakeMCChargedCorrelation(Int_t label, Int
       
       if ( part->Charge() == 0 ) continue;
       
-      TLorentzVector momentum(part->Px(),part->Py(),part->Pz(),part->E());
+      fMomentum.SetPxPyPzE(part->Px(),part->Py(),part->Pz(),part->E());
       
       //Particles in CTS acceptance, make sure to use the same selection as in the reader
-      Bool_t inCTS =  GetReader()->GetFiducialCut()->IsInFiducialCut(momentum,"CTS");
-      //printf("Accepted? %d; pt %2.2f, eta %2.2f, phi %2.2f\n",inCTS,momentum.Pt(),momentum.Eta(),momentum.Phi()*TMath::RadToDeg());
+      Bool_t inCTS =  GetReader()->GetFiducialCut()->IsInFiducialCut(fMomentum.Eta(),fMomentum.Phi(),kCTS);
+      //printf("Accepted? %d; pt %2.2f, eta %2.2f, phi %2.2f\n",inCTS,fMomentum.Pt(),fMomentum.Eta(),fMomentum.Phi()*TMath::RadToDeg());
       if( !inCTS ) continue;
       
       // Remove conversions
@@ -3999,7 +4208,7 @@ void  AliAnaParticleHadronCorrelation::MakeMCChargedCorrelation(Int_t label, Int
       Float_t phi = part->Phi();
       if(phi < 0) phi+=TMath::TwoPi();
       
-      Bool_t lead = FillChargedMCCorrelationHistograms(part->Pt(),phi,part->Eta(),ptprim,phiprim,etaprim, histoIndex);
+      Bool_t lead = FillChargedMCCorrelationHistograms(part->Pt(),phi,part->Eta(),ptprim,phiprim,etaprim, histoIndex,lostDecayPair);
       if(!lead) leadTrig = kFALSE;
       //if ( !lead && (fMakeAbsoluteLeading || fMakeNearSideLeading)) return;
       
@@ -4013,6 +4222,13 @@ void  AliAnaParticleHadronCorrelation::MakeMCChargedCorrelation(Int_t label, Int
   fhMCPhiTrigger[histoIndex]->Fill(ptprim,phiprim);
   fhMCEtaTrigger[histoIndex]->Fill(ptprim,etaprim);
   
+  if(histoIndex==2 && lostDecayPair && 7 >= fMCGenTypeMin && 7 <= fMCGenTypeMax )
+  {
+    fhMCPtTrigger [7]->Fill(ptprim);
+    fhMCPhiTrigger[7]->Fill(ptprim,phiprim);
+    fhMCEtaTrigger[7]->Fill(ptprim,etaprim);
+  }
+  
   if(!leadTrig && (fMakeAbsoluteLeading || fMakeNearSideLeading) )
   {
     if(GetDebug() > 1)
@@ -4022,6 +4238,13 @@ void  AliAnaParticleHadronCorrelation::MakeMCChargedCorrelation(Int_t label, Int
     fhMCPtTriggerNotLeading [histoIndex]->Fill(ptprim);
     fhMCPhiTriggerNotLeading[histoIndex]->Fill(ptprim,phiprim);
     fhMCEtaTriggerNotLeading[histoIndex]->Fill(ptprim,etaprim);
+    
+    if(histoIndex==2 && lostDecayPair && 7 >= fMCGenTypeMin && 7 <= fMCGenTypeMax )
+    {
+      fhMCPtTriggerNotLeading [7]->Fill(ptprim);
+      fhMCPhiTriggerNotLeading[7]->Fill(ptprim,phiprim);
+      fhMCEtaTriggerNotLeading[7]->Fill(ptprim,etaprim);
+    }
   }
 }