diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp
index bbbd3fefb9..1f6a8010a5 100644
--- a/src/tallies/tally_scoring.cpp
+++ b/src/tallies/tally_scoring.cpp
@@ -187,6 +187,11 @@ double get_nuc_fission_q(const Nuclide& nuc, const Particle* p, int score_bin)
   return 0.0;
 }
 
+//! Helper function to score fission energy
+//
+//! Pulled out to support both the fission_q scores and energy deposition
+//! score
+
 double score_fission_q(const Particle* p, int score_bin, const Tally& tally, 
   double flux, int i_nuclide, double atom_density)
 {
@@ -236,6 +241,69 @@ double score_fission_q(const Particle* p, int score_bin, const Tally& tally,
   return 0.0;
 }
 
+//! Helper function to obtain the kerma coefficient for a given nuclide
+
+double get_nuclide_neutron_heating(const Particle* p, const Nuclide& nuc,
+    int rxn_index, int i_nuclide)
+{
+  size_t mt = nuc.reaction_index_[rxn_index];
+  if (mt == C_NONE) return 0.0;
+  auto i_temp = p->neutron_xs_[i_nuclide].index_temp;
+  if (i_temp < 0) return 0.0; // Can be true due to multipole
+  const auto& rxn {*nuc.reactions_[mt]};
+  const auto& xs {rxn.xs_[i_temp]};
+  auto i_grid = p->neutron_xs_[i_nuclide].index_grid;
+  if (i_grid < xs.threshold) return 0.0;
+  auto f = p->neutron_xs_[i_nuclide].interp_factor;
+  return (1.0 - f) * xs.value[i_grid-xs.threshold]
+    + f * xs.value[i_grid-xs.threshold+1];
+}
+
+//! Helper function to obtain neutron heating [eV]
+
+double score_neutron_heating(const Particle* p, const Tally& tally, double flux,
+    int rxn_bin, int i_nuclide, double atom_density)
+{
+  double score;
+  // Get heating macroscopic "cross section"
+  double heating_xs;
+  if (i_nuclide >= 0) {
+    const Nuclide& nuc {*data::nuclides[i_nuclide]};    
+    heating_xs = get_nuclide_neutron_heating(p, nuc, rxn_bin, i_nuclide);
+    if (tally.estimator_ == ESTIMATOR_ANALOG) {
+      heating_xs /= p->neutron_xs_[i_nuclide].total;
+    } else {
+      heating_xs *= atom_density;
+    }
+  } else {
+    if (p->material_ != MATERIAL_VOID) {
+      heating_xs = 0.0;
+      const Material& material {*model::materials[p->material_]};
+      for (auto i = 0; i< material.nuclide_.size(); ++i) {
+        int j_nuclide = material.nuclide_[i];
+        double atom_density {material.atom_density_(i)};
+        const Nuclide& nuc {*data::nuclides[j_nuclide]};
+        heating_xs += atom_density * get_nuclide_neutron_heating(p, nuc, rxn_bin, j_nuclide);
+      }
+      if (tally.estimator_ == ESTIMATOR_ANALOG) {
+        heating_xs /= p->macro_xs_.total;
+      }
+    }
+  }
+  score = heating_xs * flux;
+  if (tally.estimator_ == ESTIMATOR_ANALOG) {
+    // All events score to a heating tally bin. We actually use a
+    // collision estimator in place of an analog one since there is no
+    // reaction-wise heating cross section
+    if (settings::survival_biasing) {
+      // Account for the fact that some weight has been absorbed
+      score *= p->wgt_last_ + p->wgt_absorb_;
+    } else {
+      score *= p->wgt_last_;
+    }
+  }
+  return score;
+}
 
 //! Helper function for nu-fission tallies with energyout filters.
 //
@@ -1145,105 +1213,8 @@ score_general_ce(Particle* p, int i_tally, int start_index,
     case SCORE_HEATING:
       score = 0.;
       if (p->type_ == Particle::Type::neutron) {
-        if (tally.estimator_ == ESTIMATOR_ANALOG) {
-          // All events score to a heating tally bin. We actually use a
-          // collision estimator in place of an analog one since there is no
-          // reaction-wise heating cross section
-          if (settings::survival_biasing) {
-            // We need to account for the fact that some weight was already
-            // absorbed
-            score = p->wgt_last_ + p->wgt_absorb_;
-          } else {
-            score = p->wgt_last_;
-          }
-          if (i_nuclide >= 0) {
-            // Calculate nuclide heating cross section
-            double macro_heating = 0.;
-            const auto& nuc {*data::nuclides[i_nuclide]};
-            auto m = nuc.reaction_index_[NEUTRON_HEATING];
-            if (m == C_NONE) continue;
-            const auto& rxn {*nuc.reactions_[m]};
-            auto i_temp = p->neutron_xs_[i_nuclide].index_temp;
-            if (i_temp >= 0) { // Can be false due to multipole
-              auto i_grid = p->neutron_xs_[i_nuclide].index_grid;
-              auto f = p->neutron_xs_[i_nuclide].interp_factor;
-              const auto& xs {rxn.xs_[i_temp]};
-              if (i_grid >= xs.threshold) {
-                macro_heating = ((1.0 - f) * xs.value[i_grid-xs.threshold]
-                  + f * xs.value[i_grid-xs.threshold+1]);
-              }
-            }
-            score *= macro_heating * flux / p->neutron_xs_[i_nuclide].total;
-          } else {
-            if (p->material_ != MATERIAL_VOID) {
-              // Calculate material heating cross section
-              double macro_heating = 0.;
-              const Material& material {*model::materials[p->material_]};
-              for (auto i = 0; i < material.nuclide_.size(); ++i) {
-                auto j_nuclide = material.nuclide_[i];
-                auto atom_density = material.atom_density_(i);
-                const auto& nuc {*data::nuclides[j_nuclide]};
-                auto m = nuc.reaction_index_[NEUTRON_HEATING];
-                if (m == C_NONE) continue;
-                const auto& rxn {*nuc.reactions_[m]};
-                auto i_temp = p->neutron_xs_[j_nuclide].index_temp;
-                if (i_temp >= 0) { // Can be false due to multipole
-                  auto i_grid = p->neutron_xs_[j_nuclide].index_grid;
-                  auto f = p->neutron_xs_[j_nuclide].interp_factor;
-                  const auto& xs {rxn.xs_[i_temp]};
-                  if (i_grid >= xs.threshold) {
-                    macro_heating += ((1.0 - f) * xs.value[i_grid-xs.threshold]
-                      + f * xs.value[i_grid-xs.threshold+1]) * atom_density;
-                  }
-                }
-              }
-              score *= macro_heating * flux / p->macro_xs_.total;
-            } else {
-              score = 0.;
-            }
-          }
-        } else {
-          // Calculate neutron heating cross section on-the-fly
-          if (i_nuclide >= 0) {
-            const auto& nuc {*data::nuclides[i_nuclide]};
-            auto m = nuc.reaction_index_[NEUTRON_HEATING];
-            if (m == C_NONE) continue;
-            const auto& rxn {*nuc.reactions_[m]};
-            auto i_temp = p->neutron_xs_[i_nuclide].index_temp;
-            if (i_temp >= 0) { // Can be false due to multipole
-              auto i_grid = p->neutron_xs_[i_nuclide].index_grid;
-              auto f = p->neutron_xs_[i_nuclide].interp_factor;
-              const auto& xs {rxn.xs_[i_temp]};
-              if (i_grid >= xs.threshold) {
-                score = ((1.0 - f) * xs.value[i_grid-xs.threshold]
-                  + f * xs.value[i_grid-xs.threshold+1]) * atom_density * flux;
-              }
-            }
-          } else {
-            if (p->material_ != MATERIAL_VOID) {
-              const Material& material {*model::materials[p->material_]};
-              for (auto i = 0; i < material.nuclide_.size(); ++i) {
-                auto j_nuclide = material.nuclide_[i];
-                auto atom_density = material.atom_density_(i);
-                const auto& nuc {*data::nuclides[j_nuclide]};
-                auto m = nuc.reaction_index_[NEUTRON_HEATING];
-                if (m == C_NONE) continue;
-                const auto& rxn {*nuc.reactions_[m]};
-                auto i_temp = p->neutron_xs_[j_nuclide].index_temp;
-                if (i_temp >= 0) { // Can be false due to multipole
-                  auto i_grid = p->neutron_xs_[j_nuclide].index_grid;
-                  auto f = p->neutron_xs_[j_nuclide].interp_factor;
-                  const auto& xs {rxn.xs_[i_temp]};
-                  if (i_grid >= xs.threshold) {
-                    score += ((1.0 - f) * xs.value[i_grid-xs.threshold]
-                      + f * xs.value[i_grid-xs.threshold+1]) * atom_density
-                      * flux;
-                  }
-                }
-              }
-            }
-          }
-        }
+        score = score_neutron_heating(p, tally, flux, NEUTRON_HEATING,
+            i_nuclide, atom_density);
       } else if (p->type_ == Particle::Type::photon) {
         if (tally.estimator_ == ESTIMATOR_ANALOG) {
           // Score direct energy deposition in the collision