diff --git a/src/distribution_energy.cpp b/src/distribution_energy.cpp
index 2f8e6cf1a99..7712c0d763f 100644
--- a/src/distribution_energy.cpp
+++ b/src/distribution_energy.cpp
@@ -211,59 +211,70 @@ double ContinuousTabular::sample(double E, uint64_t* seed) const
   }
 
   double E_l_k = distribution_[l].e_out[k];
-  double p_l_k = distribution_[l].p[k];
-  double E_out = E_l_k;
-  if (distribution_[l].interpolation == Interpolation::histogram) {
-    // Histogram interpolation
-    if (p_l_k > 0.0 && k >= n_discrete) {
-      E_out = E_l_k + (r1 - c_k) / p_l_k;
-    }
-
-  } else if (distribution_[l].interpolation == Interpolation::lin_lin) {
-    // Linear-linear interpolation
-    double E_l_k1 = distribution_[l].e_out[k + 1];
-    double p_l_k1 = distribution_[l].p[k + 1];
 
-    if (E_l_k != E_l_k1) {
-      double frac = (p_l_k1 - p_l_k) / (E_l_k1 - E_l_k);
-      if (frac == 0.0) {
+  if (k < n_discrete) {
+    // Discrete case
+    return E_l_k;
+  } else {
+    // Continuous case
+    double p_l_k = distribution_[l].p[k];
+    double E_out;
+    if (distribution_[l].interpolation == Interpolation::histogram) {
+      // Histogram interpolation
+      if (p_l_k > 0.0) {
         E_out = E_l_k + (r1 - c_k) / p_l_k;
       } else {
-        E_out =
-          E_l_k +
-          (std::sqrt(std::max(0.0, p_l_k * p_l_k + 2.0 * frac * (r1 - c_k))) -
-            p_l_k) /
-            frac;
+        E_out = E_l_k;
+      }
+    } else if (distribution_[l].interpolation == Interpolation::lin_lin) {
+      // Linear-linear interpolation
+      double E_l_k1 = distribution_[l].e_out[k + 1];
+      double p_l_k1 = distribution_[l].p[k + 1];
+
+      if (E_l_k != E_l_k1) {
+        double frac = (p_l_k1 - p_l_k) / (E_l_k1 - E_l_k);
+        if (frac == 0.0) {
+          E_out = E_l_k + (r1 - c_k) / p_l_k;
+        } else {
+          E_out =
+            E_l_k +
+            (std::sqrt(std::max(0.0, p_l_k * p_l_k + 2.0 * frac * (r1 - c_k))) -
+              p_l_k) /
+              frac;
+        }
+      } else {
+        E_out = E_l_k;
       }
+    } else {
+      throw std::runtime_error {
+        "Unexpected interpolation for continuous energy "
+        "distribution."};
     }
-  } else {
-    throw std::runtime_error {"Unexpected interpolation for continuous energy "
-                              "distribution."};
-  }
 
-  // Now interpolate between incident energy bins i and i + 1
-  if (!histogram_interp && n_energy_out > 1 && k >= n_discrete) {
-    // Interpolation for energy E1 and EK
-    n_energy_out = distribution_[i].e_out.size();
-    n_discrete = distribution_[i].n_discrete;
-    const double E_i_1 = distribution_[i].e_out[n_discrete];
-    const double E_i_K = distribution_[i].e_out[n_energy_out - 1];
+    // Now interpolate between incident energy bins i and i + 1
+    if (!histogram_interp && n_energy_out > 1) {
+      // Interpolation for energy E1 and EK
+      n_energy_out = distribution_[i].e_out.size();
+      n_discrete = distribution_[i].n_discrete;
+      const double E_i_1 = distribution_[i].e_out[n_discrete];
+      const double E_i_K = distribution_[i].e_out[n_energy_out - 1];
 
-    n_energy_out = distribution_[i + 1].e_out.size();
-    n_discrete = distribution_[i + 1].n_discrete;
-    const double E_i1_1 = distribution_[i + 1].e_out[n_discrete];
-    const double E_i1_K = distribution_[i + 1].e_out[n_energy_out - 1];
+      n_energy_out = distribution_[i + 1].e_out.size();
+      n_discrete = distribution_[i + 1].n_discrete;
+      const double E_i1_1 = distribution_[i + 1].e_out[n_discrete];
+      const double E_i1_K = distribution_[i + 1].e_out[n_energy_out - 1];
 
-    const double E_1 = E_i_1 + r * (E_i1_1 - E_i_1);
-    const double E_K = E_i_K + r * (E_i1_K - E_i_K);
+      const double E_1 = E_i_1 + r * (E_i1_1 - E_i_1);
+      const double E_K = E_i_K + r * (E_i1_K - E_i_K);
 
-    if (l == i) {
-      return E_1 + (E_out - E_i_1) * (E_K - E_1) / (E_i_K - E_i_1);
+      if (l == i) {
+        return E_1 + (E_out - E_i_1) * (E_K - E_1) / (E_i_K - E_i_1);
+      } else {
+        return E_1 + (E_out - E_i1_1) * (E_K - E_1) / (E_i1_K - E_i1_1);
+      }
     } else {
-      return E_1 + (E_out - E_i1_1) * (E_K - E_1) / (E_i1_K - E_i1_1);
+      return E_out;
     }
-  } else {
-    return E_out;
   }
 }
 
diff --git a/src/secondary_correlated.cpp b/src/secondary_correlated.cpp
index 42ef67ced42..642ecd82b78 100644
--- a/src/secondary_correlated.cpp
+++ b/src/secondary_correlated.cpp
@@ -210,34 +210,37 @@ Distribution& CorrelatedAngleEnergy::sample_dist(
   }
 
   double E_l_k = distribution_[l].e_out[k];
-  double p_l_k = distribution_[l].p[k];
-  if (distribution_[l].interpolation == Interpolation::histogram) {
-    // Histogram interpolation
-    if (p_l_k > 0.0 && k >= n_discrete) {
-      E_out = E_l_k + (r1 - c_k) / p_l_k;
-    } else {
-      E_out = E_l_k;
-    }
-
-  } else if (distribution_[l].interpolation == Interpolation::lin_lin) {
-    // Linear-linear interpolation
-    double E_l_k1 = distribution_[l].e_out[k + 1];
-    double p_l_k1 = distribution_[l].p[k + 1];
-
-    double frac = (p_l_k1 - p_l_k) / (E_l_k1 - E_l_k);
-    if (frac == 0.0) {
-      E_out = E_l_k + (r1 - c_k) / p_l_k;
-    } else {
-      E_out =
-        E_l_k +
-        (std::sqrt(std::max(0.0, p_l_k * p_l_k + 2.0 * frac * (r1 - c_k))) -
-          p_l_k) /
-          frac;
+  if (k < n_discrete) {
+    // Discrete case
+    E_out = E_l_k;
+  } else {
+    // Continuous case
+    double p_l_k = distribution_[l].p[k];
+    if (distribution_[l].interpolation == Interpolation::histogram) {
+      // Histogram interpolation
+      if (p_l_k > 0.0) {
+        E_out = E_l_k + (r1 - c_k) / p_l_k;
+      } else {
+        E_out = E_l_k;
+      }
+    } else if (distribution_[l].interpolation == Interpolation::lin_lin) {
+      // Linear-linear interpolation
+      double E_l_k1 = distribution_[l].e_out[k + 1];
+      double p_l_k1 = distribution_[l].p[k + 1];
+
+      double frac = (p_l_k1 - p_l_k) / (E_l_k1 - E_l_k);
+      if (frac == 0.0) {
+        E_out = E_l_k + (r1 - c_k) / p_l_k;
+      } else {
+        E_out =
+          E_l_k +
+          (std::sqrt(std::max(0.0, p_l_k * p_l_k + 2.0 * frac * (r1 - c_k))) -
+            p_l_k) /
+            frac;
+      }
     }
-  }
 
-  // Now interpolate between incident energy bins i and i + 1
-  if (k >= n_discrete) {
+    // Now interpolate between incident energy bins i and i + 1
     if (l == i) {
       E_out = E_1 + (E_out - E_i_1) * (E_K - E_1) / (E_i_K - E_i_1);
     } else {
diff --git a/src/secondary_kalbach.cpp b/src/secondary_kalbach.cpp
index 018ce1c8a9c..a5d9b9ad0b4 100644
--- a/src/secondary_kalbach.cpp
+++ b/src/secondary_kalbach.cpp
@@ -169,46 +169,53 @@ void KalbachMann::sample_params(
   }
 
   double E_l_k = distribution_[l].e_out[k];
-  double p_l_k = distribution_[l].p[k];
-  if (distribution_[l].interpolation == Interpolation::histogram) {
-    // Histogram interpolation
-    if (p_l_k > 0.0 && k >= n_discrete) {
-      E_out = E_l_k + (r1 - c_k) / p_l_k;
-    } else {
-      E_out = E_l_k;
-    }
-
-    // Determine Kalbach-Mann parameters
+  if (k < n_discrete) {
+    // Discrete case
+    E_out = E_l_k;
     km_r = distribution_[l].r[k];
     km_a = distribution_[l].a[k];
 
   } else {
-    // Linear-linear interpolation
-    double E_l_k1 = distribution_[l].e_out[k + 1];
-    double p_l_k1 = distribution_[l].p[k + 1];
+    // Continuous case
+    double p_l_k = distribution_[l].p[k];
+    if (distribution_[l].interpolation == Interpolation::histogram) {
+      // Histogram interpolation
+      if (p_l_k > 0.0) {
+        E_out = E_l_k + (r1 - c_k) / p_l_k;
+      } else {
+        E_out = E_l_k;
+      }
+
+      // Determine Kalbach-Mann parameters
+      km_r = distribution_[l].r[k];
+      km_a = distribution_[l].a[k];
 
-    double frac = (p_l_k1 - p_l_k) / (E_l_k1 - E_l_k);
-    if (frac == 0.0) {
-      E_out = E_l_k + (r1 - c_k) / p_l_k;
     } else {
-      E_out =
-        E_l_k +
-        (std::sqrt(std::max(0.0, p_l_k * p_l_k + 2.0 * frac * (r1 - c_k))) -
-          p_l_k) /
-          frac;
-    }
+      // Linear-linear interpolation
+      double E_l_k1 = distribution_[l].e_out[k + 1];
+      double p_l_k1 = distribution_[l].p[k + 1];
+
+      double frac = (p_l_k1 - p_l_k) / (E_l_k1 - E_l_k);
+      if (frac == 0.0) {
+        E_out = E_l_k + (r1 - c_k) / p_l_k;
+      } else {
+        E_out =
+          E_l_k +
+          (std::sqrt(std::max(0.0, p_l_k * p_l_k + 2.0 * frac * (r1 - c_k))) -
+            p_l_k) /
+            frac;
+      }
 
-    // Determine Kalbach-Mann parameters
-    km_r = distribution_[l].r[k] +
-           (E_out - E_l_k) / (E_l_k1 - E_l_k) *
-             (distribution_[l].r[k + 1] - distribution_[l].r[k]);
-    km_a = distribution_[l].a[k] +
-           (E_out - E_l_k) / (E_l_k1 - E_l_k) *
-             (distribution_[l].a[k + 1] - distribution_[l].a[k]);
-  }
+      // Determine Kalbach-Mann parameters
+      km_r = distribution_[l].r[k] +
+             (E_out - E_l_k) / (E_l_k1 - E_l_k) *
+               (distribution_[l].r[k + 1] - distribution_[l].r[k]);
+      km_a = distribution_[l].a[k] +
+             (E_out - E_l_k) / (E_l_k1 - E_l_k) *
+               (distribution_[l].a[k + 1] - distribution_[l].a[k]);
+    }
 
-  // Now interpolate between incident energy bins i and i + 1
-  if (k >= n_discrete) {
+    // Now interpolate between incident energy bins i and i + 1
     if (l == i) {
       E_out = E_1 + (E_out - E_i_1) * (E_K - E_1) / (E_i_K - E_i_1);
     } else {