diff --git a/source/processes/hadronic/models/radioactive_decay/include/G4BetaMinusDecay.hh b/source/processes/hadronic/models/radioactive_decay/include/G4BetaMinusDecay.hh
index 579129776c2..ec1085bb608 100644
--- a/source/processes/hadronic/models/radioactive_decay/include/G4BetaMinusDecay.hh
+++ b/source/processes/hadronic/models/radioactive_decay/include/G4BetaMinusDecay.hh
@@ -73,8 +73,8 @@ class G4BetaMinusDecay : public G4NuclearDecay
     const G4ParticleDefinition* fNeutrino;
 
     static const G4int npti{101};
+    G4double pdf[npti];
     G4double cdf[npti];
 };
 
 #endif
-
diff --git a/source/processes/hadronic/models/radioactive_decay/include/G4BetaPlusDecay.hh b/source/processes/hadronic/models/radioactive_decay/include/G4BetaPlusDecay.hh
index 6419510911d..fe832db9c34 100644
--- a/source/processes/hadronic/models/radioactive_decay/include/G4BetaPlusDecay.hh
+++ b/source/processes/hadronic/models/radioactive_decay/include/G4BetaPlusDecay.hh
@@ -72,8 +72,8 @@ class G4BetaPlusDecay : public G4NuclearDecay
     const G4ParticleDefinition* fNeutrino;
 
     static const G4int npti{101};
+    G4double pdf[npti];
     G4double cdf[npti];
 };
 
 #endif
-
diff --git a/source/processes/hadronic/models/radioactive_decay/include/G4BetaSpectrumSampler.hh b/source/processes/hadronic/models/radioactive_decay/include/G4BetaSpectrumSampler.hh
index d43fb77bf06..ea4603f1682 100644
--- a/source/processes/hadronic/models/radioactive_decay/include/G4BetaSpectrumSampler.hh
+++ b/source/processes/hadronic/models/radioactive_decay/include/G4BetaSpectrumSampler.hh
@@ -44,11 +44,10 @@ class G4BetaSpectrumSampler
 public:
 
   // sampling of any spectra using cumulative PDF function
-  // energy stepts are uniform starting from zero
-  // no check on array size is performed 
+  // energy steps are uniform starting from zero
+  // no check on array size is performed
   static G4double shoot(const G4int npoints, const G4double* aCDF,
-                        const G4double estep);
+                        const G4double* aPDF, const G4double estep);
 
 };
 #endif
-
diff --git a/source/processes/hadronic/models/radioactive_decay/src/G4BetaMinusDecay.cc b/source/processes/hadronic/models/radioactive_decay/src/G4BetaMinusDecay.cc
index b855e3a4095..98633ded2f6 100644
--- a/source/processes/hadronic/models/radioactive_decay/src/G4BetaMinusDecay.cc
+++ b/source/processes/hadronic/models/radioactive_decay/src/G4BetaMinusDecay.cc
@@ -99,7 +99,7 @@ G4DecayProducts* G4BetaMinusDecay::DecayIt(G4double)
   G4DecayProducts* products = new G4DecayProducts(prim);
 
   // Generate electron isotropic in angle, with energy from stored spectrum
-  const G4double eKE = eMass*G4BetaSpectrumSampler::shoot(npti, cdf, estep);
+  const G4double eKE = eMass*G4BetaSpectrumSampler::shoot(npti, cdf, pdf, estep);
 
   G4double eMomentum = std::sqrt(eKE*(eKE + 2.*eMass));
   G4ThreeVector dir = G4RandomDirection();
@@ -154,6 +154,7 @@ G4BetaMinusDecay::SetUpBetaSpectrumSampler(const G4int& daughterZ,
                                            const G4int& daughterA,
                                            const G4BetaDecayType& betaType)
 {
+  pdf[0] = 0.0;
   cdf[0] = 0.0;
 
   // Check for cases in which Q < 0
@@ -176,13 +177,18 @@ G4BetaMinusDecay::SetUpBetaSpectrumSampler(const G4int& daughterZ,
 
       // Apply shape factor for forbidden transitions
       f *= corrections.ShapeFactor(betaType, p, maxEnergy - ex);
+      pdf[i] = f;
       sum += f + f0;
       cdf[i] = sum;
       f0 = f;
     }
+    pdf[npti-1] = 0.0;
     cdf[npti-1] = sum + f0;
   } else {
-    for (G4int i = 1; i < npti; ++i) { cdf[i] = 0.0; }
+    for (G4int i = 1; i < npti; ++i) {
+      pdf[i] = 0.0;
+      cdf[i] = 0.0;
+    }
   }
 }
 
