From de50a7bfcac3378b1c4ac233f14d8aa9d1bcd10b Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Wed, 19 Jan 2022 09:14:41 -0600
Subject: [PATCH] Generalize russian_roulette function

---
 include/openmc/physics_common.h |  4 +++-
 src/physics.cpp                 |  6 +++---
 src/physics_common.cpp          | 14 ++++++--------
 src/physics_mg.cpp              |  6 +++---
 src/weight_windows.cpp          | 11 +++--------
 5 files changed, 18 insertions(+), 23 deletions(-)

diff --git a/include/openmc/physics_common.h b/include/openmc/physics_common.h
index 1398557dc..2c5e79dad 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_cutoff Weight below which particles are rouletted
+void russian_roulette(Particle& p, double weight_survive);
 
 } // namespace openmc
 #endif // OPENMC_PHYSICS_COMMON_H
diff --git a/src/physics.cpp b/src/physics.cpp
index 345b79947..01af1259f 100644
--- a/src/physics.cpp
+++ b/src/physics.cpp
@@ -151,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);
+    }
   }
 }
 
diff --git a/src/physics_common.cpp b/src/physics_common.cpp
index 964678fbe..7e65c9bbf 100644
--- a/src/physics_common.cpp
+++ b/src/physics_common.cpp
@@ -9,15 +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 (settings::weight_survive * prn(p.current_seed()) < p.wgt()) {
-      p.wgt() = settings::weight_survive;
-    } else {
-      p.wgt() = 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 d363c2d8e..6e89fde5e 100644
--- a/src/physics_mg.cpp
+++ b/src/physics_mg.cpp
@@ -64,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);
+    }
   }
 }
 
diff --git a/src/weight_windows.cpp b/src/weight_windows.cpp
index e742d48c3..e753fc6ac 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()