diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h
index 38c0aa1dfee..a1804f7be2f 100644
--- a/include/openmc/particle_data.h
+++ b/include/openmc/particle_data.h
@@ -244,7 +244,6 @@ class ParticleData {
   Position r_last_;         //!< previous coordinates
   Direction u_last_;        //!< previous direction coordinates
   double wgt_last_ {1.0};   //!< pre-collision particle weight
-  double wgt_absorb_ {0.0}; //!< weight absorbed for survival biasing
 
   // What event took place
   bool fission_ {false};  //!< did particle cause implicit fission
@@ -375,8 +374,6 @@ class ParticleData {
   const Position& u_last() const { return u_last_; }
   double& wgt_last() { return wgt_last_; }
   const double& wgt_last() const { return wgt_last_; }
-  double& wgt_absorb() { return wgt_absorb_; }
-  const double& wgt_absorb() const { return wgt_absorb_; }
 
   bool& fission() { return fission_; }
   TallyEvent& event() { return event_; }
diff --git a/include/openmc/physics_common.h b/include/openmc/physics_common.h
index 1398557dc58..e38a3c7f883 100644
--- a/include/openmc/physics_common.h
+++ b/include/openmc/physics_common.h
@@ -9,7 +9,9 @@
 namespace openmc {
 
 //! \brief Performs the russian roulette operation for a particle
-void russian_roulette(Particle& p);
+//! \param[in,out] p  Particle object
+//! \param[in] weight_survive Weight assigned to particles that survive
+void russian_roulette(Particle& p, double weight_survive);
 
 } // namespace openmc
 #endif // OPENMC_PHYSICS_COMMON_H
diff --git a/openmc/surface.py b/openmc/surface.py
index 146e16c2861..f2b496f72d7 100644
--- a/openmc/surface.py
+++ b/openmc/surface.py
@@ -2633,7 +2633,7 @@ def rotate(self, rotation, pivot=(0., 0., 0.), order='xyz', inplace=False,
             memo = {}
 
         # If rotated surface not in memo, add it
-        key = (self.surface, tuple(rotation), tuple(pivot), order, inplace)
+        key = (self.surface, tuple(np.ravel(rotation)), tuple(pivot), order, inplace)
         if key not in memo:
             memo[key] = self.surface.rotate(rotation, pivot=pivot, order=order,
                                             inplace=inplace)
diff --git a/src/physics.cpp b/src/physics.cpp
index 196c7dc3b8c..63d764b6779 100644
--- a/src/physics.cpp
+++ b/src/physics.cpp
@@ -134,8 +134,6 @@ void sample_neutron_reaction(Particle& p)
 
   if (p.neutron_xs(i_nuclide).absorption > 0.0) {
     absorption(p, i_nuclide);
-  } else {
-    p.wgt_absorb() = 0.0;
   }
   if (!p.alive())
     return;
@@ -153,9 +151,9 @@ void sample_neutron_reaction(Particle& p)
 
   // Play russian roulette if survival biasing is turned on
   if (settings::survival_biasing) {
-    russian_roulette(p);
-    if (!p.alive())
-      return;
+    if (p.wgt() < settings::weight_cutoff) {
+      russian_roulette(p, settings::weight_survive);
+    }
   }
 }
 
@@ -626,16 +624,15 @@ void absorption(Particle& p, int i_nuclide)
 {
   if (settings::survival_biasing) {
     // Determine weight absorbed in survival biasing
-    p.wgt_absorb() = p.wgt() * p.neutron_xs(i_nuclide).absorption /
-                     p.neutron_xs(i_nuclide).total;
+    const double wgt_absorb = p.wgt() * p.neutron_xs(i_nuclide).absorption /
+                              p.neutron_xs(i_nuclide).total;
 
     // Adjust weight of particle by probability of absorption
-    p.wgt() -= p.wgt_absorb();
-    p.wgt_last() = p.wgt();
+    p.wgt() -= wgt_absorb;
 
     // Score implicit absorption estimate of keff
     if (settings::run_mode == RunMode::EIGENVALUE) {
-      p.keff_tally_absorption() += p.wgt_absorb() *
+      p.keff_tally_absorption() += wgt_absorb *
                                    p.neutron_xs(i_nuclide).nu_fission /
                                    p.neutron_xs(i_nuclide).absorption;
     }
@@ -952,9 +949,8 @@ Direction sample_target_velocity(const Nuclide& nuc, double E, Direction u,
 
       // cdf value at upper bound attainable energy
       double m = (nuc.xs_cdf_[i_E_up + 1] - nuc.xs_cdf_[i_E_up]) /
-          (nuc.energy_0K_[i_E_up + 1] - nuc.energy_0K_[i_E_up]);
-      double cdf_up =
-        nuc.xs_cdf_[i_E_up] + m * (E_up - nuc.energy_0K_[i_E_up]);
+                 (nuc.energy_0K_[i_E_up + 1] - nuc.energy_0K_[i_E_up]);
+      double cdf_up = nuc.xs_cdf_[i_E_up] + m * (E_up - nuc.energy_0K_[i_E_up]);
 
       while (true) {
         // directly sample Maxwellian
@@ -962,9 +958,8 @@ Direction sample_target_velocity(const Nuclide& nuc, double E, Direction u,
 
         // sample a relative energy using the xs cdf
         double cdf_rel = cdf_low + prn(seed) * (cdf_up - cdf_low);
-        int i_E_rel = lower_bound_index(nuc.xs_cdf_.begin() + i_E_low, 
-                                        nuc.xs_cdf_.begin() + i_E_up+2, 
-                                        cdf_rel);
+        int i_E_rel = lower_bound_index(nuc.xs_cdf_.begin() + i_E_low,
+          nuc.xs_cdf_.begin() + i_E_up + 2, cdf_rel);
         double E_rel = nuc.energy_0K_[i_E_low + i_E_rel];
         double m = (nuc.xs_cdf_[i_E_low + i_E_rel + 1] -
                      nuc.xs_cdf_[i_E_low + i_E_rel]) /
diff --git a/src/physics_common.cpp b/src/physics_common.cpp
index f6820f81fd9..7e65c9bbfc4 100644
--- a/src/physics_common.cpp
+++ b/src/physics_common.cpp
@@ -9,17 +9,13 @@ namespace openmc {
 // RUSSIAN_ROULETTE
 //==============================================================================
 
-void russian_roulette(Particle& p)
+void russian_roulette(Particle& p, double weight_survive)
 {
-  if (p.wgt() < settings::weight_cutoff) {
-    if (prn(p.current_seed()) < p.wgt() / settings::weight_survive) {
-      p.wgt() = settings::weight_survive;
-      p.wgt_last() = p.wgt();
-    } else {
-      p.wgt() = 0.;
-      p.wgt_last() = 0.;
-      p.alive() = false;
-    }
+  if (weight_survive * prn(p.current_seed()) < p.wgt()) {
+    p.wgt() = weight_survive;
+  } else {
+    p.wgt() = 0.;
+    p.alive() = false;
   }
 }
 
diff --git a/src/physics_mg.cpp b/src/physics_mg.cpp
index eaf13b4413a..6e89fde5eac 100644
--- a/src/physics_mg.cpp
+++ b/src/physics_mg.cpp
@@ -55,8 +55,6 @@ void sample_reaction(Particle& p)
   // weight of the particle. Otherwise, it checks to see if absorption occurs.
   if (p.macro_xs().absorption > 0.) {
     absorption(p);
-  } else {
-    p.wgt_absorb() = 0.;
   }
   if (!p.alive())
     return;
@@ -66,9 +64,9 @@ void sample_reaction(Particle& p)
 
   // Play Russian roulette if survival biasing is turned on
   if (settings::survival_biasing) {
-    russian_roulette(p);
-    if (!p.alive())
-      return;
+    if (p.wgt() < settings::weight_cutoff) {
+      russian_roulette(p, settings::weight_survive);
+    }
   }
 }
 
@@ -216,15 +214,14 @@ void absorption(Particle& p)
 {
   if (settings::survival_biasing) {
     // Determine weight absorbed in survival biasing
-    p.wgt_absorb() = p.wgt() * p.macro_xs().absorption / p.macro_xs().total;
+    double wgt_absorb = p.wgt() * p.macro_xs().absorption / p.macro_xs().total;
 
     // Adjust weight of particle by the probability of absorption
-    p.wgt() -= p.wgt_absorb();
-    p.wgt_last() = p.wgt();
+    p.wgt() -= wgt_absorb;
 
     // Score implicit absorpion estimate of keff
     p.keff_tally_absorption() +=
-      p.wgt_absorb() * p.macro_xs().nu_fission / p.macro_xs().absorption;
+      wgt_absorb * p.macro_xs().nu_fission / p.macro_xs().absorption;
   } else {
     if (p.macro_xs().absorption > prn(p.current_seed()) * p.macro_xs().total) {
       p.keff_tally_absorption() +=
diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp
index fd75091661b..bd516ca18ab 100644
--- a/src/tallies/tally_scoring.cpp
+++ b/src/tallies/tally_scoring.cpp
@@ -203,10 +203,10 @@ double score_fission_q(const Particle& p, int score_bin, const Tally& tally,
       // No fission events occur if survival biasing is on -- need to
       // calculate fraction of absorptions that would have resulted in
       // fission scaled by the Q-value
-      if (p.neutron_xs(p.event_nuclide()).absorption > 0) {
-        return p.wgt_absorb() * get_nuc_fission_q(nuc, p, score_bin) *
+      if (p.neutron_xs(p.event_nuclide()).total > 0) {
+        return p.wgt_last() * get_nuc_fission_q(nuc, p, score_bin) *
                p.neutron_xs(p.event_nuclide()).fission * flux /
-               p.neutron_xs(p.event_nuclide()).absorption;
+               p.neutron_xs(p.event_nuclide()).total;
       }
     } else {
       // Skip any non-absorption events
@@ -291,7 +291,6 @@ double get_nuclide_neutron_heating(
 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) {
@@ -318,17 +317,12 @@ double score_neutron_heating(const Particle& p, const Tally& tally, double flux,
       }
     }
   }
-  score = heating_xs * flux;
+  double score = heating_xs * flux;
   if (tally.estimator_ == TallyEstimator::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();
-    }
+    score *= p.wgt_last();
   }
   return score;
 }
@@ -552,6 +546,14 @@ void score_general_ce(Particle& p, int i_tally, int start_index,
   // Get the pre-collision energy of the particle.
   auto E = p.E_last();
 
+  // Determine how much weight was absorbed due to survival biasing
+  double wgt_absorb = 0.0;
+  if (tally.estimator_ == TallyEstimator::ANALOG &&
+      settings::survival_biasing) {
+    wgt_absorb = p.wgt_last() * p.neutron_xs(p.event_nuclide()).absorption /
+                 p.neutron_xs(p.event_nuclide()).total;
+  }
+
   using Type = ParticleType;
 
   for (auto i = 0; i < tally.scores_.size(); ++i) {
@@ -565,16 +567,8 @@ void score_general_ce(Particle& p, int i_tally, int start_index,
         // All events score to a flux bin. We actually use a collision estimator
         // in place of an analog one since there is no way to count 'events'
         // exactly for the flux
-        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 (p.type() == Type::neutron || p.type() == Type::photon) {
-          score *= flux / p.macro_xs().total;
+          score = flux * p.wgt_last() / p.macro_xs().total;
         } else {
           score = 0.;
         }
@@ -587,13 +581,7 @@ void score_general_ce(Particle& p, int i_tally, int start_index,
       if (tally.estimator_ == TallyEstimator::ANALOG) {
         // All events will score to the total reaction rate. We can just use
         // use the weight of the particle entering the collision as the score
-        if (settings::survival_biasing) {
-          // We need to account for the fact that some weight was already
-          // absorbed
-          score = (p.wgt_last() + p.wgt_absorb()) * flux;
-        } else {
-          score = p.wgt_last() * flux;
-        }
+        score = p.wgt_last() * flux;
 
       } else {
         if (i_nuclide >= 0) {
@@ -616,14 +604,7 @@ void score_general_ce(Particle& p, int i_tally, int start_index,
         // All events score to an inverse velocity bin. We actually use a
         // collision estimator in place of an analog one since there is no way
         // to count 'events' exactly for the inverse velocity
-        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();
-        }
-        score *= flux / p.macro_xs().total;
+        score = flux * p.wgt_last() / p.macro_xs().total;
       } else {
         score = flux;
       }
@@ -641,7 +622,7 @@ void score_general_ce(Particle& p, int i_tally, int start_index,
           continue;
         // Since only scattering events make it here, again we can use the
         // weight entering the collision as the estimator for the reaction rate
-        score = p.wgt_last() * flux;
+        score = (p.wgt_last() - wgt_absorb) * flux;
       } else {
         if (i_nuclide >= 0) {
           if (p.type() == Type::neutron) {
@@ -672,17 +653,17 @@ void score_general_ce(Particle& p, int i_tally, int start_index,
       // For scattering production, we need to use the pre-collision weight
       // times the yield as the estimate for the number of neutrons exiting a
       // reaction with neutrons in the exit channel
-      if (p.event_mt() == ELASTIC || p.event_mt() == N_LEVEL ||
-          (p.event_mt() >= N_N1 && p.event_mt() <= N_NC)) {
-        // Don't waste time on very common reactions we know have
-        // multiplicities of one.
-        score = p.wgt_last() * flux;
-      } else {
+      score = (p.wgt_last() - wgt_absorb) * flux;
+
+      // Don't waste time on very common reactions we know have multiplicities
+      // of one.
+      if (p.event_mt() != ELASTIC && p.event_mt() != N_LEVEL &&
+          !(p.event_mt() >= N_N1 && p.event_mt() <= N_NC)) {
         // Get yield and apply to score
         auto m =
           data::nuclides[p.event_nuclide()]->reaction_index_[p.event_mt()];
         const auto& rxn {*data::nuclides[p.event_nuclide()]->reactions_[m]};
-        score = p.wgt_last() * flux * (*rxn.products_[0].yield_)(E);
+        score *= (*rxn.products_[0].yield_)(E);
       }
       break;
 
@@ -694,7 +675,7 @@ void score_general_ce(Particle& p, int i_tally, int start_index,
         if (settings::survival_biasing) {
           // No absorption events actually occur if survival biasing is on --
           // just use weight absorbed in survival biasing
-          score = p.wgt_absorb() * flux;
+          score = wgt_absorb * flux;
         } else {
           // Skip any event where the particle wasn't absorbed
           if (p.event() == TallyEvent::SCATTER)
@@ -725,16 +706,15 @@ void score_general_ce(Particle& p, int i_tally, int start_index,
       break;
 
     case SCORE_FISSION:
-      if (p.macro_xs().absorption == 0)
+      if (p.macro_xs().fission == 0)
         continue;
       if (tally.estimator_ == TallyEstimator::ANALOG) {
         if (settings::survival_biasing) {
-          // No fission events occur if survival biasing is on -- need to
-          // calculate fraction of absorptions that would have resulted in
-          // fission
-          if (p.neutron_xs(p.event_nuclide()).absorption > 0) {
-            score = p.wgt_absorb() * p.neutron_xs(p.event_nuclide()).fission /
-                    p.neutron_xs(p.event_nuclide()).absorption * flux;
+          // No fission events occur if survival biasing is on -- use collision
+          // estimator instead
+          if (p.neutron_xs(p.event_nuclide()).total > 0) {
+            score = p.wgt_last() * p.neutron_xs(p.event_nuclide()).fission /
+                    p.neutron_xs(p.event_nuclide()).total * flux;
           } else {
             score = 0.;
           }
@@ -758,7 +738,7 @@ void score_general_ce(Particle& p, int i_tally, int start_index,
       break;
 
     case SCORE_NU_FISSION:
-      if (p.macro_xs().absorption == 0)
+      if (p.macro_xs().fission == 0)
         continue;
       if (tally.estimator_ == TallyEstimator::ANALOG) {
         if (settings::survival_biasing || p.fission()) {
@@ -770,13 +750,11 @@ void score_general_ce(Particle& p, int i_tally, int start_index,
           }
         }
         if (settings::survival_biasing) {
-          // No fission events occur if survival biasing is on -- need to
-          // calculate fraction of absorptions that would have resulted in
-          // nu-fission
-          if (p.neutron_xs(p.event_nuclide()).absorption > 0) {
-            score = p.wgt_absorb() *
-                    p.neutron_xs(p.event_nuclide()).nu_fission /
-                    p.neutron_xs(p.event_nuclide()).absorption * flux;
+          // No fission events occur if survival biasing is on -- use collision
+          // estimator instead
+          if (p.neutron_xs(p.event_nuclide()).total > 0) {
+            score = p.wgt_last() * p.neutron_xs(p.event_nuclide()).nu_fission /
+                    p.neutron_xs(p.event_nuclide()).total * flux;
           } else {
             score = 0.;
           }
@@ -801,7 +779,7 @@ void score_general_ce(Particle& p, int i_tally, int start_index,
       break;
 
     case SCORE_PROMPT_NU_FISSION:
-      if (p.macro_xs().absorption == 0)
+      if (p.macro_xs().fission == 0)
         continue;
       if (tally.estimator_ == TallyEstimator::ANALOG) {
         if (settings::survival_biasing || p.fission()) {
@@ -816,11 +794,11 @@ void score_general_ce(Particle& p, int i_tally, int start_index,
           // No fission events occur if survival biasing is on -- need to
           // calculate fraction of absorptions that would have resulted in
           // prompt-nu-fission
-          if (p.neutron_xs(p.event_nuclide()).absorption > 0) {
-            score = p.wgt_absorb() * p.neutron_xs(p.event_nuclide()).fission *
+          if (p.neutron_xs(p.event_nuclide()).total > 0) {
+            score = p.wgt_last() * p.neutron_xs(p.event_nuclide()).fission *
                     data::nuclides[p.event_nuclide()]->nu(
                       E, ReactionProduct::EmissionMode::prompt) /
-                    p.neutron_xs(p.event_nuclide()).absorption * flux;
+                    p.neutron_xs(p.event_nuclide()).total * flux;
           } else {
             score = 0.;
           }
@@ -863,7 +841,7 @@ void score_general_ce(Particle& p, int i_tally, int start_index,
       break;
 
     case SCORE_DELAYED_NU_FISSION:
-      if (p.macro_xs().absorption == 0)
+      if (p.macro_xs().fission == 0)
         continue;
       if (tally.estimator_ == TallyEstimator::ANALOG) {
         if (settings::survival_biasing || p.fission()) {
@@ -878,7 +856,7 @@ void score_general_ce(Particle& p, int i_tally, int start_index,
           // No fission events occur if survival biasing is on -- need to
           // calculate fraction of absorptions that would have resulted in
           // delayed-nu-fission
-          if (p.neutron_xs(p.event_nuclide()).absorption > 0 &&
+          if (p.neutron_xs(p.event_nuclide()).total > 0 &&
               data::nuclides[p.event_nuclide()]->fissionable_) {
             if (tally.delayedgroup_filter_ != C_NONE) {
               auto i_dg_filt = tally.filters()[tally.delayedgroup_filter_];
@@ -890,9 +868,9 @@ void score_general_ce(Particle& p, int i_tally, int start_index,
                 auto dg = filt.groups()[d_bin];
                 auto yield = data::nuclides[p.event_nuclide()]->nu(
                   E, ReactionProduct::EmissionMode::delayed, dg);
-                score = p.wgt_absorb() * yield *
+                score = p.wgt_last() * yield *
                         p.neutron_xs(p.event_nuclide()).fission /
-                        p.neutron_xs(p.event_nuclide()).absorption * flux;
+                        p.neutron_xs(p.event_nuclide()).total * flux;
                 score_fission_delayed_dg(
                   i_tally, d_bin, score, score_index, p.filter_matches());
               }
@@ -901,10 +879,10 @@ void score_general_ce(Particle& p, int i_tally, int start_index,
               // If the delayed group filter is not present, compute the score
               // by multiplying the absorbed weight by the fraction of the
               // delayed-nu-fission xs to the absorption xs
-              score = p.wgt_absorb() * p.neutron_xs(p.event_nuclide()).fission *
+              score = p.wgt_last() * p.neutron_xs(p.event_nuclide()).fission *
                       data::nuclides[p.event_nuclide()]->nu(
                         E, ReactionProduct::EmissionMode::delayed) /
-                      p.neutron_xs(p.event_nuclide()).absorption * flux;
+                      p.neutron_xs(p.event_nuclide()).total * flux;
             }
           }
         } else {
@@ -1007,7 +985,7 @@ void score_general_ce(Particle& p, int i_tally, int start_index,
       break;
 
     case SCORE_DECAY_RATE:
-      if (p.macro_xs().absorption == 0)
+      if (p.macro_xs().fission == 0)
         continue;
       if (tally.estimator_ == TallyEstimator::ANALOG) {
         if (settings::survival_biasing) {
@@ -1015,8 +993,7 @@ void score_general_ce(Particle& p, int i_tally, int start_index,
           // calculate fraction of absorptions that would have resulted in
           // delayed-nu-fission
           const auto& nuc {*data::nuclides[p.event_nuclide()]};
-          if (p.neutron_xs(p.event_nuclide()).absorption > 0 &&
-              nuc.fissionable_) {
+          if (p.neutron_xs(p.event_nuclide()).total > 0 && nuc.fissionable_) {
             const auto& rxn {*nuc.fission_rx_[0]};
             if (tally.delayedgroup_filter_ != C_NONE) {
               auto i_dg_filt = tally.filters()[tally.delayedgroup_filter_];
@@ -1029,10 +1006,9 @@ void score_general_ce(Particle& p, int i_tally, int start_index,
                 auto yield =
                   nuc.nu(E, ReactionProduct::EmissionMode::delayed, d);
                 auto rate = rxn.products_[d].decay_rate_;
-                score = p.wgt_absorb() * yield *
+                score = p.wgt_last() * yield *
                         p.neutron_xs(p.event_nuclide()).fission /
-                        p.neutron_xs(p.event_nuclide()).absorption * rate *
-                        flux;
+                        p.neutron_xs(p.event_nuclide()).total * rate * flux;
                 score_fission_delayed_dg(
                   i_tally, d_bin, score, score_index, p.filter_matches());
               }
@@ -1048,13 +1024,17 @@ void score_general_ce(Particle& p, int i_tally, int start_index,
               // rxn.products_ array to be exceeded. Hence, we use the size
               // of this array and not the MAX_DELAYED_GROUPS constant for
               // this loop.
-              for (auto d = 0; d < rxn.products_.size() - 2; ++d) {
+              for (auto d = 1; d < rxn.products_.size(); ++d) {
+                const auto& product = rxn.products_[d];
+                if (product.particle_ != Type::neutron)
+                  continue;
+
                 auto yield =
-                  nuc.nu(E, ReactionProduct::EmissionMode::delayed, d + 1);
-                auto rate = rxn.products_[d + 1].decay_rate_;
-                score += rate * p.wgt_absorb() *
+                  nuc.nu(E, ReactionProduct::EmissionMode::delayed, d);
+                auto rate = product.decay_rate_;
+                score += rate * p.wgt_last() *
                          p.neutron_xs(p.event_nuclide()).fission * yield /
-                         p.neutron_xs(p.event_nuclide()).absorption * flux;
+                         p.neutron_xs(p.event_nuclide()).total * flux;
               }
             }
           }
@@ -1123,10 +1103,13 @@ void score_general_ce(Particle& p, int i_tally, int start_index,
             // rxn.products_ array to be exceeded. Hence, we use the size
             // of this array and not the MAX_DELAYED_GROUPS constant for
             // this loop.
-            for (auto d = 0; d < rxn.products_.size() - 2; ++d) {
-              auto yield =
-                nuc.nu(E, ReactionProduct::EmissionMode::delayed, d + 1);
-              auto rate = rxn.products_[d + 1].decay_rate_;
+            for (auto d = 1; d < rxn.products_.size(); ++d) {
+              const auto& product = rxn.products_[d];
+              if (product.particle_ != Type::neutron)
+                continue;
+
+              auto yield = nuc.nu(E, ReactionProduct::EmissionMode::delayed, d);
+              auto rate = product.decay_rate_;
               score += p.neutron_xs(i_nuclide).fission * flux * yield *
                        atom_density * rate;
             }
@@ -1174,10 +1157,14 @@ void score_general_ce(Particle& p, int i_tally, int start_index,
                   // rxn.products_ array to be exceeded. Hence, we use the size
                   // of this array and not the MAX_DELAYED_GROUPS constant for
                   // this loop.
-                  for (auto d = 0; d < rxn.products_.size() - 2; ++d) {
+                  for (auto d = 1; d < rxn.products_.size(); ++d) {
+                    const auto& product = rxn.products_[d];
+                    if (product.particle_ != Type::neutron)
+                      continue;
+
                     auto yield =
-                      nuc.nu(E, ReactionProduct::EmissionMode::delayed, d + 1);
-                    auto rate = rxn.products_[d + 1].decay_rate_;
+                      nuc.nu(E, ReactionProduct::EmissionMode::delayed, d);
+                    auto rate = product.decay_rate_;
                     score += p.neutron_xs(j_nuclide).fission * yield *
                              atom_density * flux * rate;
                   }
@@ -1190,7 +1177,7 @@ void score_general_ce(Particle& p, int i_tally, int start_index,
       break;
 
     case SCORE_KAPPA_FISSION:
-      if (p.macro_xs().absorption == 0.)
+      if (p.macro_xs().fission == 0.)
         continue;
       score = 0.;
       // Kappa-fission values are determined from the Q-value listed for the
@@ -1201,12 +1188,11 @@ void score_general_ce(Particle& p, int i_tally, int start_index,
           // calculate fraction of absorptions that would have resulted in
           // fission scaled by the Q-value
           const auto& nuc {*data::nuclides[p.event_nuclide()]};
-          if (p.neutron_xs(p.event_nuclide()).absorption > 0 &&
-              nuc.fissionable_) {
+          if (p.neutron_xs(p.event_nuclide()).total > 0 && nuc.fissionable_) {
             const auto& rxn {*nuc.fission_rx_[0]};
-            score = p.wgt_absorb() * rxn.q_value_ *
+            score = p.wgt_last() * rxn.q_value_ *
                     p.neutron_xs(p.event_nuclide()).fission /
-                    p.neutron_xs(p.event_nuclide()).absorption * flux;
+                    p.neutron_xs(p.event_nuclide()).total * flux;
           }
         } else {
           // Skip any non-absorption events
@@ -1262,7 +1248,7 @@ void score_general_ce(Particle& p, int i_tally, int start_index,
         // Check if event MT matches
         if (p.event_mt() != ELASTIC)
           continue;
-        score = p.wgt_last() * flux;
+        score = (p.wgt_last() - wgt_absorb) * flux;
       } else {
         if (i_nuclide >= 0) {
           if (p.neutron_xs(i_nuclide).elastic == CACHE_INVALID)
@@ -1286,7 +1272,7 @@ void score_general_ce(Particle& p, int i_tally, int start_index,
 
     case SCORE_FISS_Q_PROMPT:
     case SCORE_FISS_Q_RECOV:
-      if (p.macro_xs().absorption == 0.)
+      if (p.macro_xs().fission == 0.)
         continue;
       score =
         score_fission_q(p, score_bin, tally, flux, i_nuclide, atom_density);
@@ -1311,7 +1297,7 @@ void score_general_ce(Particle& p, int i_tally, int start_index,
         // Check if the event MT matches
         if (p.event_mt() != score_bin)
           continue;
-        score = p.wgt_last() * flux;
+        score = (p.wgt_last() - wgt_absorb) * flux;
       } else {
         int m;
         switch (score_bin) {
@@ -1423,12 +1409,11 @@ void score_general_ce(Particle& p, int i_tally, int start_index,
         continue;
 
       if (tally.estimator_ == TallyEstimator::ANALOG) {
-
         // Any other score is assumed to be a MT number. Thus, we just need
         // to check if it matches the MT number of the event
         if (p.event_mt() != score_bin)
           continue;
-        score = p.wgt_last() * flux;
+        score = (p.wgt_last() - wgt_absorb) * flux;
       } else {
         // Any other cross section has to be calculated on-the-fly
         if (score_bin < 2)
@@ -1475,10 +1460,13 @@ void score_general_mg(Particle& p, int i_tally, int start_index,
   // Set the direction and group to use with get_xs
   Direction p_u;
   int p_g;
+  double wgt_absorb = 0.0;
   if (tally.estimator_ == TallyEstimator::ANALOG ||
       tally.estimator_ == TallyEstimator::COLLISION) {
-
     if (settings::survival_biasing) {
+      // Determine weight that was absorbed
+      wgt_absorb = p.wgt_last() * p.neutron_xs(p.event_nuclide()).absorption /
+                   p.neutron_xs(p.event_nuclide()).total;
 
       // Then we either are alive and had a scatter (and so g changed),
       // or are dead and g did not change
@@ -1532,14 +1520,7 @@ void score_general_mg(Particle& p, int i_tally, int start_index,
         // All events score to a flux bin. We actually use a collision estimator
         // in place of an analog one since there is no way to count 'events'
         // exactly for the flux
-        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();
-        }
-        score *= flux / p.macro_xs().total;
+        score = flux * p.wgt_last() / p.macro_xs().total;
       } else {
         score = flux;
       }
@@ -1549,16 +1530,9 @@ void score_general_mg(Particle& p, int i_tally, int start_index,
       if (tally.estimator_ == TallyEstimator::ANALOG) {
         // All events will score to the total reaction rate. We can just use
         // use the weight of the particle entering the collision as the score
-        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();
-        }
-        // TODO: should flux be multiplied in above instead of below?
+        score = flux * p.wgt_last();
         if (i_nuclide >= 0) {
-          score *= flux * atom_density * nuc_xs.get_xs(MgxsType::TOTAL, p_g) /
+          score *= atom_density * nuc_xs.get_xs(MgxsType::TOTAL, p_g) /
                    macro_xs.get_xs(MgxsType::TOTAL, p_g);
         }
       } else {
@@ -1576,18 +1550,12 @@ void score_general_mg(Particle& p, int i_tally, int start_index,
         // All events score to an inverse velocity bin. We actually use a
         // collision estimator in place of an analog one since there is no way
         // to count 'events' exactly for the inverse velocity
-        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();
-        }
+        score = flux * p.wgt_last();
         if (i_nuclide >= 0) {
-          score *= flux * nuc_xs.get_xs(MgxsType::INVERSE_VELOCITY, p_g) /
+          score *= nuc_xs.get_xs(MgxsType::INVERSE_VELOCITY, p_g) /
                    macro_xs.get_xs(MgxsType::TOTAL, p_g);
         } else {
-          score *= flux * macro_xs.get_xs(MgxsType::INVERSE_VELOCITY, p_g) /
+          score *= macro_xs.get_xs(MgxsType::INVERSE_VELOCITY, p_g) /
                    macro_xs.get_xs(MgxsType::TOTAL, p_g);
         }
       } else {
@@ -1606,7 +1574,7 @@ void score_general_mg(Particle& p, int i_tally, int start_index,
           continue;
         // Since only scattering events make it here, again we can use the
         // weight entering the collision as the estimator for the reaction rate
-        score = p.wgt_last() * flux;
+        score = (p.wgt_last() - wgt_absorb) * flux;
         if (i_nuclide >= 0) {
           score *= atom_density *
                    nuc_xs.get_xs(MgxsType::SCATTER_FMU, p.g_last(), &p.g(),
@@ -1634,7 +1602,7 @@ void score_general_mg(Particle& p, int i_tally, int start_index,
         // For scattering production, we need to use the pre-collision weight
         // times the multiplicity as the estimate for the number of neutrons
         // exiting a reaction with neutrons in the exit channel
-        score = p.wgt() * flux;
+        score = (p.wgt_last() - wgt_absorb) * flux;
         // Since we transport based on material data, the angle selected
         // was not selected from the f(mu) for the nuclide.  Therefore
         // adjust the score by the actual probability for that nuclide.
@@ -1660,7 +1628,7 @@ void score_general_mg(Particle& p, int i_tally, int start_index,
         if (settings::survival_biasing) {
           // No absorption events actually occur if survival biasing is on --
           // just use weight absorbed in survival biasing
-          score = p.wgt_absorb() * flux;
+          score = wgt_absorb * flux;
         } else {
           // Skip any event where the particle wasn't absorbed
           if (p.event() == TallyEvent::SCATTER)
@@ -1689,7 +1657,7 @@ void score_general_mg(Particle& p, int i_tally, int start_index,
           // No fission events occur if survival biasing is on -- need to
           // calculate fraction of absorptions that would have resulted in
           // fission
-          score = p.wgt_absorb() * flux;
+          score = wgt_absorb * flux;
         } else {
           // Skip any non-absorption events
           if (p.event() == TallyEvent::SCATTER)
@@ -1729,7 +1697,7 @@ void score_general_mg(Particle& p, int i_tally, int start_index,
           // No fission events occur if survival biasing is on -- need to
           // calculate fraction of absorptions that would have resulted in
           // nu-fission
-          score = p.wgt_absorb() * flux;
+          score = wgt_absorb * flux;
           if (i_nuclide >= 0) {
             score *= atom_density * nuc_xs.get_xs(MgxsType::NU_FISSION, p_g) /
                      macro_xs.get_xs(MgxsType::ABSORPTION, p_g);
@@ -1776,7 +1744,7 @@ void score_general_mg(Particle& p, int i_tally, int start_index,
           // No fission events occur if survival biasing is on -- need to
           // calculate fraction of absorptions that would have resulted in
           // prompt-nu-fission
-          score = p.wgt_absorb() * flux;
+          score = wgt_absorb * flux;
           if (i_nuclide >= 0) {
             score *= atom_density *
                      nuc_xs.get_xs(MgxsType::PROMPT_NU_FISSION, p_g) /
@@ -1837,7 +1805,7 @@ void score_general_mg(Particle& p, int i_tally, int start_index,
               // Tally each delayed group bin individually
               for (auto d_bin = 0; d_bin < filt.n_bins(); ++d_bin) {
                 auto d = filt.groups()[d_bin] - 1;
-                score = p.wgt_absorb() * flux;
+                score = wgt_absorb * flux;
                 if (i_nuclide >= 0) {
                   score *= nuc_xs.get_xs(MgxsType::DELAYED_NU_FISSION, p_g,
                              nullptr, nullptr, &d) /
@@ -1855,7 +1823,7 @@ void score_general_mg(Particle& p, int i_tally, int start_index,
               // If the delayed group filter is not present, compute the score
               // by multiplying the absorbed weight by the fraction of the
               // delayed-nu-fission xs to the absorption xs
-              score = p.wgt_absorb() * flux;
+              score = wgt_absorb * flux;
               if (i_nuclide >= 0) {
                 score *=
                   nuc_xs.get_xs(MgxsType::DELAYED_NU_FISSION, p_g) / abs_xs;
@@ -1952,7 +1920,7 @@ void score_general_mg(Particle& p, int i_tally, int start_index,
               // Tally each delayed group bin individually
               for (auto d_bin = 0; d_bin < filt.n_bins(); ++d_bin) {
                 auto d = filt.groups()[d_bin] - 1;
-                score = p.wgt_absorb() * flux;
+                score = wgt_absorb * flux;
                 if (i_nuclide >= 0) {
                   score *= nuc_xs.get_xs(
                              MgxsType::DECAY_RATE, p_g, nullptr, nullptr, &d) *
@@ -1978,14 +1946,14 @@ void score_general_mg(Particle& p, int i_tally, int start_index,
               score = 0.;
               for (auto d = 0; d < data::mg.num_delayed_groups_; ++d) {
                 if (i_nuclide >= 0) {
-                  score += p.wgt_absorb() * flux *
+                  score += wgt_absorb * flux *
                            nuc_xs.get_xs(
                              MgxsType::DECAY_RATE, p_g, nullptr, nullptr, &d) *
                            nuc_xs.get_xs(MgxsType::DELAYED_NU_FISSION, p_g,
                              nullptr, nullptr, &d) /
                            abs_xs;
                 } else {
-                  score += p.wgt_absorb() * flux *
+                  score += wgt_absorb * flux *
                            macro_xs.get_xs(
                              MgxsType::DECAY_RATE, p_g, nullptr, nullptr, &d) *
                            macro_xs.get_xs(MgxsType::DELAYED_NU_FISSION, p_g,
@@ -2093,7 +2061,7 @@ void score_general_mg(Particle& p, int i_tally, int start_index,
           // No fission events occur if survival biasing is on -- need to
           // calculate fraction of absorptions that would have resulted in
           // fission scaled by the Q-value
-          score = p.wgt_absorb() * flux;
+          score = wgt_absorb * flux;
         } else {
           // Skip any non-absorption events
           if (p.event() == TallyEvent::SCATTER)
@@ -2362,11 +2330,7 @@ void score_collision_tally(Particle& p)
   // Determine the collision estimate of the flux
   double flux = 0.0;
   if (p.type() == ParticleType::neutron || p.type() == ParticleType::photon) {
-    if (!settings::survival_biasing) {
-      flux = p.wgt_last() / p.macro_xs().total;
-    } else {
-      flux = (p.wgt_last() + p.wgt_absorb()) / p.macro_xs().total;
-    }
+    flux = p.wgt_last() / p.macro_xs().total;
   }
 
   for (auto i_tally : model::active_collision_tallies) {
diff --git a/src/weight_windows.cpp b/src/weight_windows.cpp
index e742d48c3f8..e753fc6acf5 100644
--- a/src/weight_windows.cpp
+++ b/src/weight_windows.cpp
@@ -5,6 +5,7 @@
 #include "openmc/hdf5_interface.h"
 #include "openmc/particle.h"
 #include "openmc/particle_data.h"
+#include "openmc/physics_common.h"
 #include "openmc/search.h"
 #include "openmc/xml_interface.h"
 
@@ -94,14 +95,8 @@ void apply_weight_windows(Particle& p)
     // if the particle weight is below the window, play Russian roulette
     double weight_survive =
       std::min(weight * weight_window.max_split, weight_window.survival_weight);
-    if (weight_survive * prn(p.current_seed()) <= weight) {
-      p.wgt() = weight_survive;
-    } else {
-      p.alive() = false;
-      p.wgt() = 0.0;
-    }
-    // else particle is in the window, continue as normal
-  }
+    russian_roulette(p, weight_survive);
+  } // else particle is in the window, continue as normal
 }
 
 void free_memory_weight_windows()
diff --git a/tests/regression_tests/survival_biasing/results_true.dat b/tests/regression_tests/survival_biasing/results_true.dat
index b9868f3d3ab..edd96834016 100644
--- a/tests/regression_tests/survival_biasing/results_true.dat
+++ b/tests/regression_tests/survival_biasing/results_true.dat
@@ -1,10 +1,10 @@
 k-combined:
 9.826269E-01 1.626276E-02
 tally 1:
-4.212303E+01
-3.550855E+02
-1.760388E+01
-6.205773E+01
+4.209907E+01
+3.546281E+02
+1.759415E+01
+6.197412E+01
 2.165888E+00
 9.392975E-01
 1.873459E+00
@@ -16,5 +16,5 @@ tally 1:
 3.628554E+08
 2.635633E+16
 tally 2:
-1.760388E+01
-6.205773E+01
+1.759415E+01
+6.197412E+01
diff --git a/tests/regression_tests/tallies/results_true.dat b/tests/regression_tests/tallies/results_true.dat
index c73ff33a432..6ff90033120 100644
--- a/tests/regression_tests/tallies/results_true.dat
+++ b/tests/regression_tests/tallies/results_true.dat
@@ -1 +1 @@
-de09325b517aad9f58940e5fcd53b005c57ea008a15699769ab91aba723de141014101e55be2621e34c369beefb6db8c6cf847e241ee26e73e394c0f155138b0
\ No newline at end of file
+683a82a10c4254c7e503f9c86e192070224c0d19c1817e51da7acd74b9529a475020a6cbd658644f7cd18b244cc477d1810521aa720c4fa9354cb6955bce13f6
\ No newline at end of file