@@ -194,4 +200,3 @@ void G4BetaMinusDecay::DumpNuclearInfo()
 	 << " + " << fNeutrino->GetParticleName() << " Eemax(MeV)="
 	 << maxEnergy*eMass << " BR=" << GetBR() << "%" << G4endl;
 }
-
diff --git a/source/processes/hadronic/models/radioactive_decay/src/G4BetaPlusDecay.cc b/source/processes/hadronic/models/radioactive_decay/src/G4BetaPlusDecay.cc
index 055222bbe81..b1b638e88ec 100644
--- a/source/processes/hadronic/models/radioactive_decay/src/G4BetaPlusDecay.cc
+++ b/source/processes/hadronic/models/radioactive_decay/src/G4BetaPlusDecay.cc
@@ -101,7 +101,7 @@ G4DecayProducts* G4BetaPlusDecay::DecayIt(G4double)
   G4DecayProducts* products = new G4DecayProducts(prim);
 
   // Generate positron isotropic in angle, with energy from stored spectrum
-  const G4double eKE = eMass*G4BetaSpectrumSampler::shoot(npti, cdf, estep);
+  const G4double eKE = eMass*G4BetaSpectrumSampler::shoot(npti, cdf, pdf, estep);
 
   G4double eMomentum = std::sqrt(eKE*(eKE + 2.*eMass));
   G4ThreeVector dir = G4RandomDirection();
@@ -155,6 +155,7 @@ G4BetaPlusDecay::SetUpBetaSpectrumSampler(const G4int& daughterZ,
                                           const G4int& daughterA,
                                           const G4BetaDecayType& betaType)
 {
+  pdf[0] = 0.0;
   cdf[0] = 0.0;
 
   // Check for cases in which Q < 2Me (e.g. z67.a162) 
@@ -177,13 +178,18 @@ G4BetaPlusDecay::SetUpBetaSpectrumSampler(const G4int& daughterZ,
 
       // Apply shape factor for forbidden transitions
       f *= corrections.ShapeFactor(betaType, p, maxEnergy - ex);
+      pdf[i] = f;
       sum += f + f0;
       cdf[i] = sum;
       f0 = f;
     }
+    pdf[npti-1] = 0.0;
     cdf[npti-1] = sum + f0;
   } else {
-    for (G4int i = 0; i < npti; ++i) { cdf[i] = 0.0; }
+    for (G4int i = 1; i < npti; ++i) {
+      pdf[i] = 0.0;
+      cdf[i] = 0.0;
+    }
   }
 }
 
diff --git a/source/processes/hadronic/models/radioactive_decay/src/G4BetaSpectrumSampler.cc b/source/processes/hadronic/models/radioactive_decay/src/G4BetaSpectrumSampler.cc
index 71175935902..00e315c729e 100644
--- a/source/processes/hadronic/models/radioactive_decay/src/G4BetaSpectrumSampler.cc
+++ b/source/processes/hadronic/models/radioactive_decay/src/G4BetaSpectrumSampler.cc
@@ -37,16 +37,42 @@
 #include "G4BetaSpectrumSampler.hh"
 #include "Randomize.hh"
 
+#include <cmath>
+
+namespace
+{
+G4double SampleSlopedLine(G4double y1, G4double y2, G4double rand)
+{
+  const G4double ysum = y1 + y2;
+  if (ysum <= 0.0) {
+    return 0.0;
+  }
+  if (std::abs(2.0 * (y2 - y1) / ysum) < 1.e-12) {
+    return rand;
+  }
+  return (std::sqrt(y1 * y1 + rand * (y2 * y2 - y1 * y1)) - y1) / (y2 - y1);
+}
+}
+
 G4double G4BetaSpectrumSampler::shoot(const G4int npoints, const G4double* aCDF,
+                                      const G4double* aPDF,
                                       const G4double estep)
 {
   G4double prob = aCDF[npoints - 1]*G4UniformRand();
   G4int i = 0;
   for (; i<npoints; ++i) { if (prob <= aCDF[i]) { break; } }
+  if (i <= 0) {
+    return 0.0;
+  }
+  if (i >= npoints) {
+    i = npoints - 1;
+  }
   const G4double p1 = (i > 0) ? aCDF[i - 1] : aCDF[0];
   const G4double p2 = aCDF[i];
   const G4double delta = p2 - p1;
-  const G4double x = (delta > 0.0) ? estep*i - estep*(p2 - prob)/delta : estep*i;
+  const G4double x = (delta > 0.0)
+                       ? estep * (i - 1 + SampleSlopedLine(aPDF[i - 1], aPDF[i],
+                                                           (prob - p1) / delta))
+                       : estep * i;
   return x;
 }
-