From 529067f95841a9410e15d9b2a8dc906ed9f24ac2 Mon Sep 17 00:00:00 2001
From: Nathan Glaser <glaserna@oregonstate.edu>
Date: Wed, 1 Apr 2026 18:44:43 -0700
Subject: [PATCH 01/53] add split method to techniques

---
 mcdc/transport/simulation.py |  3 ++-
 mcdc/transport/technique.py  | 17 +++++++++++++++++
 2 files changed, 19 insertions(+), 1 deletion(-)

diff --git a/mcdc/transport/simulation.py b/mcdc/transport/simulation.py
index ba0be7639..796e32beb 100644
--- a/mcdc/transport/simulation.py
+++ b/mcdc/transport/simulation.py
@@ -305,8 +305,9 @@ def step_particle(particle_container, mcdc, data):
     if particle["event"] & EVENT_TIME_BOUNDARY:
         particle["alive"] = False
 
-    # Weight roulette
+    # Weight splitting / rouletting
     if particle["alive"]:
+        technique.weight_split(particle_container, mcdc)
         technique.weight_roulette(particle_container, mcdc)
 
 
diff --git a/mcdc/transport/technique.py b/mcdc/transport/technique.py
index f8b15c05f..340e73c14 100644
--- a/mcdc/transport/technique.py
+++ b/mcdc/transport/technique.py
@@ -10,6 +10,23 @@
 import mcdc.transport.particle_bank as particle_bank_module
 import mcdc.transport.rng as rng
 
+# ======================================================================================
+# Weight Splitting
+# ======================================================================================
+
+
+@njit
+def weight_split(particle_container, mcdc):
+    particle = particle_container[0]
+    weight = particle["w"]
+    threshold = mcdc["weight_split"]["weight_threshold"]
+    if weight > threshold:
+        num_split = math.ceil(weight / threshold)
+        particle["w"] = weight / num_split
+        for _ in range(num_split - 1):
+            particle_bank_module.bank_active_particle(particle_container, mcdc)
+
+
 # ======================================================================================
 # Weight Roulette
 # ======================================================================================

From 023cbbec0071e45f708622d53a84496a48a62f57 Mon Sep 17 00:00:00 2001
From: Nathan Glaser <glaserna@oregonstate.edu>
Date: Wed, 1 Apr 2026 18:45:34 -0700
Subject: [PATCH 02/53] add split method to simulation

---
 mcdc/object_/simulation.py |  3 +++
 mcdc/object_/technique.py  | 18 ++++++++++++++++++
 2 files changed, 21 insertions(+)

diff --git a/mcdc/object_/simulation.py b/mcdc/object_/simulation.py
index f7a1156e9..3acd888a4 100644
--- a/mcdc/object_/simulation.py
+++ b/mcdc/object_/simulation.py
@@ -4,6 +4,7 @@
 from mcdc.object_.technique import (
     ImplicitCapture,
     PopulationControl,
+    WeightSplit,
     WeightRoulette,
     WeightedEmission,
 )
@@ -77,6 +78,7 @@ class Simulation(ObjectSingleton):
     # Techniques
     implicit_capture: ImplicitCapture
     weighted_emission: WeightedEmission
+    weight_split: WeightSplit
     weight_roulette: WeightRoulette
     population_control: PopulationControl
 
@@ -159,6 +161,7 @@ def __init__(self):
         # Techniques
         self.implicit_capture = ImplicitCapture()
         self.weighted_emission = WeightedEmission()
+        self.weight_split = WeightSplit()
         self.weight_roulette = WeightRoulette()
         self.population_control = PopulationControl()
 
diff --git a/mcdc/object_/technique.py b/mcdc/object_/technique.py
index 1c3bacc1b..0a7a051bc 100644
--- a/mcdc/object_/technique.py
+++ b/mcdc/object_/technique.py
@@ -39,6 +39,24 @@ def __call__(self, active: bool = True, weight_target: float = 1.0):
         self.weight_target = weight_target
 
 
+# ======================================================================================
+# Weight split 
+# ======================================================================================
+
+
+class WeightSplit(ObjectSingleton):
+    # Annotations for Numba mode
+    label: str = "weight_split"
+
+    weight_threshold: float
+
+    def __init__(self):
+        self.weight_threshold = 1.0
+
+    def __call__(self, weight_threshold: float = 1.0):
+        self.weight_threshold = weight_threshold
+
+
 # ======================================================================================
 # Weight roulette
 # ======================================================================================

From 1ce5146a194c204c6aeffa5ee92cd9a8c4a26c22 Mon Sep 17 00:00:00 2001
From: Nathan Glaser <glaserna@oregonstate.edu>
Date: Wed, 1 Apr 2026 21:52:11 -0700
Subject: [PATCH 03/53] remove initial split, add weight window mockup

---
 mcdc/object_/technique.py   | 54 ++++++++++++++++++++++++-------------
 mcdc/transport/technique.py | 53 +++++++++++++++++++++++++++---------
 2 files changed, 75 insertions(+), 32 deletions(-)

diff --git a/mcdc/object_/technique.py b/mcdc/object_/technique.py
index 0a7a051bc..f8848863a 100644
--- a/mcdc/object_/technique.py
+++ b/mcdc/object_/technique.py
@@ -1,4 +1,7 @@
-from mcdc.object_.base import ObjectSingleton
+from dataclasses import dataclass
+from numpy.typing import NDArray
+from mcdc.object_.base import ObjectBase, ObjectSingleton
+from mcdc.object_.mesh import MeshBase
 from mcdc.print_ import print_error
 
 # ======================================================================================
@@ -39,24 +42,6 @@ def __call__(self, active: bool = True, weight_target: float = 1.0):
         self.weight_target = weight_target
 
 
-# ======================================================================================
-# Weight split 
-# ======================================================================================
-
-
-class WeightSplit(ObjectSingleton):
-    # Annotations for Numba mode
-    label: str = "weight_split"
-
-    weight_threshold: float
-
-    def __init__(self):
-        self.weight_threshold = 1.0
-
-    def __call__(self, weight_threshold: float = 1.0):
-        self.weight_threshold = weight_threshold
-
-
 # ======================================================================================
 # Weight roulette
 # ======================================================================================
@@ -82,6 +67,37 @@ def __call__(self, weight_threshold: float = 0.0, weight_target: float = 1.0):
         self.weight_target = weight_target
 
 
+# ======================================================================================
+# Weight Windows
+# ======================================================================================
+
+
+@dataclass
+class WeightWindow(ObjectBase):
+    label: str = "weight_window"
+    upper_bound: float = 1.0
+    lower_bound: float = 0.0
+    target_weight: float = 1.0
+
+
+class WeightWindows(ObjectSingleton):
+    label: str = "weight_windows"
+
+    mesh: MeshBase
+    # 3d array of datatype WeightWindow
+    weight_windows: NDArray[tuple[int, int, int], WeightWindow]
+
+    def __call__(self, mesh, weight_windows):
+        ww_shape = weight_windows.shape
+        mesh_shape = (mesh.Nx, mesh.Ny, mesh.Nz)
+        if ww_shape != mesh_shape:
+            print_error(
+                f"Weight window array has shape {ww_shape}, but expected {mesh_shape}"
+            )
+        self.mesh = mesh
+        self.weight_windows = weight_windows
+
+
 # ======================================================================================
 # Population control
 # ======================================================================================
diff --git a/mcdc/transport/technique.py b/mcdc/transport/technique.py
index 340e73c14..26891590b 100644
--- a/mcdc/transport/technique.py
+++ b/mcdc/transport/technique.py
@@ -6,37 +6,64 @@
 ####
 
 import mcdc.numba_types as type_
+from mcdc.transport.mesh import get_indices as get_mesh_indices
 import mcdc.transport.particle as particle_module
 import mcdc.transport.particle_bank as particle_bank_module
 import mcdc.transport.rng as rng
 
+
 # ======================================================================================
-# Weight Splitting
+# Weight Roulette
 # ======================================================================================
 
 
 @njit
-def weight_split(particle_container, mcdc):
+def weight_roulette(particle_container, mcdc):
     particle = particle_container[0]
-    weight = particle["w"]
-    threshold = mcdc["weight_split"]["weight_threshold"]
-    if weight > threshold:
-        num_split = math.ceil(weight / threshold)
-        particle["w"] = weight / num_split
-        for _ in range(num_split - 1):
-            particle_bank_module.bank_active_particle(particle_container, mcdc)
+    if particle["w"] < mcdc["weight_roulette"]["weight_threshold"]:
+        w_target = mcdc["weight_roulette"]["weight_target"]
+        survival_probability = particle["w"] / w_target
+        if rng.lcg(particle_container) < survival_probability:
+            particle["w"] = w_target
+        else:
+            particle["alive"] = False
 
 
 # ======================================================================================
-# Weight Roulette
+# Weight Windows
 # ======================================================================================
 
 
 @njit
-def weight_roulette(particle_container, mcdc):
+def weight_windows(particle_container, mcdc, data):
+    weight_window = query_weight_window(particle_container, mcdc, data)
+    split_from_weight_window(particle_container, weight_window["upper_bound"], mcdc)
+    roulette_from_weight_window(particle_container, weight_window["lower_bound"], weight_window["target_weight"])
+
+
+@njit
+def query_weight_window(particle_container, mcdc, data):
+    mesh = mcdc["weight_windows"]["mesh"]
+    ww_array = mcdc["weight_windows"]["weight_windows"]
+    idx, idy, idz = get_mesh_indices(particle_container, mesh, mcdc, data)
+    return ww_array[idx, idy, idz]
+
+
+@njit
+def split_from_weight_window(particle_container, threshold_weight, mcdc):
     particle = particle_container[0]
-    if particle["w"] < mcdc["weight_roulette"]["weight_threshold"]:
-        w_target = mcdc["weight_roulette"]["weight_target"]
+    weight = particle["w"]
+    if weight > threshold_weight:
+        num_split = math.ceil(weight / threshold_weight)
+        particle["w"] = weight / num_split
+        for _ in range(num_split - 1):
+            particle_bank_module.bank_active_particle(particle_container, mcdc)
+
+
+@njit
+def roulette_from_weight_window(particle_container, w_threshold, w_target):
+    particle = particle_container[0]
+    if particle["w"] < w_threshold:
         survival_probability = particle["w"] / w_target
         if rng.lcg(particle_container) < survival_probability:
             particle["w"] = w_target

From 3d3e239bd9e93494ceeca4ae892c3d36b5143e9f Mon Sep 17 00:00:00 2001
From: Nathan Glaser <glaserna@oregonstate.edu>
Date: Wed, 1 Apr 2026 22:00:02 -0700
Subject: [PATCH 04/53] add notion of active for weight windows, fully removed
 initial split mockup

---
 mcdc/object_/simulation.py   | 6 +++---
 mcdc/object_/technique.py    | 5 +++++
 mcdc/transport/simulation.py | 5 +++--
 3 files changed, 11 insertions(+), 5 deletions(-)

diff --git a/mcdc/object_/simulation.py b/mcdc/object_/simulation.py
index 3acd888a4..6b75fc438 100644
--- a/mcdc/object_/simulation.py
+++ b/mcdc/object_/simulation.py
@@ -4,8 +4,8 @@
 from mcdc.object_.technique import (
     ImplicitCapture,
     PopulationControl,
-    WeightSplit,
     WeightRoulette,
+    WeightWindows,
     WeightedEmission,
 )
 
@@ -78,8 +78,8 @@ class Simulation(ObjectSingleton):
     # Techniques
     implicit_capture: ImplicitCapture
     weighted_emission: WeightedEmission
-    weight_split: WeightSplit
     weight_roulette: WeightRoulette
+    weight_windows: WeightWindows
     population_control: PopulationControl
 
     # Particle banks
@@ -161,8 +161,8 @@ def __init__(self):
         # Techniques
         self.implicit_capture = ImplicitCapture()
         self.weighted_emission = WeightedEmission()
-        self.weight_split = WeightSplit()
         self.weight_roulette = WeightRoulette()
+        self.weight_windows = WeightWindows()
         self.population_control = PopulationControl()
 
         # ==============================================================================
diff --git a/mcdc/object_/technique.py b/mcdc/object_/technique.py
index f8848863a..2ed660e4d 100644
--- a/mcdc/object_/technique.py
+++ b/mcdc/object_/technique.py
@@ -83,10 +83,14 @@ class WeightWindow(ObjectBase):
 class WeightWindows(ObjectSingleton):
     label: str = "weight_windows"
 
+    active: bool
     mesh: MeshBase
     # 3d array of datatype WeightWindow
     weight_windows: NDArray[tuple[int, int, int], WeightWindow]
 
+    def __init__(self):
+        self.active = False
+
     def __call__(self, mesh, weight_windows):
         ww_shape = weight_windows.shape
         mesh_shape = (mesh.Nx, mesh.Ny, mesh.Nz)
@@ -94,6 +98,7 @@ def __call__(self, mesh, weight_windows):
             print_error(
                 f"Weight window array has shape {ww_shape}, but expected {mesh_shape}"
             )
+        self.active = True
         self.mesh = mesh
         self.weight_windows = weight_windows
 
diff --git a/mcdc/transport/simulation.py b/mcdc/transport/simulation.py
index 796e32beb..cafd5a008 100644
--- a/mcdc/transport/simulation.py
+++ b/mcdc/transport/simulation.py
@@ -306,8 +306,9 @@ def step_particle(particle_container, mcdc, data):
         particle["alive"] = False
 
     # Weight splitting / rouletting
-    if particle["alive"]:
-        technique.weight_split(particle_container, mcdc)
+    if particle["alive"] & mcdc.weight_windows["active"]:
+        technique.weight_windows(particle_container, mcdc, data)
+    elif particle["alive"]:
         technique.weight_roulette(particle_container, mcdc)
 
 

From 520986de371b5122b58984dd090c83e745e0300d Mon Sep 17 00:00:00 2001
From: Nathan Glaser <glaserna@oregonstate.edu>
Date: Sun, 12 Apr 2026 18:22:12 -0700
Subject: [PATCH 05/53] switch to 1darray from nested, fix mesh get

---
 mcdc/object_/technique.py    | 41 ++++++++++++++++++++++++------------
 mcdc/transport/simulation.py |  2 +-
 mcdc/transport/technique.py  | 15 +++++++------
 3 files changed, 36 insertions(+), 22 deletions(-)

diff --git a/mcdc/object_/technique.py b/mcdc/object_/technique.py
index 2ed660e4d..535dbb1ff 100644
--- a/mcdc/object_/technique.py
+++ b/mcdc/object_/technique.py
@@ -1,5 +1,6 @@
 from dataclasses import dataclass
 from numpy.typing import NDArray
+from numpy import float64
 from mcdc.object_.base import ObjectBase, ObjectSingleton
 from mcdc.object_.mesh import MeshBase
 from mcdc.print_ import print_error
@@ -72,35 +73,47 @@ def __call__(self, weight_threshold: float = 0.0, weight_target: float = 1.0):
 # ======================================================================================
 
 
-@dataclass
-class WeightWindow(ObjectBase):
-    label: str = "weight_window"
-    upper_bound: float = 1.0
-    lower_bound: float = 0.0
-    target_weight: float = 1.0
-
-
 class WeightWindows(ObjectSingleton):
     label: str = "weight_windows"
 
     active: bool
     mesh: MeshBase
-    # 3d array of datatype WeightWindow
-    weight_windows: NDArray[tuple[int, int, int], WeightWindow]
+    Nx: int
+    Ny: int
+    Nz: int
+
+    # flattened array of ww
+    weight_windows: NDArray[float64]
 
     def __init__(self):
         self.active = False
 
     def __call__(self, mesh, weight_windows):
+        # get mesh size
+        match mesh.label:
+            case "uniform_mesh":
+                nx, ny, nz = mesh.Nx, mesh.Ny, mesh.Nz
+            case "structured_mesh":
+                nx, ny, nz = mesh.x.shape[0], mesh.y.shape[0], mesh.z.shape[0]
+            case _:
+                print_error(
+                    f"{type(mesh).__name__} is not supported for weight windows"
+                )
+        
+        mesh_shape = (nx, ny, nz)
         ww_shape = weight_windows.shape
-        mesh_shape = (mesh.Nx, mesh.Ny, mesh.Nz)
-        if ww_shape != mesh_shape:
+        expected_shape = (*mesh_shape, 3)
+        if ww_shape != expected_shape:
             print_error(
-                f"Weight window array has shape {ww_shape}, but expected {mesh_shape}"
+                f"Weight window array has shape {ww_shape}, but expected {expected_shape}"
             )
+
         self.active = True
         self.mesh = mesh
-        self.weight_windows = weight_windows
+        self.Nx = nx
+        self.Ny = ny
+        self.Nz = nz
+        self.weight_windows = weight_windows.reshape(-1)
 
 
 # ======================================================================================
diff --git a/mcdc/transport/simulation.py b/mcdc/transport/simulation.py
index cafd5a008..89b18fef8 100644
--- a/mcdc/transport/simulation.py
+++ b/mcdc/transport/simulation.py
@@ -306,7 +306,7 @@ def step_particle(particle_container, mcdc, data):
         particle["alive"] = False
 
     # Weight splitting / rouletting
-    if particle["alive"] & mcdc.weight_windows["active"]:
+    if particle["alive"] & mcdc["weight_windows"]["active"]:
         technique.weight_windows(particle_container, mcdc, data)
     elif particle["alive"]:
         technique.weight_roulette(particle_container, mcdc)
diff --git a/mcdc/transport/technique.py b/mcdc/transport/technique.py
index 26891590b..56a9b49ec 100644
--- a/mcdc/transport/technique.py
+++ b/mcdc/transport/technique.py
@@ -4,7 +4,7 @@
 from numba import njit
 
 ####
-
+from mcdc.mcdc_get.weight_windows import weight_windows_all, weight_windows_chunk
 import mcdc.numba_types as type_
 from mcdc.transport.mesh import get_indices as get_mesh_indices
 import mcdc.transport.particle as particle_module
@@ -36,17 +36,18 @@ def weight_roulette(particle_container, mcdc):
 
 @njit
 def weight_windows(particle_container, mcdc, data):
-    weight_window = query_weight_window(particle_container, mcdc, data)
-    split_from_weight_window(particle_container, weight_window["upper_bound"], mcdc)
-    roulette_from_weight_window(particle_container, weight_window["lower_bound"], weight_window["target_weight"])
+    [lower, target, upper] = query_weight_window(particle_container, mcdc, data)
+    split_from_weight_window(particle_container, upper, mcdc)
+    roulette_from_weight_window(particle_container, lower, target)
 
 
 @njit
 def query_weight_window(particle_container, mcdc, data):
-    mesh = mcdc["weight_windows"]["mesh"]
-    ww_array = mcdc["weight_windows"]["weight_windows"]
+    ww_obj = mcdc["weight_windows"]
+    mesh = mcdc["meshes"][ww_obj["mesh_ID"]]
     idx, idy, idz = get_mesh_indices(particle_container, mesh, mcdc, data)
-    return ww_array[idx, idy, idz]
+    start = 3*(idx * mesh["Nx"] + idy * mesh["Ny"] + idz * mesh["Nz"])
+    return weight_windows_chunk(start, 3, ww_obj, data)
 
 
 @njit

From f9d6940cc6a4680aaf09f82c02bfe6b8a5ceb959 Mon Sep 17 00:00:00 2001
From: Nathan Glaser <glaserna@oregonstate.edu>
Date: Sun, 12 Apr 2026 18:22:30 -0700
Subject: [PATCH 06/53] generate mcdc_get and mcdc_set

---
 mcdc/mcdc_get/__init__.py       |  2 ++
 mcdc/mcdc_get/weight_windows.py | 32 ++++++++++++++++++++++++++++++++
 mcdc/mcdc_set/__init__.py       |  2 ++
 mcdc/mcdc_set/weight_windows.py | 32 ++++++++++++++++++++++++++++++++
 4 files changed, 68 insertions(+)
 create mode 100644 mcdc/mcdc_get/weight_windows.py
 create mode 100644 mcdc/mcdc_set/weight_windows.py

diff --git a/mcdc/mcdc_get/__init__.py b/mcdc/mcdc_get/__init__.py
index 0f51357ba..8138c4763 100644
--- a/mcdc/mcdc_get/__init__.py
+++ b/mcdc/mcdc_get/__init__.py
@@ -78,6 +78,8 @@
 
 import mcdc.mcdc_get.weight_roulette as weight_roulette
 
+import mcdc.mcdc_get.weight_windows as weight_windows
+
 import mcdc.mcdc_get.weighted_emission as weighted_emission
 
 import mcdc.mcdc_get.source as source
diff --git a/mcdc/mcdc_get/weight_windows.py b/mcdc/mcdc_get/weight_windows.py
new file mode 100644
index 000000000..143d55fc2
--- /dev/null
+++ b/mcdc/mcdc_get/weight_windows.py
@@ -0,0 +1,32 @@
+# The following is automatically generated by code_factory.py
+
+from numba import njit
+
+
+@njit
+def weight_windows(index, weight_windows, data):
+    offset = weight_windows["weight_windows_offset"]
+    return data[offset + index]
+
+
+@njit
+def weight_windows_all(weight_windows, data):
+    start = weight_windows["weight_windows_offset"]
+    size = weight_windows["weight_windows_length"]
+    end = start + size
+    return data[start:end]
+
+
+@njit
+def weight_windows_last(weight_windows, data):
+    start = weight_windows["weight_windows_offset"]
+    size = weight_windows["weight_windows_length"]
+    end = start + size
+    return data[end - 1]
+
+
+@njit
+def weight_windows_chunk(start, length, weight_windows, data):
+    start += weight_windows["weight_windows_offset"]
+    end = start + length
+    return data[start:end]
diff --git a/mcdc/mcdc_set/__init__.py b/mcdc/mcdc_set/__init__.py
index 840de390b..1d2486ed7 100644
--- a/mcdc/mcdc_set/__init__.py
+++ b/mcdc/mcdc_set/__init__.py
@@ -78,6 +78,8 @@
 
 import mcdc.mcdc_set.weight_roulette as weight_roulette
 
+import mcdc.mcdc_set.weight_windows as weight_windows
+
 import mcdc.mcdc_set.weighted_emission as weighted_emission
 
 import mcdc.mcdc_set.source as source
diff --git a/mcdc/mcdc_set/weight_windows.py b/mcdc/mcdc_set/weight_windows.py
new file mode 100644
index 000000000..9f599f3c1
--- /dev/null
+++ b/mcdc/mcdc_set/weight_windows.py
@@ -0,0 +1,32 @@
+# The following is automatically generated by code_factory.py
+
+from numba import njit
+
+
+@njit
+def weight_windows(index, weight_windows, data, value):
+    offset = weight_windows["weight_windows_offset"]
+    data[offset + index] = value
+
+
+@njit
+def weight_windows_all(weight_windows, data, value):
+    start = weight_windows["weight_windows_offset"]
+    size = weight_windows["weight_windows_length"]
+    end = start + size
+    data[start:end] = value
+
+
+@njit
+def weight_windows_last(weight_windows, data, value):
+    start = weight_windows["weight_windows_offset"]
+    size = weight_windows["weight_windows_length"]
+    end = start + size
+    data[end - 1] = value
+
+
+@njit
+def weight_windows_chunk(start, length, weight_windows, data, value):
+    start += weight_windows["weight_windows_offset"]
+    end = start + length
+    data[start:end] = value

From 30357378ccc400b75dfa1045520784798c59f87c Mon Sep 17 00:00:00 2001
From: Nathan Glaser <glaserna@oregonstate.edu>
Date: Sun, 12 Apr 2026 18:46:03 -0700
Subject: [PATCH 07/53] switch to tracking different arrays of lower, target,
 upper

---
 mcdc/mcdc_get/weight_windows.py | 78 ++++++++++++++++++++++++++++-----
 mcdc/mcdc_set/weight_windows.py | 78 ++++++++++++++++++++++++++++-----
 mcdc/object_/technique.py       | 16 ++++---
 mcdc/transport/technique.py     |  9 ++--
 4 files changed, 151 insertions(+), 30 deletions(-)

diff --git a/mcdc/mcdc_get/weight_windows.py b/mcdc/mcdc_get/weight_windows.py
index 143d55fc2..00252e542 100644
--- a/mcdc/mcdc_get/weight_windows.py
+++ b/mcdc/mcdc_get/weight_windows.py
@@ -4,29 +4,87 @@
 
 
 @njit
-def weight_windows(index, weight_windows, data):
-    offset = weight_windows["weight_windows_offset"]
+def lower_weights(index, weight_windows, data):
+    offset = weight_windows["lower_weights_offset"]
     return data[offset + index]
 
 
 @njit
-def weight_windows_all(weight_windows, data):
-    start = weight_windows["weight_windows_offset"]
-    size = weight_windows["weight_windows_length"]
+def lower_weights_all(weight_windows, data):
+    start = weight_windows["lower_weights_offset"]
+    size = weight_windows["lower_weights_length"]
     end = start + size
     return data[start:end]
 
 
 @njit
-def weight_windows_last(weight_windows, data):
-    start = weight_windows["weight_windows_offset"]
-    size = weight_windows["weight_windows_length"]
+def lower_weights_last(weight_windows, data):
+    start = weight_windows["lower_weights_offset"]
+    size = weight_windows["lower_weights_length"]
     end = start + size
     return data[end - 1]
 
 
 @njit
-def weight_windows_chunk(start, length, weight_windows, data):
-    start += weight_windows["weight_windows_offset"]
+def lower_weights_chunk(start, length, weight_windows, data):
+    start += weight_windows["lower_weights_offset"]
+    end = start + length
+    return data[start:end]
+
+
+@njit
+def target_weights(index, weight_windows, data):
+    offset = weight_windows["target_weights_offset"]
+    return data[offset + index]
+
+
+@njit
+def target_weights_all(weight_windows, data):
+    start = weight_windows["target_weights_offset"]
+    size = weight_windows["target_weights_length"]
+    end = start + size
+    return data[start:end]
+
+
+@njit
+def target_weights_last(weight_windows, data):
+    start = weight_windows["target_weights_offset"]
+    size = weight_windows["target_weights_length"]
+    end = start + size
+    return data[end - 1]
+
+
+@njit
+def target_weights_chunk(start, length, weight_windows, data):
+    start += weight_windows["target_weights_offset"]
+    end = start + length
+    return data[start:end]
+
+
+@njit
+def upper_weights(index, weight_windows, data):
+    offset = weight_windows["upper_weights_offset"]
+    return data[offset + index]
+
+
+@njit
+def upper_weights_all(weight_windows, data):
+    start = weight_windows["upper_weights_offset"]
+    size = weight_windows["upper_weights_length"]
+    end = start + size
+    return data[start:end]
+
+
+@njit
+def upper_weights_last(weight_windows, data):
+    start = weight_windows["upper_weights_offset"]
+    size = weight_windows["upper_weights_length"]
+    end = start + size
+    return data[end - 1]
+
+
+@njit
+def upper_weights_chunk(start, length, weight_windows, data):
+    start += weight_windows["upper_weights_offset"]
     end = start + length
     return data[start:end]
diff --git a/mcdc/mcdc_set/weight_windows.py b/mcdc/mcdc_set/weight_windows.py
index 9f599f3c1..95613786f 100644
--- a/mcdc/mcdc_set/weight_windows.py
+++ b/mcdc/mcdc_set/weight_windows.py
@@ -4,29 +4,87 @@
 
 
 @njit
-def weight_windows(index, weight_windows, data, value):
-    offset = weight_windows["weight_windows_offset"]
+def lower_weights(index, weight_windows, data, value):
+    offset = weight_windows["lower_weights_offset"]
     data[offset + index] = value
 
 
 @njit
-def weight_windows_all(weight_windows, data, value):
-    start = weight_windows["weight_windows_offset"]
-    size = weight_windows["weight_windows_length"]
+def lower_weights_all(weight_windows, data, value):
+    start = weight_windows["lower_weights_offset"]
+    size = weight_windows["lower_weights_length"]
     end = start + size
     data[start:end] = value
 
 
 @njit
-def weight_windows_last(weight_windows, data, value):
-    start = weight_windows["weight_windows_offset"]
-    size = weight_windows["weight_windows_length"]
+def lower_weights_last(weight_windows, data, value):
+    start = weight_windows["lower_weights_offset"]
+    size = weight_windows["lower_weights_length"]
     end = start + size
     data[end - 1] = value
 
 
 @njit
-def weight_windows_chunk(start, length, weight_windows, data, value):
-    start += weight_windows["weight_windows_offset"]
+def lower_weights_chunk(start, length, weight_windows, data, value):
+    start += weight_windows["lower_weights_offset"]
+    end = start + length
+    data[start:end] = value
+
+
+@njit
+def target_weights(index, weight_windows, data, value):
+    offset = weight_windows["target_weights_offset"]
+    data[offset + index] = value
+
+
+@njit
+def target_weights_all(weight_windows, data, value):
+    start = weight_windows["target_weights_offset"]
+    size = weight_windows["target_weights_length"]
+    end = start + size
+    data[start:end] = value
+
+
+@njit
+def target_weights_last(weight_windows, data, value):
+    start = weight_windows["target_weights_offset"]
+    size = weight_windows["target_weights_length"]
+    end = start + size
+    data[end - 1] = value
+
+
+@njit
+def target_weights_chunk(start, length, weight_windows, data, value):
+    start += weight_windows["target_weights_offset"]
+    end = start + length
+    data[start:end] = value
+
+
+@njit
+def upper_weights(index, weight_windows, data, value):
+    offset = weight_windows["upper_weights_offset"]
+    data[offset + index] = value
+
+
+@njit
+def upper_weights_all(weight_windows, data, value):
+    start = weight_windows["upper_weights_offset"]
+    size = weight_windows["upper_weights_length"]
+    end = start + size
+    data[start:end] = value
+
+
+@njit
+def upper_weights_last(weight_windows, data, value):
+    start = weight_windows["upper_weights_offset"]
+    size = weight_windows["upper_weights_length"]
+    end = start + size
+    data[end - 1] = value
+
+
+@njit
+def upper_weights_chunk(start, length, weight_windows, data, value):
+    start += weight_windows["upper_weights_offset"]
     end = start + length
     data[start:end] = value
diff --git a/mcdc/object_/technique.py b/mcdc/object_/technique.py
index 535dbb1ff..5ebec3d96 100644
--- a/mcdc/object_/technique.py
+++ b/mcdc/object_/technique.py
@@ -82,8 +82,10 @@ class WeightWindows(ObjectSingleton):
     Ny: int
     Nz: int
 
-    # flattened array of ww
-    weight_windows: NDArray[float64]
+    # flattened arrays of ww params
+    lower_weights: NDArray[float64]
+    target_weights: NDArray[float64]
+    upper_weights: NDArray[float64]
 
     def __init__(self):
         self.active = False
@@ -94,7 +96,7 @@ def __call__(self, mesh, weight_windows):
             case "uniform_mesh":
                 nx, ny, nz = mesh.Nx, mesh.Ny, mesh.Nz
             case "structured_mesh":
-                nx, ny, nz = mesh.x.shape[0], mesh.y.shape[0], mesh.z.shape[0]
+                nx, ny, nz = mesh.x.shape[0]-1, mesh.y.shape[0]-1, mesh.z.shape[0]-1
             case _:
                 print_error(
                     f"{type(mesh).__name__} is not supported for weight windows"
@@ -110,10 +112,10 @@ def __call__(self, mesh, weight_windows):
 
         self.active = True
         self.mesh = mesh
-        self.Nx = nx
-        self.Ny = ny
-        self.Nz = nz
-        self.weight_windows = weight_windows.reshape(-1)
+        self.Nx, self.Ny, self.Nz = mesh_shape
+        self.lower_weights = weight_windows[..., 0].reshape(-1)
+        self.target_weights = weight_windows[..., 1].reshape(-1)
+        self.upper_weights = weight_windows[..., 2].reshape(-1)
 
 
 # ======================================================================================
diff --git a/mcdc/transport/technique.py b/mcdc/transport/technique.py
index 56a9b49ec..b44a6af85 100644
--- a/mcdc/transport/technique.py
+++ b/mcdc/transport/technique.py
@@ -4,7 +4,7 @@
 from numba import njit
 
 ####
-from mcdc.mcdc_get.weight_windows import weight_windows_all, weight_windows_chunk
+import mcdc.mcdc_get.weight_windows as ww_get
 import mcdc.numba_types as type_
 from mcdc.transport.mesh import get_indices as get_mesh_indices
 import mcdc.transport.particle as particle_module
@@ -46,8 +46,11 @@ def query_weight_window(particle_container, mcdc, data):
     ww_obj = mcdc["weight_windows"]
     mesh = mcdc["meshes"][ww_obj["mesh_ID"]]
     idx, idy, idz = get_mesh_indices(particle_container, mesh, mcdc, data)
-    start = 3*(idx * mesh["Nx"] + idy * mesh["Ny"] + idz * mesh["Nz"])
-    return weight_windows_chunk(start, 3, ww_obj, data)
+    index = ((idx * ww_obj["Ny"]) + idy) * ww_obj["Nz"] + idz
+    lower = ww_get.lower_weights(index, ww_obj, data)
+    target = ww_get.target_weights(index, ww_obj, data)
+    upper = ww_get.upper_weights(index, ww_obj, data) 
+    return lower, target, upper
 
 
 @njit

From 5707fe187852022951e6fe089989b68a33e40999 Mon Sep 17 00:00:00 2001
From: Nathan Glaser <glaserna@oregonstate.edu>
Date: Sun, 12 Apr 2026 18:57:35 -0700
Subject: [PATCH 08/53] add checks for weight window validity

---
 mcdc/object_/technique.py | 17 +++++++++++++++++
 1 file changed, 17 insertions(+)

diff --git a/mcdc/object_/technique.py b/mcdc/object_/technique.py
index 5ebec3d96..137363d0e 100644
--- a/mcdc/object_/technique.py
+++ b/mcdc/object_/technique.py
@@ -1,6 +1,7 @@
 from dataclasses import dataclass
 from numpy.typing import NDArray
 from numpy import float64
+import numpy as np
 from mcdc.object_.base import ObjectBase, ObjectSingleton
 from mcdc.object_.mesh import MeshBase
 from mcdc.print_ import print_error
@@ -102,6 +103,7 @@ def __call__(self, mesh, weight_windows):
                     f"{type(mesh).__name__} is not supported for weight windows"
                 )
         
+        # check correct shape
         mesh_shape = (nx, ny, nz)
         ww_shape = weight_windows.shape
         expected_shape = (*mesh_shape, 3)
@@ -117,6 +119,21 @@ def __call__(self, mesh, weight_windows):
         self.target_weights = weight_windows[..., 1].reshape(-1)
         self.upper_weights = weight_windows[..., 2].reshape(-1)
 
+        # check weight windows are valid
+        if (self.lower_weights <= 0.0).any():
+            print_error(
+                "Lower bound weights must be strictly positive to avoid invalid roulette behavior"
+            )
+        if (self.lower_weights > self.target_weights).any():
+            print_error(
+                f"Lower bound weight can not be greater than the target weight for any weight window"
+            )
+        if (self.target_weights > self.upper_weights).any():
+            print_error(
+                f"Target weight can not be greater than the upper bound weight for any weight window"
+            )
+
+
 
 # ======================================================================================
 # Population control

From e9a6afa13384c3f6450eabca3a903ac24d7441d7 Mon Sep 17 00:00:00 2001
From: Nathan Glaser <glaserna@oregonstate.edu>
Date: Sun, 12 Apr 2026 19:49:42 -0700
Subject: [PATCH 09/53] add test for instantiation

---
 mcdc/object_/technique.py                     |  4 +-
 .../transport/technique/weight_windows.py     | 84 +++++++++++++++++++
 2 files changed, 86 insertions(+), 2 deletions(-)
 create mode 100644 test/unit/transport/technique/weight_windows.py

diff --git a/mcdc/object_/technique.py b/mcdc/object_/technique.py
index 137363d0e..4017f3b73 100644
--- a/mcdc/object_/technique.py
+++ b/mcdc/object_/technique.py
@@ -126,11 +126,11 @@ def __call__(self, mesh, weight_windows):
             )
         if (self.lower_weights > self.target_weights).any():
             print_error(
-                f"Lower bound weight can not be greater than the target weight for any weight window"
+                "Lower bound weight can not be greater than the target weight for any weight window"
             )
         if (self.target_weights > self.upper_weights).any():
             print_error(
-                f"Target weight can not be greater than the upper bound weight for any weight window"
+                "Target weight can not be greater than the upper bound weight for any weight window"
             )
 
 
diff --git a/test/unit/transport/technique/weight_windows.py b/test/unit/transport/technique/weight_windows.py
new file mode 100644
index 000000000..fa39a0d9a
--- /dev/null
+++ b/test/unit/transport/technique/weight_windows.py
@@ -0,0 +1,84 @@
+import mcdc
+import numpy as np
+import pytest
+import os
+
+from mcdc.main import preparation
+
+def make_ww_model(lower=0.1, target=1.0, upper=1.0, mess_up_size=False):
+  # Geometry
+  # surfaces
+  cylinder = mcdc.Surface.CylinderZ(radius=0.5)
+  pitch = 2.0
+  height = 10.0
+  x0 = mcdc.Surface.PlaneX(x=-pitch / 2, boundary_condition="reflective")
+  x1 = mcdc.Surface.PlaneX(x=pitch / 2, boundary_condition="reflective")
+  y0 = mcdc.Surface.PlaneY(y=-pitch / 2, boundary_condition="reflective")
+  y1 = mcdc.Surface.PlaneY(y=pitch / 2, boundary_condition="reflective")
+  z0 = mcdc.Surface.PlaneZ(z=0.0, boundary_condition="reflective")
+  z1 = mcdc.Surface.PlaneZ(z=height, boundary_condition="reflective")
+
+  # cells
+  mcdc.Cell(-cylinder)
+  mcdc.Cell(+x0 & -x1 & +y0 & -y1 & +z0 & -z1 & +cylinder)
+
+  # Source
+  mcdc.Source(position=[0.0, 0.0, 0.0], isotropic=True, time=0.0, energy=14.1e6)
+
+  # Setting
+  mcdc.settings.N_particle = 20
+  mcdc.settings.N_batch = 2
+  mcdc.settings.time_boundary = 1.0
+  mcdc.settings.active_bank_buffer = 1000
+
+  # weight windows
+  N = 3
+  mesh = mcdc.MeshUniform(
+    "mesh",
+    x=(-pitch/2, pitch/2, N),
+    y=(-pitch/2, pitch/2, N),
+    z=(0.0, height, N)
+  )
+  if mess_up_size:
+    shape = (N, N, 4, 3)
+  else:
+    shape = (N, N, N, 3)
+  ww_array = np.ones(shape)
+  ww_array[:, :, :, 0] = lower
+  ww_array[:, :, :, 1] = target
+  ww_array[:, :, :, 2] = upper
+  mcdc.simulation.weight_windows(mesh, ww_array)
+
+  return preparation()
+
+@pytest.mark.parametrize(
+    "kwargs, expected_msg",
+    [
+        # incorrect size
+        (
+            {"mess_up_size": True},
+            "Weight window array has shape (3, 3, 4, 3), but expected (3, 3, 3, 3)",
+        ),
+        # negative lower
+        (
+            {"lower": -1.0},
+            "Lower bound weights must be strictly positive",
+        ),
+        # lower > target
+        (
+            {"lower": 1.0, "target": 0.5},
+            "Lower bound weight can not be greater than the target weight",
+        ),
+        # target > upper
+        (
+            {"target": 1.5},
+            "Target weight can not be greater than the upper bound weight",
+        ),
+    ],
+)
+def test_error_throw(capsys, kwargs, expected_msg):
+  with pytest.raises(SystemExit):
+    make_ww_model(**kwargs)
+
+  out = capsys.readouterr().out
+  assert expected_msg in out
\ No newline at end of file

From 28d9f2cfe219ee09790b81f3258ce31044b41d1c Mon Sep 17 00:00:00 2001
From: Nathan Glaser <glaserna@oregonstate.edu>
Date: Sun, 12 Apr 2026 20:22:54 -0700
Subject: [PATCH 10/53] add test for ww roulette method

---
 .../transport/technique/weight_windows.py     | 41 ++++++++++++++++++-
 1 file changed, 40 insertions(+), 1 deletion(-)

diff --git a/test/unit/transport/technique/weight_windows.py b/test/unit/transport/technique/weight_windows.py
index fa39a0d9a..20da77cab 100644
--- a/test/unit/transport/technique/weight_windows.py
+++ b/test/unit/transport/technique/weight_windows.py
@@ -4,6 +4,21 @@
 import os
 
 from mcdc.main import preparation
+import mcdc.numba_types as type_
+from mcdc.transport.technique import (
+  roulette_from_weight_window,
+  split_from_weight_window,
+  query_weight_window,
+  weight_windows
+)
+from mcdc.constant import (
+  TINY
+)
+
+# =========================================================================== #
+# Helper method for creating a dummy model
+# =========================================================================== #
+
 
 def make_ww_model(lower=0.1, target=1.0, upper=1.0, mess_up_size=False):
   # Geometry
@@ -51,6 +66,12 @@ def make_ww_model(lower=0.1, target=1.0, upper=1.0, mess_up_size=False):
 
   return preparation()
 
+
+# =========================================================================== #
+# Test error throwing in object creation
+# =========================================================================== #
+
+
 @pytest.mark.parametrize(
     "kwargs, expected_msg",
     [
@@ -81,4 +102,22 @@ def test_error_throw(capsys, kwargs, expected_msg):
     make_ww_model(**kwargs)
 
   out = capsys.readouterr().out
-  assert expected_msg in out
\ No newline at end of file
+  assert expected_msg in out
+
+
+# =========================================================================== #
+# Tests for helper methods 
+# =========================================================================== #
+
+
+def test_roullete_from_weight_window():
+  # because of rng, want to loop over to hit both branches
+  for i in range(10):
+    particles = np.zeros(1, type_.particle_data)
+    particles[0]["w"] = 0.1
+    target = 0.2
+    threshold = 0.1 + TINY
+    
+    roulette_from_weight_window(particles, threshold, target)
+    p = particles[0]
+    assert (p["w"] == target or not p["alive"])

From e832317492534c207b274c889ede07f41ec13e7b Mon Sep 17 00:00:00 2001
From: Nathan Glaser <glaserna@oregonstate.edu>
Date: Sun, 12 Apr 2026 20:50:11 -0700
Subject: [PATCH 11/53] Add split test

---
 .../transport/technique/weight_windows.py     | 33 +++++++++++++++++--
 1 file changed, 31 insertions(+), 2 deletions(-)

diff --git a/test/unit/transport/technique/weight_windows.py b/test/unit/transport/technique/weight_windows.py
index 20da77cab..63933d7cf 100644
--- a/test/unit/transport/technique/weight_windows.py
+++ b/test/unit/transport/technique/weight_windows.py
@@ -9,7 +9,8 @@
   roulette_from_weight_window,
   split_from_weight_window,
   query_weight_window,
-  weight_windows
+  weight_windows,
+  particle_bank_module
 )
 from mcdc.constant import (
   TINY
@@ -64,7 +65,8 @@ def make_ww_model(lower=0.1, target=1.0, upper=1.0, mess_up_size=False):
   ww_array[:, :, :, 2] = upper
   mcdc.simulation.weight_windows(mesh, ww_array)
 
-  return preparation()
+  mcdc_container, data = preparation()
+  return mcdc_container[0], data
 
 
 # =========================================================================== #
@@ -121,3 +123,30 @@ def test_roullete_from_weight_window():
     roulette_from_weight_window(particles, threshold, target)
     p = particles[0]
     assert (p["w"] == target or not p["alive"])
+
+
+def test_split_from_weight_window():
+  particles = np.zeros(1, type_.particle_data)
+  init_weight = 2.0 + TINY
+  particles[0]["w"] = init_weight
+  threshold = 1.0
+  mcdc_obj, data = make_ww_model()
+
+  # get bank and init size
+  bank = mcdc_obj["bank_active"]
+  init_bank_size = particle_bank_module.get_bank_size(bank)
+
+  # split
+  split_from_weight_window(particles, threshold, mcdc_obj)
+
+  p1 = particles[0]
+  num_split = np.ceil(init_weight / threshold)
+  num_new = num_split - 1
+
+  # check weight of original particle
+  assert p1["w"] == init_weight / num_split
+  assert particle_bank_module.get_bank_size(bank) == init_bank_size + num_new
+  for i in range(2):
+    pnew = bank["particles"][init_bank_size + i]
+    assert pnew["w"] == p1["w"]
+

From eb037e9392332accfac49e0095e996e02e13ba09 Mon Sep 17 00:00:00 2001
From: Nathan Glaser <glaserna@oregonstate.edu>
Date: Sun, 12 Apr 2026 21:37:39 -0700
Subject: [PATCH 12/53] Add query of weight windows test

---
 .../transport/technique/weight_windows.py     | 159 ++++++++++++------
 1 file changed, 111 insertions(+), 48 deletions(-)

diff --git a/test/unit/transport/technique/weight_windows.py b/test/unit/transport/technique/weight_windows.py
index 63933d7cf..7974d6910 100644
--- a/test/unit/transport/technique/weight_windows.py
+++ b/test/unit/transport/technique/weight_windows.py
@@ -12,6 +12,7 @@
   weight_windows,
   particle_bank_module
 )
+from mcdc.transport.mesh.interface import get_indices
 from mcdc.constant import (
   TINY
 )
@@ -20,53 +21,82 @@
 # Helper method for creating a dummy model
 # =========================================================================== #
 
+def make_base_model():
 
-def make_ww_model(lower=0.1, target=1.0, upper=1.0, mess_up_size=False):
-  # Geometry
-  # surfaces
-  cylinder = mcdc.Surface.CylinderZ(radius=0.5)
-  pitch = 2.0
-  height = 10.0
-  x0 = mcdc.Surface.PlaneX(x=-pitch / 2, boundary_condition="reflective")
-  x1 = mcdc.Surface.PlaneX(x=pitch / 2, boundary_condition="reflective")
-  y0 = mcdc.Surface.PlaneY(y=-pitch / 2, boundary_condition="reflective")
-  y1 = mcdc.Surface.PlaneY(y=pitch / 2, boundary_condition="reflective")
-  z0 = mcdc.Surface.PlaneZ(z=0.0, boundary_condition="reflective")
-  z1 = mcdc.Surface.PlaneZ(z=height, boundary_condition="reflective")
-
-  # cells
-  mcdc.Cell(-cylinder)
-  mcdc.Cell(+x0 & -x1 & +y0 & -y1 & +z0 & -z1 & +cylinder)
-
-  # Source
-  mcdc.Source(position=[0.0, 0.0, 0.0], isotropic=True, time=0.0, energy=14.1e6)
-
-  # Setting
-  mcdc.settings.N_particle = 20
-  mcdc.settings.N_batch = 2
-  mcdc.settings.time_boundary = 1.0
-  mcdc.settings.active_bank_buffer = 1000
-
-  # weight windows
-  N = 3
-  mesh = mcdc.MeshUniform(
-    "mesh",
-    x=(-pitch/2, pitch/2, N),
-    y=(-pitch/2, pitch/2, N),
-    z=(0.0, height, N)
-  )
-  if mess_up_size:
-    shape = (N, N, 4, 3)
-  else:
-    shape = (N, N, N, 3)
-  ww_array = np.ones(shape)
-  ww_array[:, :, :, 0] = lower
-  ww_array[:, :, :, 1] = target
-  ww_array[:, :, :, 2] = upper
-  mcdc.simulation.weight_windows(mesh, ww_array)
-
-  mcdc_container, data = preparation()
-  return mcdc_container[0], data
+    # Geometry
+    cylinder = mcdc.Surface.CylinderZ(radius=0.5)
+    pitch = 2.0
+    height = 10.0
+
+    x0 = mcdc.Surface.PlaneX(x=-pitch / 2, boundary_condition="reflective")
+    x1 = mcdc.Surface.PlaneX(x=pitch / 2, boundary_condition="reflective")
+    y0 = mcdc.Surface.PlaneY(y=-pitch / 2, boundary_condition="reflective")
+    y1 = mcdc.Surface.PlaneY(y=pitch / 2, boundary_condition="reflective")
+    z0 = mcdc.Surface.PlaneZ(z=0.0, boundary_condition="reflective")
+    z1 = mcdc.Surface.PlaneZ(z=height, boundary_condition="reflective")
+
+    mcdc.Cell(-cylinder)
+    mcdc.Cell(+x0 & -x1 & +y0 & -y1 & +z0 & -z1 & +cylinder)
+
+    # Source
+    mcdc.Source(position=[0.0, 0.0, 0.0], isotropic=True, time=0.0, energy=14.1e6)
+
+    # Settings
+    mcdc.settings.N_particle = 20
+    mcdc.settings.N_batch = 2
+    mcdc.settings.time_boundary = 1.0
+    mcdc.settings.active_bank_buffer = 1000
+
+    # Mesh
+    N = 3
+    mesh = mcdc.MeshUniform(
+        "mesh",
+        x=(-pitch/2, pitch/N, N),
+        y=(-pitch/2, pitch/N, N),
+        z=(0.0, height/N, N)
+    )
+
+    return mesh, N
+
+
+def make_ww_model_params(lower=0.1, target=1.0, upper=1.0, mess_up_size=False):
+    import mcdc
+
+    mesh, N = make_base_model()
+
+    if mess_up_size:
+        ww_array = np.ones((N, N, 4, 3))
+    else:
+        ww_array = np.ones((N, N, N, 3))
+        ww_array[..., 0] = lower
+        ww_array[..., 1] = target
+        ww_array[..., 2] = upper
+
+    mcdc.simulation.weight_windows(mesh, ww_array)
+
+    mcdc_container, data = preparation()
+    return mcdc_container[0], data
+
+
+def make_ww_model_distinct():
+    import mcdc
+
+    mesh, N = make_base_model()
+
+    ww_array = np.empty((N, N, N, 3))
+
+    for i in range(N):
+        for j in range(N):
+            for k in range(N):
+                val = 100*i + 10*j + k + 1
+                ww_array[i, j, k, 0] = val
+                ww_array[i, j, k, 1] = 1000 + val
+                ww_array[i, j, k, 2] = 2000 + val
+
+    mcdc.simulation.weight_windows(mesh, ww_array)
+
+    mcdc_container, data = preparation()
+    return mcdc_container[0], data
 
 
 # =========================================================================== #
@@ -101,7 +131,7 @@ def make_ww_model(lower=0.1, target=1.0, upper=1.0, mess_up_size=False):
 )
 def test_error_throw(capsys, kwargs, expected_msg):
   with pytest.raises(SystemExit):
-    make_ww_model(**kwargs)
+    make_ww_model_params(**kwargs)
 
   out = capsys.readouterr().out
   assert expected_msg in out
@@ -130,7 +160,7 @@ def test_split_from_weight_window():
   init_weight = 2.0 + TINY
   particles[0]["w"] = init_weight
   threshold = 1.0
-  mcdc_obj, data = make_ww_model()
+  mcdc_obj, data = make_ww_model_distinct()
 
   # get bank and init size
   bank = mcdc_obj["bank_active"]
@@ -150,3 +180,36 @@ def test_split_from_weight_window():
     pnew = bank["particles"][init_bank_size + i]
     assert pnew["w"] == p1["w"]
 
+
+def test_query_weight_window():
+  p = np.zeros(1, type_.particle_data)
+
+  mcdc_obj, data = make_ww_model_distinct()
+  ww_obj = mcdc_obj["weight_windows"]
+
+  mesh = mcdc_obj["meshes"][ww_obj["mesh_ID"]]
+  # hardcode mesh params
+  nx, ny, nz = 3, 3, 3
+  xmin, ymin, zmin = -1.0, -1.0, 0.0
+  dx, dy, dz = 2/3, 2/3, 10/3
+
+
+  # loop over all bins, check query against expected ww
+  for ix in range(nx):
+    for iy in range(ny):
+      for iz in range(nz):
+        p[0]["x"] = xmin + dx * (ix + 0.5)
+        p[0]["y"] = ymin + dy * (iy + 0.5)
+        p[0]["z"] = zmin + dz * (iz + 0.5)
+
+        lower, target, upper = query_weight_window(
+          p, mcdc_obj, data
+        )
+        exp_lower = 100*ix + 10*iy + iz + 1
+        exp_target = 1000 + exp_lower
+        exp_upper = 2000 + exp_lower
+
+        assert lower == exp_lower
+        assert target == exp_target
+        assert upper == exp_upper
+

From 8c2e62661fb5b6acff7992703d8ddb5e3723f845 Mon Sep 17 00:00:00 2001
From: Nathan Glaser <glaserna@oregonstate.edu>
Date: Sun, 12 Apr 2026 21:42:36 -0700
Subject: [PATCH 13/53] Clean up tests a bit, removed needless setup

---
 .../transport/technique/weight_windows.py     | 41 ++++---------------
 1 file changed, 9 insertions(+), 32 deletions(-)

diff --git a/test/unit/transport/technique/weight_windows.py b/test/unit/transport/technique/weight_windows.py
index 7974d6910..0687c26e0 100644
--- a/test/unit/transport/technique/weight_windows.py
+++ b/test/unit/transport/technique/weight_windows.py
@@ -21,34 +21,12 @@
 # Helper method for creating a dummy model
 # =========================================================================== #
 
-def make_base_model():
-
-    # Geometry
-    cylinder = mcdc.Surface.CylinderZ(radius=0.5)
+def make_mesh():
+    # constants
     pitch = 2.0
     height = 10.0
-
-    x0 = mcdc.Surface.PlaneX(x=-pitch / 2, boundary_condition="reflective")
-    x1 = mcdc.Surface.PlaneX(x=pitch / 2, boundary_condition="reflective")
-    y0 = mcdc.Surface.PlaneY(y=-pitch / 2, boundary_condition="reflective")
-    y1 = mcdc.Surface.PlaneY(y=pitch / 2, boundary_condition="reflective")
-    z0 = mcdc.Surface.PlaneZ(z=0.0, boundary_condition="reflective")
-    z1 = mcdc.Surface.PlaneZ(z=height, boundary_condition="reflective")
-
-    mcdc.Cell(-cylinder)
-    mcdc.Cell(+x0 & -x1 & +y0 & -y1 & +z0 & -z1 & +cylinder)
-
-    # Source
-    mcdc.Source(position=[0.0, 0.0, 0.0], isotropic=True, time=0.0, energy=14.1e6)
-
-    # Settings
-    mcdc.settings.N_particle = 20
-    mcdc.settings.N_batch = 2
-    mcdc.settings.time_boundary = 1.0
-    mcdc.settings.active_bank_buffer = 1000
-
-    # Mesh
     N = 3
+
     mesh = mcdc.MeshUniform(
         "mesh",
         x=(-pitch/2, pitch/N, N),
@@ -62,7 +40,7 @@ def make_base_model():
 def make_ww_model_params(lower=0.1, target=1.0, upper=1.0, mess_up_size=False):
     import mcdc
 
-    mesh, N = make_base_model()
+    mesh, N = make_mesh()
 
     if mess_up_size:
         ww_array = np.ones((N, N, 4, 3))
@@ -81,7 +59,7 @@ def make_ww_model_params(lower=0.1, target=1.0, upper=1.0, mess_up_size=False):
 def make_ww_model_distinct():
     import mcdc
 
-    mesh, N = make_base_model()
+    mesh, N = make_mesh()
 
     ww_array = np.empty((N, N, N, 3))
 
@@ -185,13 +163,12 @@ def test_query_weight_window():
   p = np.zeros(1, type_.particle_data)
 
   mcdc_obj, data = make_ww_model_distinct()
-  ww_obj = mcdc_obj["weight_windows"]
 
-  mesh = mcdc_obj["meshes"][ww_obj["mesh_ID"]]
   # hardcode mesh params
-  nx, ny, nz = 3, 3, 3
-  xmin, ymin, zmin = -1.0, -1.0, 0.0
-  dx, dy, dz = 2/3, 2/3, 10/3
+  pitch, height, N = 2.0, 10.0, 3
+  nx, ny, nz = N, N, N
+  xmin, ymin, zmin = -pitch/2, -pitch/2, 0.0
+  dx, dy, dz = pitch/N, pitch/N, height/N
 
 
   # loop over all bins, check query against expected ww

From 4933873c499cc5d7f8e0e9b94dd55f96a00e550b Mon Sep 17 00:00:00 2001
From: Nathan Glaser <glaserna@oregonstate.edu>
Date: Sun, 12 Apr 2026 21:49:32 -0700
Subject: [PATCH 14/53] integrate the changed names into ww tests

---
 test/unit/transport/technique/weight_windows.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/unit/transport/technique/weight_windows.py b/test/unit/transport/technique/weight_windows.py
index 0687c26e0..32d5c0712 100644
--- a/test/unit/transport/technique/weight_windows.py
+++ b/test/unit/transport/technique/weight_windows.py
@@ -155,7 +155,7 @@ def test_split_from_weight_window():
   assert p1["w"] == init_weight / num_split
   assert particle_bank_module.get_bank_size(bank) == init_bank_size + num_new
   for i in range(2):
-    pnew = bank["particles"][init_bank_size + i]
+    pnew = bank["particle_data"][init_bank_size + i]
     assert pnew["w"] == p1["w"]
 
 

From 6e8cb5b4bfafab9b3c566f14da66624329afdabb Mon Sep 17 00:00:00 2001
From: Nathan Glaser <glaserna@oregonstate.edu>
Date: Sun, 12 Apr 2026 21:51:23 -0700
Subject: [PATCH 15/53] black format

---
 mcdc/object_/technique.py                     |   9 +-
 mcdc/transport/technique.py                   |   3 +-
 .../transport/technique/weight_windows.py     | 147 +++++++++---------
 3 files changed, 78 insertions(+), 81 deletions(-)

diff --git a/mcdc/object_/technique.py b/mcdc/object_/technique.py
index 4017f3b73..5509774ed 100644
--- a/mcdc/object_/technique.py
+++ b/mcdc/object_/technique.py
@@ -97,12 +97,16 @@ def __call__(self, mesh, weight_windows):
             case "uniform_mesh":
                 nx, ny, nz = mesh.Nx, mesh.Ny, mesh.Nz
             case "structured_mesh":
-                nx, ny, nz = mesh.x.shape[0]-1, mesh.y.shape[0]-1, mesh.z.shape[0]-1
+                nx, ny, nz = (
+                    mesh.x.shape[0] - 1,
+                    mesh.y.shape[0] - 1,
+                    mesh.z.shape[0] - 1,
+                )
             case _:
                 print_error(
                     f"{type(mesh).__name__} is not supported for weight windows"
                 )
-        
+
         # check correct shape
         mesh_shape = (nx, ny, nz)
         ww_shape = weight_windows.shape
@@ -134,7 +138,6 @@ def __call__(self, mesh, weight_windows):
             )
 
 
-
 # ======================================================================================
 # Population control
 # ======================================================================================
diff --git a/mcdc/transport/technique.py b/mcdc/transport/technique.py
index b14d25b01..5ca54c1a2 100644
--- a/mcdc/transport/technique.py
+++ b/mcdc/transport/technique.py
@@ -12,7 +12,6 @@
 import mcdc.transport.rng as rng
 import mcdc.transport.util as util
 
-
 # ======================================================================================
 # Weight Roulette
 # ======================================================================================
@@ -50,7 +49,7 @@ def query_weight_window(particle_container, mcdc, data):
     index = ((idx * ww_obj["Ny"]) + idy) * ww_obj["Nz"] + idz
     lower = ww_get.lower_weights(index, ww_obj, data)
     target = ww_get.target_weights(index, ww_obj, data)
-    upper = ww_get.upper_weights(index, ww_obj, data) 
+    upper = ww_get.upper_weights(index, ww_obj, data)
     return lower, target, upper
 
 
diff --git a/test/unit/transport/technique/weight_windows.py b/test/unit/transport/technique/weight_windows.py
index 32d5c0712..b90b6a257 100644
--- a/test/unit/transport/technique/weight_windows.py
+++ b/test/unit/transport/technique/weight_windows.py
@@ -6,21 +6,20 @@
 from mcdc.main import preparation
 import mcdc.numba_types as type_
 from mcdc.transport.technique import (
-  roulette_from_weight_window,
-  split_from_weight_window,
-  query_weight_window,
-  weight_windows,
-  particle_bank_module
+    roulette_from_weight_window,
+    split_from_weight_window,
+    query_weight_window,
+    weight_windows,
+    particle_bank_module,
 )
 from mcdc.transport.mesh.interface import get_indices
-from mcdc.constant import (
-  TINY
-)
+from mcdc.constant import TINY
 
 # =========================================================================== #
 # Helper method for creating a dummy model
 # =========================================================================== #
 
+
 def make_mesh():
     # constants
     pitch = 2.0
@@ -29,9 +28,9 @@ def make_mesh():
 
     mesh = mcdc.MeshUniform(
         "mesh",
-        x=(-pitch/2, pitch/N, N),
-        y=(-pitch/2, pitch/N, N),
-        z=(0.0, height/N, N)
+        x=(-pitch / 2, pitch / N, N),
+        y=(-pitch / 2, pitch / N, N),
+        z=(0.0, height / N, N),
     )
 
     return mesh, N
@@ -66,7 +65,7 @@ def make_ww_model_distinct():
     for i in range(N):
         for j in range(N):
             for k in range(N):
-                val = 100*i + 10*j + k + 1
+                val = 100 * i + 10 * j + k + 1
                 ww_array[i, j, k, 0] = val
                 ww_array[i, j, k, 1] = 1000 + val
                 ww_array[i, j, k, 2] = 2000 + val
@@ -108,85 +107,81 @@ def make_ww_model_distinct():
     ],
 )
 def test_error_throw(capsys, kwargs, expected_msg):
-  with pytest.raises(SystemExit):
-    make_ww_model_params(**kwargs)
+    with pytest.raises(SystemExit):
+        make_ww_model_params(**kwargs)
 
-  out = capsys.readouterr().out
-  assert expected_msg in out
+    out = capsys.readouterr().out
+    assert expected_msg in out
 
 
 # =========================================================================== #
-# Tests for helper methods 
+# Tests for helper methods
 # =========================================================================== #
 
 
 def test_roullete_from_weight_window():
-  # because of rng, want to loop over to hit both branches
-  for i in range(10):
-    particles = np.zeros(1, type_.particle_data)
-    particles[0]["w"] = 0.1
-    target = 0.2
-    threshold = 0.1 + TINY
-    
-    roulette_from_weight_window(particles, threshold, target)
-    p = particles[0]
-    assert (p["w"] == target or not p["alive"])
+    # because of rng, want to loop over to hit both branches
+    for i in range(10):
+        particles = np.zeros(1, type_.particle_data)
+        particles[0]["w"] = 0.1
+        target = 0.2
+        threshold = 0.1 + TINY
+
+        roulette_from_weight_window(particles, threshold, target)
+        p = particles[0]
+        assert p["w"] == target or not p["alive"]
 
 
 def test_split_from_weight_window():
-  particles = np.zeros(1, type_.particle_data)
-  init_weight = 2.0 + TINY
-  particles[0]["w"] = init_weight
-  threshold = 1.0
-  mcdc_obj, data = make_ww_model_distinct()
+    particles = np.zeros(1, type_.particle_data)
+    init_weight = 2.0 + TINY
+    particles[0]["w"] = init_weight
+    threshold = 1.0
+    mcdc_obj, data = make_ww_model_distinct()
 
-  # get bank and init size
-  bank = mcdc_obj["bank_active"]
-  init_bank_size = particle_bank_module.get_bank_size(bank)
+    # get bank and init size
+    bank = mcdc_obj["bank_active"]
+    init_bank_size = particle_bank_module.get_bank_size(bank)
 
-  # split
-  split_from_weight_window(particles, threshold, mcdc_obj)
+    # split
+    split_from_weight_window(particles, threshold, mcdc_obj)
 
-  p1 = particles[0]
-  num_split = np.ceil(init_weight / threshold)
-  num_new = num_split - 1
+    p1 = particles[0]
+    num_split = np.ceil(init_weight / threshold)
+    num_new = num_split - 1
 
-  # check weight of original particle
-  assert p1["w"] == init_weight / num_split
-  assert particle_bank_module.get_bank_size(bank) == init_bank_size + num_new
-  for i in range(2):
-    pnew = bank["particle_data"][init_bank_size + i]
-    assert pnew["w"] == p1["w"]
+    # check weight of original particle
+    assert p1["w"] == init_weight / num_split
+    assert particle_bank_module.get_bank_size(bank) == init_bank_size + num_new
+    for i in range(2):
+        pnew = bank["particle_data"][init_bank_size + i]
+        assert pnew["w"] == p1["w"]
 
 
 def test_query_weight_window():
-  p = np.zeros(1, type_.particle_data)
-
-  mcdc_obj, data = make_ww_model_distinct()
-
-  # hardcode mesh params
-  pitch, height, N = 2.0, 10.0, 3
-  nx, ny, nz = N, N, N
-  xmin, ymin, zmin = -pitch/2, -pitch/2, 0.0
-  dx, dy, dz = pitch/N, pitch/N, height/N
-
-
-  # loop over all bins, check query against expected ww
-  for ix in range(nx):
-    for iy in range(ny):
-      for iz in range(nz):
-        p[0]["x"] = xmin + dx * (ix + 0.5)
-        p[0]["y"] = ymin + dy * (iy + 0.5)
-        p[0]["z"] = zmin + dz * (iz + 0.5)
-
-        lower, target, upper = query_weight_window(
-          p, mcdc_obj, data
-        )
-        exp_lower = 100*ix + 10*iy + iz + 1
-        exp_target = 1000 + exp_lower
-        exp_upper = 2000 + exp_lower
-
-        assert lower == exp_lower
-        assert target == exp_target
-        assert upper == exp_upper
-
+    p = np.zeros(1, type_.particle_data)
+
+    mcdc_obj, data = make_ww_model_distinct()
+
+    # hardcode mesh params
+    pitch, height, N = 2.0, 10.0, 3
+    nx, ny, nz = N, N, N
+    xmin, ymin, zmin = -pitch / 2, -pitch / 2, 0.0
+    dx, dy, dz = pitch / N, pitch / N, height / N
+
+    # loop over all bins, check query against expected ww
+    for ix in range(nx):
+        for iy in range(ny):
+            for iz in range(nz):
+                p[0]["x"] = xmin + dx * (ix + 0.5)
+                p[0]["y"] = ymin + dy * (iy + 0.5)
+                p[0]["z"] = zmin + dz * (iz + 0.5)
+
+                lower, target, upper = query_weight_window(p, mcdc_obj, data)
+                exp_lower = 100 * ix + 10 * iy + iz + 1
+                exp_target = 1000 + exp_lower
+                exp_upper = 2000 + exp_lower
+
+                assert lower == exp_lower
+                assert target == exp_target
+                assert upper == exp_upper

From a57e72e6d70ed74dc06c4d51c1ef8a46c700b813 Mon Sep 17 00:00:00 2001
From: Nathan Glaser <glaserna@oregonstate.edu>
Date: Sun, 12 Apr 2026 22:15:05 -0700
Subject: [PATCH 16/53] add oneliner on how to add weight windows to a
 simulation

---
 docs/source/pythonapi/index.rst | 1 +
 1 file changed, 1 insertion(+)

diff --git a/docs/source/pythonapi/index.rst b/docs/source/pythonapi/index.rst
index 3b108838a..c16a2fa5f 100644
--- a/docs/source/pythonapi/index.rst
+++ b/docs/source/pythonapi/index.rst
@@ -89,6 +89,7 @@ Techniques are enabled by calling methods on the ``mcdc.simulation`` singleton:
 - ``mcdc.simulation.implicit_capture(active=True)``
 - ``mcdc.simulation.weighted_emission(active=True, weight_target=1.0)``
 - ``mcdc.simulation.weight_roulette(weight_threshold=0.0, weight_target=1.0)``
+- ``mcdc.simulation.weight_windows(mesh, weight_windows)``
 - ``mcdc.simulation.population_control(active=True)``
 
 Running

From 4bc9c774fbeb43aad0b0cb5e1992b0f6b308caf0 Mon Sep 17 00:00:00 2001
From: Nathan Glaser <glaserna@oregonstate.edu>
Date: Sun, 12 Apr 2026 22:25:29 -0700
Subject: [PATCH 17/53] add docstrings and some comments for clarity

---
 mcdc/transport/technique.py | 83 +++++++++++++++++++++++++++++++++----
 1 file changed, 74 insertions(+), 9 deletions(-)

diff --git a/mcdc/transport/technique.py b/mcdc/transport/technique.py
index 5ca54c1a2..a9c72fac3 100644
--- a/mcdc/transport/technique.py
+++ b/mcdc/transport/technique.py
@@ -35,18 +35,55 @@ def weight_roulette(particle_container, simulation):
 
 
 @njit
-def weight_windows(particle_container, mcdc, data):
-    [lower, target, upper] = query_weight_window(particle_container, mcdc, data)
-    split_from_weight_window(particle_container, upper, mcdc)
+def weight_windows(particle_container, simulation, data):
+    """
+    Apply weight window splitting and rouletting to a particle.
+
+    Parameters
+    ----------
+    particle_container : ndarray
+        Container holding the particle.
+    simulation : object
+        Simulation state containing weight window and mesh data.
+    data : object
+        Simulation data for array access.
+    """
+    [lower, target, upper] = query_weight_window(particle_container, simulation, data)
+    split_from_weight_window(particle_container, upper, simulation)
     roulette_from_weight_window(particle_container, lower, target)
 
 
 @njit
-def query_weight_window(particle_container, mcdc, data):
-    ww_obj = mcdc["weight_windows"]
-    mesh = mcdc["meshes"][ww_obj["mesh_ID"]]
-    idx, idy, idz = get_mesh_indices(particle_container, mesh, mcdc, data)
+def query_weight_window(particle_container, simulation, data):
+    """
+    Query weight window bounds for the particle.
+
+    Parameters
+    ----------
+    particle_container : ndarray
+        Container holding the particle.
+    simulation : object
+        Simulation state containing weight window and mesh data.
+    data : object
+        Simulation data for array access.
+
+    Returns
+    -------
+    lower : float
+        Lower weight bound.
+    target : float
+        Target weight.
+    upper : float
+        Upper weight bound.
+    """
+    # grab objects
+    ww_obj = simulation["weight_windows"]
+    mesh = simulation["meshes"][ww_obj["mesh_ID"]]
+    # find spatial indices from mesh
+    idx, idy, idz = get_mesh_indices(particle_container, mesh, simulation, data)
+    # get index of ww in ww arrays
     index = ((idx * ww_obj["Ny"]) + idy) * ww_obj["Nz"] + idz
+    # grab the actual ww parameters
     lower = ww_get.lower_weights(index, ww_obj, data)
     target = ww_get.target_weights(index, ww_obj, data)
     upper = ww_get.upper_weights(index, ww_obj, data)
@@ -54,21 +91,49 @@ def query_weight_window(particle_container, mcdc, data):
 
 
 @njit
-def split_from_weight_window(particle_container, threshold_weight, mcdc):
+def split_from_weight_window(particle_container, threshold_weight, simulation):
+    """
+    Split a particle if its weight exceeds the threshold.
+
+    Parameters
+    ----------
+    particle_container : ndarray
+        Container holding the particle.
+    threshold_weight : float
+        Upper weight bound triggering splitting.
+    simulation : object
+        Simulation state used for banking split particles.
+    """
     particle = particle_container[0]
     weight = particle["w"]
     if weight > threshold_weight:
+        # determine how many to split into
         num_split = math.ceil(weight / threshold_weight)
+        # distribute weight
         particle["w"] = weight / num_split
         for _ in range(num_split - 1):
-            particle_bank_module.bank_active_particle(particle_container, mcdc)
+            # bank split particles into the active bank
+            particle_bank_module.bank_active_particle(particle_container, simulation)
 
 
 @njit
 def roulette_from_weight_window(particle_container, w_threshold, w_target):
+    """
+    Russian roulette particle if weight is below threshold.
+
+    Parameters
+    ----------
+    particle_container : ndarray
+        Container holding the particle.
+    w_threshold : float
+        Lower weight bound triggering roulette.
+    w_target : float
+        Target weight assigned upon survival.
+    """
     particle = particle_container[0]
     if particle["w"] < w_threshold:
         survival_probability = particle["w"] / w_target
+        # sample random number to determine survival
         if rng.lcg(particle_container) < survival_probability:
             particle["w"] = w_target
         else:

From 33438fe261bf1350e42529574821fee009b023ca Mon Sep 17 00:00:00 2001
From: Nathan Glaser <glaserna@oregonstate.edu>
Date: Sun, 12 Apr 2026 23:33:56 -0700
Subject: [PATCH 18/53] clean up imports, a few more comments for clarity

---
 mcdc/object_/technique.py                       | 4 +---
 mcdc/transport/technique.py                     | 1 +
 test/unit/transport/technique/weight_windows.py | 7 ++++++-
 3 files changed, 8 insertions(+), 4 deletions(-)

diff --git a/mcdc/object_/technique.py b/mcdc/object_/technique.py
index 5509774ed..f18fcafd3 100644
--- a/mcdc/object_/technique.py
+++ b/mcdc/object_/technique.py
@@ -1,8 +1,6 @@
-from dataclasses import dataclass
 from numpy.typing import NDArray
 from numpy import float64
-import numpy as np
-from mcdc.object_.base import ObjectBase, ObjectSingleton
+from mcdc.object_.base import ObjectSingleton
 from mcdc.object_.mesh import MeshBase
 from mcdc.print_ import print_error
 
diff --git a/mcdc/transport/technique.py b/mcdc/transport/technique.py
index a9c72fac3..cb586e694 100644
--- a/mcdc/transport/technique.py
+++ b/mcdc/transport/technique.py
@@ -4,6 +4,7 @@
 from numba import njit
 
 ####
+
 import mcdc.mcdc_get.weight_windows as ww_get
 import mcdc.numba_types as type_
 from mcdc.transport.mesh import get_indices as get_mesh_indices
diff --git a/test/unit/transport/technique/weight_windows.py b/test/unit/transport/technique/weight_windows.py
index b90b6a257..7a12eb5d0 100644
--- a/test/unit/transport/technique/weight_windows.py
+++ b/test/unit/transport/technique/weight_windows.py
@@ -44,6 +44,7 @@ def make_ww_model_params(lower=0.1, target=1.0, upper=1.0, mess_up_size=False):
     if mess_up_size:
         ww_array = np.ones((N, N, 4, 3))
     else:
+        # global assign for simplicity
         ww_array = np.ones((N, N, N, 3))
         ww_array[..., 0] = lower
         ww_array[..., 1] = target
@@ -62,6 +63,7 @@ def make_ww_model_distinct():
 
     ww_array = np.empty((N, N, N, 3))
 
+    # value at index is related to index, easy to predict during later test
     for i in range(N):
         for j in range(N):
             for k in range(N):
@@ -121,7 +123,7 @@ def test_error_throw(capsys, kwargs, expected_msg):
 
 def test_roullete_from_weight_window():
     # because of rng, want to loop over to hit both branches
-    for i in range(10):
+    for _ in range(10):
         particles = np.zeros(1, type_.particle_data)
         particles[0]["w"] = 0.1
         target = 0.2
@@ -152,6 +154,7 @@ def test_split_from_weight_window():
 
     # check weight of original particle
     assert p1["w"] == init_weight / num_split
+    # check particles were created and check their weights
     assert particle_bank_module.get_bank_size(bank) == init_bank_size + num_new
     for i in range(2):
         pnew = bank["particle_data"][init_bank_size + i]
@@ -173,10 +176,12 @@ def test_query_weight_window():
     for ix in range(nx):
         for iy in range(ny):
             for iz in range(nz):
+                # put particle in center of current mesh bin
                 p[0]["x"] = xmin + dx * (ix + 0.5)
                 p[0]["y"] = ymin + dy * (iy + 0.5)
                 p[0]["z"] = zmin + dz * (iz + 0.5)
 
+                # query and predict
                 lower, target, upper = query_weight_window(p, mcdc_obj, data)
                 exp_lower = 100 * ix + 10 * iy + iz + 1
                 exp_target = 1000 + exp_lower

From c01de007ec06e2406a3c4a70768fc5ae8a04aa59 Mon Sep 17 00:00:00 2001
From: Nathan Glaser <glaserna@oregonstate.edu>
Date: Tue, 14 Apr 2026 00:07:48 -0700
Subject: [PATCH 19/53] add ww numba type

---
 mcdc/numba_types.py | 15 +++++++++++++++
 1 file changed, 15 insertions(+)

diff --git a/mcdc/numba_types.py b/mcdc/numba_types.py
index 8ca7ba4f4..46aa629d6 100644
--- a/mcdc/numba_types.py
+++ b/mcdc/numba_types.py
@@ -574,6 +574,20 @@
     ('weight_target', float64),
 ])
 
+weight_windows = into_dtype([
+    ('active', bool),
+    ('mesh_ID', int64),
+    ('Nx', int64),
+    ('Ny', int64),
+    ('Nz', int64),
+    ('lower_weights_offset', int64),
+    ('lower_weights_length', int64),
+    ('target_weights_offset', int64),
+    ('target_weights_length', int64),
+    ('upper_weights_offset', int64),
+    ('upper_weights_length', int64),
+])
+
 weighted_emission = into_dtype([
     ('active', bool),
     ('weight_target', float64),
@@ -815,6 +829,7 @@ def set_simulation(N: dict):
         ('implicit_capture', implicit_capture),
         ('weighted_emission', weighted_emission),
         ('weight_roulette', weight_roulette),
+        ('weight_windows', weight_windows),
         ('population_control', population_control),
         ('gpu_meta', gpu_meta),
         ('bank_future', bank_future),

From b8a769bf46b6ec1f0b8db83d3967c60bc9742f73 Mon Sep 17 00:00:00 2001
From: Nathan Glaser <glaserna@oregonstate.edu>
Date: Tue, 14 Apr 2026 01:07:47 -0700
Subject: [PATCH 20/53] properly initialize weight_windows singleton

---
 mcdc/object_/technique.py | 9 +++++++--
 1 file changed, 7 insertions(+), 2 deletions(-)

diff --git a/mcdc/object_/technique.py b/mcdc/object_/technique.py
index f18fcafd3..6e93a845e 100644
--- a/mcdc/object_/technique.py
+++ b/mcdc/object_/technique.py
@@ -1,7 +1,7 @@
 from numpy.typing import NDArray
-from numpy import float64
+from numpy import float64, array
 from mcdc.object_.base import ObjectSingleton
-from mcdc.object_.mesh import MeshBase
+from mcdc.object_.mesh import MeshBase, MeshUniform
 from mcdc.print_ import print_error
 
 # ======================================================================================
@@ -88,6 +88,11 @@ class WeightWindows(ObjectSingleton):
 
     def __init__(self):
         self.active = False
+        self.mesh_ID = -1  # skirt around having to create a MeshBase instance
+        self.Nx, self.Ny, self.Nz = 1, 1, 1
+        self.lower_weights = array([1.0])
+        self.target_weights = array([1.0])
+        self.upper_weights = array([1.0])
 
     def __call__(self, mesh, weight_windows):
         # get mesh size

From e7aa22a367ed03835fa832aeff51f88431164d55 Mon Sep 17 00:00:00 2001
From: Nathan Glaser <glaserna@oregonstate.edu>
Date: Tue, 14 Apr 2026 11:42:56 -0700
Subject: [PATCH 21/53] add energy bins and indexing to technique for ww

---
 mcdc/object_/technique.py   | 39 ++++++++++++++++++++++++++++---------
 mcdc/transport/technique.py | 24 ++++++++++++++++++-----
 2 files changed, 49 insertions(+), 14 deletions(-)

diff --git a/mcdc/object_/technique.py b/mcdc/object_/technique.py
index 6e93a845e..381b031e9 100644
--- a/mcdc/object_/technique.py
+++ b/mcdc/object_/technique.py
@@ -1,5 +1,6 @@
 from numpy.typing import NDArray
-from numpy import float64, array
+import numpy as np
+from mcdc.constant import INF
 from mcdc.object_.base import ObjectSingleton
 from mcdc.object_.mesh import MeshBase, MeshUniform
 from mcdc.print_ import print_error
@@ -76,25 +77,32 @@ class WeightWindows(ObjectSingleton):
     label: str = "weight_windows"
 
     active: bool
+
+    # energy
+    energy_bounds: NDArray[np.float64]
+    Ne: int
+    # space
     mesh: MeshBase
     Nx: int
     Ny: int
     Nz: int
 
     # flattened arrays of ww params
-    lower_weights: NDArray[float64]
-    target_weights: NDArray[float64]
-    upper_weights: NDArray[float64]
+    lower_weights: NDArray[np.float64]
+    target_weights: NDArray[np.float64]
+    upper_weights: NDArray[np.float64]
 
     def __init__(self):
         self.active = False
+        self.energy_bounds = np.array([0.0, 1.0])
+        self.Ne = 1
         self.mesh_ID = -1  # skirt around having to create a MeshBase instance
         self.Nx, self.Ny, self.Nz = 1, 1, 1
-        self.lower_weights = array([1.0])
-        self.target_weights = array([1.0])
-        self.upper_weights = array([1.0])
+        self.lower_weights = np.array([1.0])
+        self.target_weights = np.array([1.0])
+        self.upper_weights = np.array([1.0])
 
-    def __call__(self, mesh, weight_windows):
+    def __call__(self, mesh, weight_windows, energy=np.array([-0.5, INF])):
         # get mesh size
         match mesh.label:
             case "uniform_mesh":
@@ -109,17 +117,30 @@ def __call__(self, mesh, weight_windows):
                 print_error(
                     f"{type(mesh).__name__} is not supported for weight windows"
                 )
+        # validate energy as strictly increasing
+        if not (np.diff(energy) > 0).all():
+            print_error(
+                "Energy bounds must be strictly increasing"
+            )
+        # get energy size
+        if len(energy.shape) != 1:
+            print_error(
+                f"Invalid shape for energy; expected 1D got {len(energy.shape)}D"
+            )
+        ne = energy.shape[0] - 1
 
         # check correct shape
         mesh_shape = (nx, ny, nz)
         ww_shape = weight_windows.shape
-        expected_shape = (*mesh_shape, 3)
+        expected_shape = (ne, *mesh_shape, 3)
         if ww_shape != expected_shape:
             print_error(
                 f"Weight window array has shape {ww_shape}, but expected {expected_shape}"
             )
 
         self.active = True
+        self.energy_bounds = energy
+        self.Ne = ne
         self.mesh = mesh
         self.Nx, self.Ny, self.Nz = mesh_shape
         self.lower_weights = weight_windows[..., 0].reshape(-1)
diff --git a/mcdc/transport/technique.py b/mcdc/transport/technique.py
index cb586e694..9798744e4 100644
--- a/mcdc/transport/technique.py
+++ b/mcdc/transport/technique.py
@@ -79,17 +79,31 @@ def query_weight_window(particle_container, simulation, data):
     """
     # grab objects
     ww_obj = simulation["weight_windows"]
-    mesh = simulation["meshes"][ww_obj["mesh_ID"]]
-    # find spatial indices from mesh
-    idx, idy, idz = get_mesh_indices(particle_container, mesh, simulation, data)
-    # get index of ww in ww arrays
-    index = ((idx * ww_obj["Ny"]) + idy) * ww_obj["Nz"] + idz
+    index = get_ww_index(particle_container, ww_obj, simulation, data) 
     # grab the actual ww parameters
     lower = ww_get.lower_weights(index, ww_obj, data)
     target = ww_get.target_weights(index, ww_obj, data)
     upper = ww_get.upper_weights(index, ww_obj, data)
     return lower, target, upper
 
+@njit
+def get_ww_index(particle_container, ww_obj, simulation, data):
+    particle = particle_container[0]
+
+    # get energy index
+    energy_bounds = ww_obj["energy_bounds"]
+    if simulation["settings"]["neutron_multigroup_mode"]:
+        energy = particle["g"]
+    else:
+        energy = particle["E"]
+    ie = util.find_bin(energy, energy_bounds)
+
+    # get spatial index
+    mesh = simulation["meshes"][ww_obj["mesh_ID"]]
+    idx, idy, idz = get_mesh_indices(particle_container, mesh, simulation, data)
+
+    index = (((ie * ww_obj["Nx"]) + idx) * ww_obj["Ny"] + idy) * ww_obj["Nz"] + idz
+    return index
 
 @njit
 def split_from_weight_window(particle_container, threshold_weight, simulation):

From 38e4e213d87714d0d458487b9de09fe83eea710f Mon Sep 17 00:00:00 2001
From: Nathan Glaser <glaserna@oregonstate.edu>
Date: Tue, 14 Apr 2026 12:04:16 -0700
Subject: [PATCH 22/53] correct unit tests for ww

---
 mcdc/transport/technique.py                   |  2 +-
 .../transport/technique/weight_windows.py     | 69 +++++++++++--------
 2 files changed, 41 insertions(+), 30 deletions(-)

diff --git a/mcdc/transport/technique.py b/mcdc/transport/technique.py
index 9798744e4..0a4edb80c 100644
--- a/mcdc/transport/technique.py
+++ b/mcdc/transport/technique.py
@@ -91,7 +91,7 @@ def get_ww_index(particle_container, ww_obj, simulation, data):
     particle = particle_container[0]
 
     # get energy index
-    energy_bounds = ww_obj["energy_bounds"]
+    energy_bounds = ww_get.energy_bounds_all(ww_obj, data) 
     if simulation["settings"]["neutron_multigroup_mode"]:
         energy = particle["g"]
     else:
diff --git a/test/unit/transport/technique/weight_windows.py b/test/unit/transport/technique/weight_windows.py
index 7a12eb5d0..9991c37c9 100644
--- a/test/unit/transport/technique/weight_windows.py
+++ b/test/unit/transport/technique/weight_windows.py
@@ -40,12 +40,13 @@ def make_ww_model_params(lower=0.1, target=1.0, upper=1.0, mess_up_size=False):
     import mcdc
 
     mesh, N = make_mesh()
+    Ne = 1
 
     if mess_up_size:
-        ww_array = np.ones((N, N, 4, 3))
+        ww_array = np.ones((Ne, N, N, 4, 3))
     else:
         # global assign for simplicity
-        ww_array = np.ones((N, N, N, 3))
+        ww_array = np.ones((Ne, N, N, N, 3))
         ww_array[..., 0] = lower
         ww_array[..., 1] = target
         ww_array[..., 2] = upper
@@ -60,19 +61,22 @@ def make_ww_model_distinct():
     import mcdc
 
     mesh, N = make_mesh()
+    energy = np.linspace(0.0, 6.0, 7)
+    Ne = 6
 
-    ww_array = np.empty((N, N, N, 3))
+    ww_array = np.empty((Ne, N, N, N, 3))
 
     # value at index is related to index, easy to predict during later test
-    for i in range(N):
-        for j in range(N):
-            for k in range(N):
-                val = 100 * i + 10 * j + k + 1
-                ww_array[i, j, k, 0] = val
-                ww_array[i, j, k, 1] = 1000 + val
-                ww_array[i, j, k, 2] = 2000 + val
+    for e in range(Ne):
+        for i in range(N):
+            for j in range(N):
+                for k in range(N):
+                    val = 1000 * e + 100 * i + 10 * j + k + 1
+                    ww_array[e, i, j, k, 0] = val
+                    ww_array[e, i, j, k, 1] = 10000 + val
+                    ww_array[e, i, j, k, 2] = 20000 + val
 
-    mcdc.simulation.weight_windows(mesh, ww_array)
+    mcdc.simulation.weight_windows(mesh, ww_array, energy)
 
     mcdc_container, data = preparation()
     return mcdc_container[0], data
@@ -89,7 +93,7 @@ def make_ww_model_distinct():
         # incorrect size
         (
             {"mess_up_size": True},
-            "Weight window array has shape (3, 3, 4, 3), but expected (3, 3, 3, 3)",
+            "Weight window array has shape (1, 3, 3, 4, 3), but expected (1, 3, 3, 3, 3)",
         ),
         # negative lower
         (
@@ -165,6 +169,7 @@ def test_query_weight_window():
     p = np.zeros(1, type_.particle_data)
 
     mcdc_obj, data = make_ww_model_distinct()
+    mcdc_obj["settings"]["neutron_multigroup_mode"] = False
 
     # hardcode mesh params
     pitch, height, N = 2.0, 10.0, 3
@@ -172,21 +177,27 @@ def test_query_weight_window():
     xmin, ymin, zmin = -pitch / 2, -pitch / 2, 0.0
     dx, dy, dz = pitch / N, pitch / N, height / N
 
+    # hardcode energy params
+    energies = np.linspace(0.5, 5.5, 6)
+
     # loop over all bins, check query against expected ww
-    for ix in range(nx):
-        for iy in range(ny):
-            for iz in range(nz):
-                # put particle in center of current mesh bin
-                p[0]["x"] = xmin + dx * (ix + 0.5)
-                p[0]["y"] = ymin + dy * (iy + 0.5)
-                p[0]["z"] = zmin + dz * (iz + 0.5)
-
-                # query and predict
-                lower, target, upper = query_weight_window(p, mcdc_obj, data)
-                exp_lower = 100 * ix + 10 * iy + iz + 1
-                exp_target = 1000 + exp_lower
-                exp_upper = 2000 + exp_lower
-
-                assert lower == exp_lower
-                assert target == exp_target
-                assert upper == exp_upper
+    for ne, energy in enumerate(energies):
+        for ix in range(nx):
+            for iy in range(ny):
+                for iz in range(nz):
+                    # put particle in center of current mesh bin
+                    p[0]["x"] = xmin + dx * (ix + 0.5)
+                    p[0]["y"] = ymin + dy * (iy + 0.5)
+                    p[0]["z"] = zmin + dz * (iz + 0.5)
+                    # assign energy to be in center of bins
+                    p[0]["E"] = energy
+
+                    # query and predict
+                    lower, target, upper = query_weight_window(p, mcdc_obj, data)
+                    exp_lower = 1000* ne + 100 * ix + 10 * iy + iz + 1
+                    exp_target = 10000 + exp_lower
+                    exp_upper = 20000 + exp_lower
+
+                    assert lower == exp_lower
+                    assert target == exp_target
+                    assert upper == exp_upper

From 6bc20f3883e25ec97ca50abb51186e3a4558eadc Mon Sep 17 00:00:00 2001
From: Nathan Glaser <glaserna@oregonstate.edu>
Date: Tue, 14 Apr 2026 12:06:37 -0700
Subject: [PATCH 23/53] add code factory gen for ww

---
 mcdc/mcdc_get/weight_windows.py | 29 +++++++++++++++++++++++++++++
 mcdc/mcdc_set/weight_windows.py | 29 +++++++++++++++++++++++++++++
 mcdc/numba_types.py             |  3 +++
 3 files changed, 61 insertions(+)

diff --git a/mcdc/mcdc_get/weight_windows.py b/mcdc/mcdc_get/weight_windows.py
index 00252e542..8e80a196a 100644
--- a/mcdc/mcdc_get/weight_windows.py
+++ b/mcdc/mcdc_get/weight_windows.py
@@ -3,6 +3,35 @@
 from numba import njit
 
 
+@njit
+def energy_bounds(index, weight_windows, data):
+    offset = weight_windows["energy_bounds_offset"]
+    return data[offset + index]
+
+
+@njit
+def energy_bounds_all(weight_windows, data):
+    start = weight_windows["energy_bounds_offset"]
+    size = weight_windows["energy_bounds_length"]
+    end = start + size
+    return data[start:end]
+
+
+@njit
+def energy_bounds_last(weight_windows, data):
+    start = weight_windows["energy_bounds_offset"]
+    size = weight_windows["energy_bounds_length"]
+    end = start + size
+    return data[end - 1]
+
+
+@njit
+def energy_bounds_chunk(start, length, weight_windows, data):
+    start += weight_windows["energy_bounds_offset"]
+    end = start + length
+    return data[start:end]
+
+
 @njit
 def lower_weights(index, weight_windows, data):
     offset = weight_windows["lower_weights_offset"]
diff --git a/mcdc/mcdc_set/weight_windows.py b/mcdc/mcdc_set/weight_windows.py
index 95613786f..ca35418c6 100644
--- a/mcdc/mcdc_set/weight_windows.py
+++ b/mcdc/mcdc_set/weight_windows.py
@@ -3,6 +3,35 @@
 from numba import njit
 
 
+@njit
+def energy_bounds(index, weight_windows, data, value):
+    offset = weight_windows["energy_bounds_offset"]
+    data[offset + index] = value
+
+
+@njit
+def energy_bounds_all(weight_windows, data, value):
+    start = weight_windows["energy_bounds_offset"]
+    size = weight_windows["energy_bounds_length"]
+    end = start + size
+    data[start:end] = value
+
+
+@njit
+def energy_bounds_last(weight_windows, data, value):
+    start = weight_windows["energy_bounds_offset"]
+    size = weight_windows["energy_bounds_length"]
+    end = start + size
+    data[end - 1] = value
+
+
+@njit
+def energy_bounds_chunk(start, length, weight_windows, data, value):
+    start += weight_windows["energy_bounds_offset"]
+    end = start + length
+    data[start:end] = value
+
+
 @njit
 def lower_weights(index, weight_windows, data, value):
     offset = weight_windows["lower_weights_offset"]
diff --git a/mcdc/numba_types.py b/mcdc/numba_types.py
index 46aa629d6..c0e2c1871 100644
--- a/mcdc/numba_types.py
+++ b/mcdc/numba_types.py
@@ -576,6 +576,9 @@
 
 weight_windows = into_dtype([
     ('active', bool),
+    ('energy_bounds_offset', int64),
+    ('energy_bounds_length', int64),
+    ('Ne', int64),
     ('mesh_ID', int64),
     ('Nx', int64),
     ('Ny', int64),

From b2ebd3c867d9e458f9675c00f7d4be476ecb5bcc Mon Sep 17 00:00:00 2001
From: Nathan Glaser <glaserna@oregonstate.edu>
Date: Tue, 14 Apr 2026 12:08:27 -0700
Subject: [PATCH 24/53] back in black

---
 mcdc/object_/technique.py                       | 4 +---
 mcdc/transport/technique.py                     | 6 ++++--
 test/unit/transport/technique/weight_windows.py | 2 +-
 3 files changed, 6 insertions(+), 6 deletions(-)

diff --git a/mcdc/object_/technique.py b/mcdc/object_/technique.py
index 381b031e9..a75b7f3f1 100644
--- a/mcdc/object_/technique.py
+++ b/mcdc/object_/technique.py
@@ -119,9 +119,7 @@ def __call__(self, mesh, weight_windows, energy=np.array([-0.5, INF])):
                 )
         # validate energy as strictly increasing
         if not (np.diff(energy) > 0).all():
-            print_error(
-                "Energy bounds must be strictly increasing"
-            )
+            print_error("Energy bounds must be strictly increasing")
         # get energy size
         if len(energy.shape) != 1:
             print_error(
diff --git a/mcdc/transport/technique.py b/mcdc/transport/technique.py
index 0a4edb80c..1e3e2a712 100644
--- a/mcdc/transport/technique.py
+++ b/mcdc/transport/technique.py
@@ -79,19 +79,20 @@ def query_weight_window(particle_container, simulation, data):
     """
     # grab objects
     ww_obj = simulation["weight_windows"]
-    index = get_ww_index(particle_container, ww_obj, simulation, data) 
+    index = get_ww_index(particle_container, ww_obj, simulation, data)
     # grab the actual ww parameters
     lower = ww_get.lower_weights(index, ww_obj, data)
     target = ww_get.target_weights(index, ww_obj, data)
     upper = ww_get.upper_weights(index, ww_obj, data)
     return lower, target, upper
 
+
 @njit
 def get_ww_index(particle_container, ww_obj, simulation, data):
     particle = particle_container[0]
 
     # get energy index
-    energy_bounds = ww_get.energy_bounds_all(ww_obj, data) 
+    energy_bounds = ww_get.energy_bounds_all(ww_obj, data)
     if simulation["settings"]["neutron_multigroup_mode"]:
         energy = particle["g"]
     else:
@@ -105,6 +106,7 @@ def get_ww_index(particle_container, ww_obj, simulation, data):
     index = (((ie * ww_obj["Nx"]) + idx) * ww_obj["Ny"] + idy) * ww_obj["Nz"] + idz
     return index
 
+
 @njit
 def split_from_weight_window(particle_container, threshold_weight, simulation):
     """
diff --git a/test/unit/transport/technique/weight_windows.py b/test/unit/transport/technique/weight_windows.py
index 9991c37c9..347fb92d5 100644
--- a/test/unit/transport/technique/weight_windows.py
+++ b/test/unit/transport/technique/weight_windows.py
@@ -194,7 +194,7 @@ def test_query_weight_window():
 
                     # query and predict
                     lower, target, upper = query_weight_window(p, mcdc_obj, data)
-                    exp_lower = 1000* ne + 100 * ix + 10 * iy + iz + 1
+                    exp_lower = 1000 * ne + 100 * ix + 10 * iy + iz + 1
                     exp_target = 10000 + exp_lower
                     exp_upper = 20000 + exp_lower
 

From dd8929e38c9bf96fd88290ac3e68fe332a63eb2e Mon Sep 17 00:00:00 2001
From: Nathan Glaser <glaserna@oregonstate.edu>
Date: Tue, 14 Apr 2026 12:12:17 -0700
Subject: [PATCH 25/53] add docstring on ww index

---
 mcdc/transport/technique.py | 19 +++++++++++++++++++
 1 file changed, 19 insertions(+)

diff --git a/mcdc/transport/technique.py b/mcdc/transport/technique.py
index 1e3e2a712..e87bcefa9 100644
--- a/mcdc/transport/technique.py
+++ b/mcdc/transport/technique.py
@@ -89,6 +89,25 @@ def query_weight_window(particle_container, simulation, data):
 
 @njit
 def get_ww_index(particle_container, ww_obj, simulation, data):
+    """
+    Get flattened weight window index from particle information
+
+    Parameters
+    ----------
+    particle_container : ndarray
+        Container holding the particle.
+    weight_window_object : object
+        The weight window object containing index information
+    simulation : object
+        Simulation state containing weight window and mesh data.
+    data : object
+        Simulation data for array access. 
+
+    Returns
+    -------
+    index: int
+        the flattened index in the weight window array
+    """
     particle = particle_container[0]
 
     # get energy index

From 93b2c276595b7df0396968986e3b0df21a53ca24 Mon Sep 17 00:00:00 2001
From: Nathan Glaser <glaserna@oregonstate.edu>
Date: Tue, 14 Apr 2026 12:37:45 -0700
Subject: [PATCH 26/53] add none defaults to mesh and energy for ww

---
 docs/source/pythonapi/index.rst                 | 2 +-
 mcdc/object_/technique.py                       | 9 ++++++++-
 test/unit/transport/technique/weight_windows.py | 4 ++--
 3 files changed, 11 insertions(+), 4 deletions(-)

diff --git a/docs/source/pythonapi/index.rst b/docs/source/pythonapi/index.rst
index c16a2fa5f..a0f9a6445 100644
--- a/docs/source/pythonapi/index.rst
+++ b/docs/source/pythonapi/index.rst
@@ -89,7 +89,7 @@ Techniques are enabled by calling methods on the ``mcdc.simulation`` singleton:
 - ``mcdc.simulation.implicit_capture(active=True)``
 - ``mcdc.simulation.weighted_emission(active=True, weight_target=1.0)``
 - ``mcdc.simulation.weight_roulette(weight_threshold=0.0, weight_target=1.0)``
-- ``mcdc.simulation.weight_windows(mesh, weight_windows)``
+- ``mcdc.simulation.weight_windows(weight_windows, mesh=None, energy=None)``
 - ``mcdc.simulation.population_control(active=True)``
 
 Running
diff --git a/mcdc/object_/technique.py b/mcdc/object_/technique.py
index a75b7f3f1..f1ebe2922 100644
--- a/mcdc/object_/technique.py
+++ b/mcdc/object_/technique.py
@@ -102,7 +102,14 @@ def __init__(self):
         self.target_weights = np.array([1.0])
         self.upper_weights = np.array([1.0])
 
-    def __call__(self, mesh, weight_windows, energy=np.array([-0.5, INF])):
+    def __call__(self, weight_windows, mesh=None, energy=None):
+        # fill in defaults
+        if mesh is None:
+            mesh = MeshUniform()
+        if energy is None:
+            # usable for both groups and max energy
+            energy = np.array([-0.5, INF])
+
         # get mesh size
         match mesh.label:
             case "uniform_mesh":
diff --git a/test/unit/transport/technique/weight_windows.py b/test/unit/transport/technique/weight_windows.py
index 347fb92d5..2bcf54ddb 100644
--- a/test/unit/transport/technique/weight_windows.py
+++ b/test/unit/transport/technique/weight_windows.py
@@ -51,7 +51,7 @@ def make_ww_model_params(lower=0.1, target=1.0, upper=1.0, mess_up_size=False):
         ww_array[..., 1] = target
         ww_array[..., 2] = upper
 
-    mcdc.simulation.weight_windows(mesh, ww_array)
+    mcdc.simulation.weight_windows(ww_array, mesh=mesh)
 
     mcdc_container, data = preparation()
     return mcdc_container[0], data
@@ -76,7 +76,7 @@ def make_ww_model_distinct():
                     ww_array[e, i, j, k, 1] = 10000 + val
                     ww_array[e, i, j, k, 2] = 20000 + val
 
-    mcdc.simulation.weight_windows(mesh, ww_array, energy)
+    mcdc.simulation.weight_windows(ww_array, mesh=mesh, energy=energy)
 
     mcdc_container, data = preparation()
     return mcdc_container[0], data

From dcb9c8d5d02844da21e113c36140180b6a35c363 Mon Sep 17 00:00:00 2001
From: Nathan Glaser <glaserna@oregonstate.edu>
Date: Tue, 14 Apr 2026 14:36:04 -0700
Subject: [PATCH 27/53] reformat with black

---
 mcdc/transport/technique.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/mcdc/transport/technique.py b/mcdc/transport/technique.py
index e87bcefa9..25fe77569 100644
--- a/mcdc/transport/technique.py
+++ b/mcdc/transport/technique.py
@@ -101,7 +101,7 @@ def get_ww_index(particle_container, ww_obj, simulation, data):
     simulation : object
         Simulation state containing weight window and mesh data.
     data : object
-        Simulation data for array access. 
+        Simulation data for array access.
 
     Returns
     -------

From a21ef700238b24d69c63fcd2a6d2f197b1ce6a79 Mon Sep 17 00:00:00 2001
From: Nathan Glaser <123763273+nglaser3@users.noreply.github.com>
Date: Wed, 15 Apr 2026 17:38:36 -0700
Subject: [PATCH 28/53] Update mcdc/transport/simulation.py

Co-authored-by: Ilham Variansyah <variansi@oregonstate.edu>
---
 mcdc/transport/simulation.py | 9 ++++++++-
 1 file changed, 8 insertions(+), 1 deletion(-)

diff --git a/mcdc/transport/simulation.py b/mcdc/transport/simulation.py
index 74d3849d2..7cf6bd247 100644
--- a/mcdc/transport/simulation.py
+++ b/mcdc/transport/simulation.py
@@ -347,7 +347,14 @@ def step_particle(particle_container, program, data):
         particle["alive"] = False
 
     # Weight splitting / rouletting
-    if particle["alive"] & simulation["weight_windows"]["active"]:
+    # Apply techniques
+    if particle["alive"]:
+        # Weight windows
+        if simulation["weight_windows"]["active"]:
+            technique.weight_windows(particle_container, simulation, data)
+        
+        # Weight roulette
+        else:
         technique.weight_windows(particle_container, simulation, data)
     elif particle["alive"]:
         technique.weight_roulette(particle_container, simulation)

From 1d039ebcc5d0da754e0b679584ccc89aafbfaa74 Mon Sep 17 00:00:00 2001
From: Nathan Glaser <glaserna@oregonstate.edu>
Date: Wed, 15 Apr 2026 17:47:06 -0700
Subject: [PATCH 29/53] Switch from simulation to program for split call

---
 mcdc/transport/simulation.py                    |  6 ++----
 mcdc/transport/technique.py                     | 15 ++++++++-------
 test/unit/transport/technique/weight_windows.py | 16 +++++++++-------
 3 files changed, 19 insertions(+), 18 deletions(-)

diff --git a/mcdc/transport/simulation.py b/mcdc/transport/simulation.py
index 7cf6bd247..8807977d1 100644
--- a/mcdc/transport/simulation.py
+++ b/mcdc/transport/simulation.py
@@ -351,13 +351,11 @@ def step_particle(particle_container, program, data):
     if particle["alive"]:
         # Weight windows
         if simulation["weight_windows"]["active"]:
-            technique.weight_windows(particle_container, simulation, data)
+            technique.weight_windows(particle_container, program, data)
         
         # Weight roulette
         else:
-        technique.weight_windows(particle_container, simulation, data)
-    elif particle["alive"]:
-        technique.weight_roulette(particle_container, simulation)
+            technique.weight_windows(particle_container, simulation, data)
 
 
 @njit
diff --git a/mcdc/transport/technique.py b/mcdc/transport/technique.py
index 25fe77569..ba4895d80 100644
--- a/mcdc/transport/technique.py
+++ b/mcdc/transport/technique.py
@@ -36,7 +36,7 @@ def weight_roulette(particle_container, simulation):
 
 
 @njit
-def weight_windows(particle_container, simulation, data):
+def weight_windows(particle_container, program, data):
     """
     Apply weight window splitting and rouletting to a particle.
 
@@ -44,11 +44,12 @@ def weight_windows(particle_container, simulation, data):
     ----------
     particle_container : ndarray
         Container holding the particle.
-    simulation : object
-        Simulation state containing weight window and mesh data.
+    program : object
+        Program object containing simulation state with weight window and mesh data.
     data : object
         Simulation data for array access.
     """
+    simulation = util.access_simulation(program)
     [lower, target, upper] = query_weight_window(particle_container, simulation, data)
     split_from_weight_window(particle_container, upper, simulation)
     roulette_from_weight_window(particle_container, lower, target)
@@ -127,7 +128,7 @@ def get_ww_index(particle_container, ww_obj, simulation, data):
 
 
 @njit
-def split_from_weight_window(particle_container, threshold_weight, simulation):
+def split_from_weight_window(particle_container, threshold_weight, program):
     """
     Split a particle if its weight exceeds the threshold.
 
@@ -137,8 +138,8 @@ def split_from_weight_window(particle_container, threshold_weight, simulation):
         Container holding the particle.
     threshold_weight : float
         Upper weight bound triggering splitting.
-    simulation : object
-        Simulation state used for banking split particles.
+    program : object
+        Program object used for banking split particles.
     """
     particle = particle_container[0]
     weight = particle["w"]
@@ -149,7 +150,7 @@ def split_from_weight_window(particle_container, threshold_weight, simulation):
         particle["w"] = weight / num_split
         for _ in range(num_split - 1):
             # bank split particles into the active bank
-            particle_bank_module.bank_active_particle(particle_container, simulation)
+            particle_bank_module.bank_active_particle(particle_container, program)
 
 
 @njit
diff --git a/test/unit/transport/technique/weight_windows.py b/test/unit/transport/technique/weight_windows.py
index 2bcf54ddb..3693cfe93 100644
--- a/test/unit/transport/technique/weight_windows.py
+++ b/test/unit/transport/technique/weight_windows.py
@@ -14,6 +14,7 @@
 )
 from mcdc.transport.mesh.interface import get_indices
 from mcdc.constant import TINY
+import mcdc.transport.util as util
 
 # =========================================================================== #
 # Helper method for creating a dummy model
@@ -143,14 +144,15 @@ def test_split_from_weight_window():
     init_weight = 2.0 + TINY
     particles[0]["w"] = init_weight
     threshold = 1.0
-    mcdc_obj, data = make_ww_model_distinct()
+    program, data = make_ww_model_distinct()
+    simulation = util.access_simulation(program)
 
     # get bank and init size
-    bank = mcdc_obj["bank_active"]
+    bank = simulation["bank_active"]
     init_bank_size = particle_bank_module.get_bank_size(bank)
 
     # split
-    split_from_weight_window(particles, threshold, mcdc_obj)
+    split_from_weight_window(particles, threshold, program)
 
     p1 = particles[0]
     num_split = np.ceil(init_weight / threshold)
@@ -168,9 +170,9 @@ def test_split_from_weight_window():
 def test_query_weight_window():
     p = np.zeros(1, type_.particle_data)
 
-    mcdc_obj, data = make_ww_model_distinct()
-    mcdc_obj["settings"]["neutron_multigroup_mode"] = False
-
+    program, data = make_ww_model_distinct()
+    simulation = util.access_simulation(program)
+    simulation["settings"]["neutron_multigroup_mode"] = False
     # hardcode mesh params
     pitch, height, N = 2.0, 10.0, 3
     nx, ny, nz = N, N, N
@@ -193,7 +195,7 @@ def test_query_weight_window():
                     p[0]["E"] = energy
 
                     # query and predict
-                    lower, target, upper = query_weight_window(p, mcdc_obj, data)
+                    lower, target, upper = query_weight_window(p, simulation, data)
                     exp_lower = 1000 * ne + 100 * ix + 10 * iy + iz + 1
                     exp_target = 10000 + exp_lower
                     exp_upper = 20000 + exp_lower

From 38004804a85f52180ac1f66c67c1c30b7ac6410d Mon Sep 17 00:00:00 2001
From: Nathan Glaser <glaserna@oregonstate.edu>
Date: Wed, 15 Apr 2026 19:52:00 -0700
Subject: [PATCH 30/53] add 4d accessor

---
 mcdc/code_factory/numba_objects_generator.py | 27 ++++++++++++++++++++
 1 file changed, 27 insertions(+)

diff --git a/mcdc/code_factory/numba_objects_generator.py b/mcdc/code_factory/numba_objects_generator.py
index d07da87ad..b03e05eb4 100644
--- a/mcdc/code_factory/numba_objects_generator.py
+++ b/mcdc/code_factory/numba_objects_generator.py
@@ -1013,6 +1013,14 @@ def generate_mcdc_access(targets):
                     object_name, attribute_name, shape[1], shape[2], True
                 )
 
+            elif len(shape) == 4:
+                text_getter += _accessor_4d_element(
+                    object_name, attribute_name, shape[1], shape[2], shape[3]
+                )
+
+                text_setter += _accessor_4d_element(
+                    object_name, attribute_name, shape[1], shape[2], shape[3], True
+                )
             text_getter += _accessor_chunk(object_name, attribute_name)
             text_setter += _accessor_chunk(object_name, attribute_name, True)
 
@@ -1157,6 +1165,25 @@ def _accessor_3d_element(object_name, attribute_name, stride_2, stride_3, setter
     return text
 
 
+def _accessor_4d_element(object_name, attribute_name, stride_2, stride_3, stride_4, setter=False):
+    text = f"@njit\n"
+    if setter:
+        text += f"def {attribute_name}(index_1, index_2, index_3, index_4, {object_name}, data, value):\n"
+    else:
+        text += (
+            f"def {attribute_name}(index_1, index_2, index_3, index_4, {object_name}, data):\n"
+        )
+    text += f'    offset = {object_name}["{attribute_name}_offset"]\n'
+    text += f'    stride_2 = {object_name}["{stride_2}"]\n'
+    text += f'    stride_3 = {object_name}["{stride_3}"]\n'
+    text += f'    stride_4 = {object_name}["{stride_4}"]\n'
+    if setter:
+        text += f"    data[offset + index_1 * stride_2 * stride_3 * stride_4 + index_2 * stride_3 * stride_4 + index_3 * stride_4 + index_4] = value\n\n\n"
+    else:
+        text += f"    return data[offset + index_1 * stride_2 * stride_3 * stride_4 + index_2 * stride_3 * stride_4 + index_3 * stride_4 + index_4]\n\n\n"
+    return text
+
+
 # ======================================================================================
 # Misc.
 # ======================================================================================

From b403037259fd31fd27b2ee7a77c622dccf6d085d Mon Sep 17 00:00:00 2001
From: Nathan Glaser <glaserna@oregonstate.edu>
Date: Wed, 15 Apr 2026 19:52:27 -0700
Subject: [PATCH 31/53] use annotate for indexing

---
 mcdc/mcdc_get/weight_windows.py | 69 +++++++--------------------------
 mcdc/mcdc_set/weight_windows.py | 69 +++++++--------------------------
 mcdc/object_/technique.py       | 24 +++++++-----
 3 files changed, 44 insertions(+), 118 deletions(-)

diff --git a/mcdc/mcdc_get/weight_windows.py b/mcdc/mcdc_get/weight_windows.py
index 8e80a196a..b1b7cbece 100644
--- a/mcdc/mcdc_get/weight_windows.py
+++ b/mcdc/mcdc_get/weight_windows.py
@@ -33,25 +33,12 @@ def energy_bounds_chunk(start, length, weight_windows, data):
 
 
 @njit
-def lower_weights(index, weight_windows, data):
+def lower_weights(index_1, index_2, index_3, index_4, weight_windows, data):
     offset = weight_windows["lower_weights_offset"]
-    return data[offset + index]
-
-
-@njit
-def lower_weights_all(weight_windows, data):
-    start = weight_windows["lower_weights_offset"]
-    size = weight_windows["lower_weights_length"]
-    end = start + size
-    return data[start:end]
-
-
-@njit
-def lower_weights_last(weight_windows, data):
-    start = weight_windows["lower_weights_offset"]
-    size = weight_windows["lower_weights_length"]
-    end = start + size
-    return data[end - 1]
+    stride_2 = weight_windows["Nx"]
+    stride_3 = weight_windows["Ny"]
+    stride_4 = weight_windows["Nz"]
+    return data[offset + index_1 * stride_2 * stride_3 * stride_4 + index_2 * stride_3 * stride_4 + index_3 * stride_4 + index_4]
 
 
 @njit
@@ -62,25 +49,12 @@ def lower_weights_chunk(start, length, weight_windows, data):
 
 
 @njit
-def target_weights(index, weight_windows, data):
+def target_weights(index_1, index_2, index_3, index_4, weight_windows, data):
     offset = weight_windows["target_weights_offset"]
-    return data[offset + index]
-
-
-@njit
-def target_weights_all(weight_windows, data):
-    start = weight_windows["target_weights_offset"]
-    size = weight_windows["target_weights_length"]
-    end = start + size
-    return data[start:end]
-
-
-@njit
-def target_weights_last(weight_windows, data):
-    start = weight_windows["target_weights_offset"]
-    size = weight_windows["target_weights_length"]
-    end = start + size
-    return data[end - 1]
+    stride_2 = weight_windows["Nx"]
+    stride_3 = weight_windows["Ny"]
+    stride_4 = weight_windows["Nz"]
+    return data[offset + index_1 * stride_2 * stride_3 * stride_4 + index_2 * stride_3 * stride_4 + index_3 * stride_4 + index_4]
 
 
 @njit
@@ -91,25 +65,12 @@ def target_weights_chunk(start, length, weight_windows, data):
 
 
 @njit
-def upper_weights(index, weight_windows, data):
+def upper_weights(index_1, index_2, index_3, index_4, weight_windows, data):
     offset = weight_windows["upper_weights_offset"]
-    return data[offset + index]
-
-
-@njit
-def upper_weights_all(weight_windows, data):
-    start = weight_windows["upper_weights_offset"]
-    size = weight_windows["upper_weights_length"]
-    end = start + size
-    return data[start:end]
-
-
-@njit
-def upper_weights_last(weight_windows, data):
-    start = weight_windows["upper_weights_offset"]
-    size = weight_windows["upper_weights_length"]
-    end = start + size
-    return data[end - 1]
+    stride_2 = weight_windows["Nx"]
+    stride_3 = weight_windows["Ny"]
+    stride_4 = weight_windows["Nz"]
+    return data[offset + index_1 * stride_2 * stride_3 * stride_4 + index_2 * stride_3 * stride_4 + index_3 * stride_4 + index_4]
 
 
 @njit
diff --git a/mcdc/mcdc_set/weight_windows.py b/mcdc/mcdc_set/weight_windows.py
index ca35418c6..e6a2ab8a2 100644
--- a/mcdc/mcdc_set/weight_windows.py
+++ b/mcdc/mcdc_set/weight_windows.py
@@ -33,25 +33,12 @@ def energy_bounds_chunk(start, length, weight_windows, data, value):
 
 
 @njit
-def lower_weights(index, weight_windows, data, value):
+def lower_weights(index_1, index_2, index_3, index_4, weight_windows, data, value):
     offset = weight_windows["lower_weights_offset"]
-    data[offset + index] = value
-
-
-@njit
-def lower_weights_all(weight_windows, data, value):
-    start = weight_windows["lower_weights_offset"]
-    size = weight_windows["lower_weights_length"]
-    end = start + size
-    data[start:end] = value
-
-
-@njit
-def lower_weights_last(weight_windows, data, value):
-    start = weight_windows["lower_weights_offset"]
-    size = weight_windows["lower_weights_length"]
-    end = start + size
-    data[end - 1] = value
+    stride_2 = weight_windows["Nx"]
+    stride_3 = weight_windows["Ny"]
+    stride_4 = weight_windows["Nz"]
+    data[offset + index_1 * stride_2 * stride_3 * stride_4 + index_2 * stride_3 * stride_4 + index_3 * stride_4 + index_4] = value
 
 
 @njit
@@ -62,25 +49,12 @@ def lower_weights_chunk(start, length, weight_windows, data, value):
 
 
 @njit
-def target_weights(index, weight_windows, data, value):
+def target_weights(index_1, index_2, index_3, index_4, weight_windows, data, value):
     offset = weight_windows["target_weights_offset"]
-    data[offset + index] = value
-
-
-@njit
-def target_weights_all(weight_windows, data, value):
-    start = weight_windows["target_weights_offset"]
-    size = weight_windows["target_weights_length"]
-    end = start + size
-    data[start:end] = value
-
-
-@njit
-def target_weights_last(weight_windows, data, value):
-    start = weight_windows["target_weights_offset"]
-    size = weight_windows["target_weights_length"]
-    end = start + size
-    data[end - 1] = value
+    stride_2 = weight_windows["Nx"]
+    stride_3 = weight_windows["Ny"]
+    stride_4 = weight_windows["Nz"]
+    data[offset + index_1 * stride_2 * stride_3 * stride_4 + index_2 * stride_3 * stride_4 + index_3 * stride_4 + index_4] = value
 
 
 @njit
@@ -91,25 +65,12 @@ def target_weights_chunk(start, length, weight_windows, data, value):
 
 
 @njit
-def upper_weights(index, weight_windows, data, value):
+def upper_weights(index_1, index_2, index_3, index_4, weight_windows, data, value):
     offset = weight_windows["upper_weights_offset"]
-    data[offset + index] = value
-
-
-@njit
-def upper_weights_all(weight_windows, data, value):
-    start = weight_windows["upper_weights_offset"]
-    size = weight_windows["upper_weights_length"]
-    end = start + size
-    data[start:end] = value
-
-
-@njit
-def upper_weights_last(weight_windows, data, value):
-    start = weight_windows["upper_weights_offset"]
-    size = weight_windows["upper_weights_length"]
-    end = start + size
-    data[end - 1] = value
+    stride_2 = weight_windows["Nx"]
+    stride_3 = weight_windows["Ny"]
+    stride_4 = weight_windows["Nz"]
+    data[offset + index_1 * stride_2 * stride_3 * stride_4 + index_2 * stride_3 * stride_4 + index_3 * stride_4 + index_4] = value
 
 
 @njit
diff --git a/mcdc/object_/technique.py b/mcdc/object_/technique.py
index f1ebe2922..12435226c 100644
--- a/mcdc/object_/technique.py
+++ b/mcdc/object_/technique.py
@@ -1,5 +1,6 @@
 from numpy.typing import NDArray
 import numpy as np
+from typing import Annotated
 from mcdc.constant import INF
 from mcdc.object_.base import ObjectSingleton
 from mcdc.object_.mesh import MeshBase, MeshUniform
@@ -87,10 +88,10 @@ class WeightWindows(ObjectSingleton):
     Ny: int
     Nz: int
 
-    # flattened arrays of ww params
-    lower_weights: NDArray[np.float64]
-    target_weights: NDArray[np.float64]
-    upper_weights: NDArray[np.float64]
+    # arrays of ww params
+    lower_weights: Annotated[NDArray[np.float64], ("Ne", "Nx", "Ny", "Nz")]
+    target_weights: Annotated[NDArray[np.float64], ("Ne", "Nx", "Ny", "Nz")]
+    upper_weights: Annotated[NDArray[np.float64], ("Ne", "Nx", "Ny", "Nz")]
 
     def __init__(self):
         self.active = False
@@ -98,9 +99,11 @@ def __init__(self):
         self.Ne = 1
         self.mesh_ID = -1  # skirt around having to create a MeshBase instance
         self.Nx, self.Ny, self.Nz = 1, 1, 1
-        self.lower_weights = np.array([1.0])
-        self.target_weights = np.array([1.0])
-        self.upper_weights = np.array([1.0])
+        self.Nt = 1
+        shape = (self.Ne, self.Nx, self.Ny, self.Nz)
+        self.lower_weights = np.array([1.0]).reshape(*shape)
+        self.target_weights = np.array([1.0]).reshape(*shape)
+        self.upper_weights = np.array([1.0]).reshape(*shape)
 
     def __call__(self, weight_windows, mesh=None, energy=None):
         # fill in defaults
@@ -148,9 +151,10 @@ def __call__(self, weight_windows, mesh=None, energy=None):
         self.Ne = ne
         self.mesh = mesh
         self.Nx, self.Ny, self.Nz = mesh_shape
-        self.lower_weights = weight_windows[..., 0].reshape(-1)
-        self.target_weights = weight_windows[..., 1].reshape(-1)
-        self.upper_weights = weight_windows[..., 2].reshape(-1)
+        shape = (self.Ne, *mesh_shape)
+        self.lower_weights = weight_windows[..., 0]
+        self.target_weights = weight_windows[..., 1]
+        self.upper_weights = weight_windows[..., 2]
 
         # check weight windows are valid
         if (self.lower_weights <= 0.0).any():

From 481737af8054e5695fc1f2899f578a7aeb6718e2 Mon Sep 17 00:00:00 2001
From: Nathan Glaser <glaserna@oregonstate.edu>
Date: Wed, 15 Apr 2026 19:53:53 -0700
Subject: [PATCH 32/53] change up indices getter

---
 mcdc/transport/technique.py | 15 +++++++--------
 1 file changed, 7 insertions(+), 8 deletions(-)

diff --git a/mcdc/transport/technique.py b/mcdc/transport/technique.py
index ba4895d80..51edf85e8 100644
--- a/mcdc/transport/technique.py
+++ b/mcdc/transport/technique.py
@@ -80,16 +80,16 @@ def query_weight_window(particle_container, simulation, data):
     """
     # grab objects
     ww_obj = simulation["weight_windows"]
-    index = get_ww_index(particle_container, ww_obj, simulation, data)
+    indices = get_ww_indices(particle_container, ww_obj, simulation, data)
     # grab the actual ww parameters
-    lower = ww_get.lower_weights(index, ww_obj, data)
-    target = ww_get.target_weights(index, ww_obj, data)
-    upper = ww_get.upper_weights(index, ww_obj, data)
+    lower = ww_get.lower_weights(*indices, ww_obj, data)
+    target = ww_get.target_weights(*indices, ww_obj, data)
+    upper = ww_get.upper_weights(*indices, ww_obj, data)
     return lower, target, upper
 
 
 @njit
-def get_ww_index(particle_container, ww_obj, simulation, data):
+def get_ww_indices(particle_container, ww_obj, simulation, data):
     """
     Get flattened weight window index from particle information
 
@@ -106,7 +106,7 @@ def get_ww_index(particle_container, ww_obj, simulation, data):
 
     Returns
     -------
-    index: int
+    indices: Tuple[int]
         the flattened index in the weight window array
     """
     particle = particle_container[0]
@@ -123,8 +123,7 @@ def get_ww_index(particle_container, ww_obj, simulation, data):
     mesh = simulation["meshes"][ww_obj["mesh_ID"]]
     idx, idy, idz = get_mesh_indices(particle_container, mesh, simulation, data)
 
-    index = (((ie * ww_obj["Nx"]) + idx) * ww_obj["Ny"] + idy) * ww_obj["Nz"] + idz
-    return index
+    return (ie, idx, idy, idz)
 
 
 @njit

From 2ca86ca91252b122987054ad9826d7586b939f6e Mon Sep 17 00:00:00 2001
From: Nathan Glaser <glaserna@oregonstate.edu>
Date: Wed, 15 Apr 2026 19:58:06 -0700
Subject: [PATCH 33/53] generalize weight roullete to use helper

---
 mcdc/transport/technique.py                     | 15 +++++----------
 test/unit/transport/technique/weight_windows.py |  6 +++---
 2 files changed, 8 insertions(+), 13 deletions(-)

diff --git a/mcdc/transport/technique.py b/mcdc/transport/technique.py
index 51edf85e8..8fdd70685 100644
--- a/mcdc/transport/technique.py
+++ b/mcdc/transport/technique.py
@@ -20,14 +20,9 @@
 
 @njit
 def weight_roulette(particle_container, simulation):
-    particle = particle_container[0]
-    if particle["w"] < simulation["weight_roulette"]["weight_threshold"]:
-        w_target = simulation["weight_roulette"]["weight_target"]
-        survival_probability = particle["w"] / w_target
-        if rng.lcg(particle_container) < survival_probability:
-            particle["w"] = w_target
-        else:
-            particle["alive"] = False
+    threshold = simulation["weight_roulette"]["weight_threshold"]
+    target = simulation["weight_roulette"]["weight_target"]
+    roulette_from_weight_bounds(particle_container, threshold, target)
 
 
 # ======================================================================================
@@ -52,7 +47,7 @@ def weight_windows(particle_container, program, data):
     simulation = util.access_simulation(program)
     [lower, target, upper] = query_weight_window(particle_container, simulation, data)
     split_from_weight_window(particle_container, upper, simulation)
-    roulette_from_weight_window(particle_container, lower, target)
+    roulette_from_weight_bounds(particle_container, lower, target)
 
 
 @njit
@@ -153,7 +148,7 @@ def split_from_weight_window(particle_container, threshold_weight, program):
 
 
 @njit
-def roulette_from_weight_window(particle_container, w_threshold, w_target):
+def roulette_from_weight_bounds(particle_container, w_threshold, w_target):
     """
     Russian roulette particle if weight is below threshold.
 
diff --git a/test/unit/transport/technique/weight_windows.py b/test/unit/transport/technique/weight_windows.py
index 3693cfe93..e29709cd2 100644
--- a/test/unit/transport/technique/weight_windows.py
+++ b/test/unit/transport/technique/weight_windows.py
@@ -6,7 +6,7 @@
 from mcdc.main import preparation
 import mcdc.numba_types as type_
 from mcdc.transport.technique import (
-    roulette_from_weight_window,
+    roulette_from_weight_bounds,
     split_from_weight_window,
     query_weight_window,
     weight_windows,
@@ -126,7 +126,7 @@ def test_error_throw(capsys, kwargs, expected_msg):
 # =========================================================================== #
 
 
-def test_roullete_from_weight_window():
+def test_roulette_from_weight_bounds():
     # because of rng, want to loop over to hit both branches
     for _ in range(10):
         particles = np.zeros(1, type_.particle_data)
@@ -134,7 +134,7 @@ def test_roullete_from_weight_window():
         target = 0.2
         threshold = 0.1 + TINY
 
-        roulette_from_weight_window(particles, threshold, target)
+        roulette_from_weight_bounds(particles, threshold, target)
         p = particles[0]
         assert p["w"] == target or not p["alive"]
 

From bbc8ef012acf0785150665b3764373b816239267 Mon Sep 17 00:00:00 2001
From: Nathan Glaser <glaserna@oregonstate.edu>
Date: Wed, 15 Apr 2026 20:25:53 -0700
Subject: [PATCH 34/53] add pincell regression test with weight windows

---
 .../regression/basic_weight_windows/answer.h5 | Bin 0 -> 39920 bytes
 test/regression/basic_weight_windows/input.py |  63 ++++++++++++++++++
 2 files changed, 63 insertions(+)
 create mode 100644 test/regression/basic_weight_windows/answer.h5
 create mode 100644 test/regression/basic_weight_windows/input.py

diff --git a/test/regression/basic_weight_windows/answer.h5 b/test/regression/basic_weight_windows/answer.h5
new file mode 100644
index 0000000000000000000000000000000000000000..93163bbc349c8468b636d0816d868c413af83ed4
GIT binary patch
literal 39920
zcmeHv30zLw_x{^}Oo?zaCnSYRq)bl=g-TbWvB}#&^CSwHk|~L(XjExXnoFWGWTp)D
zjtyp+LT2(mZ|5BD?UL*M@9%zdZ|D7d-nIAI>+G}FbJkvapLO;=n{}qm=-R1oCx+k<
z6Jv-nB$|@+-zVy*-$!!jJY8-<C4?v%K+!~DG9kne{1^-=ir<~8zqy>z-09O9gd|-*
zO{vouN}^;+jApda{}ci`)2A8;6gH!@R-))x%Fd`P!;9(WZtv{GaA$gW*gM&}(|iw0
z2M2qmJHyS>$-~}}*<_P~GC5k{8&rtVj+CVgT!=whAZ_m=gk=lTeT0OXEL+ff$TCu{
z9U(8s3k#EMF{&qs&=oc7ZwyB0H(Js<*oOLZrU?;Vg2t_EGubX0k2h3aY8%O?y-7g(
z-^)o+4k4Z5OOwS6yHUTFD{H#z3`YH2l24+_U`UHlC;X)3`*tPKWH1th!lcd`O`3kO
zC*dZ@zu1!mg(!^>4NNqhXw~#v<3w{iq5VKH)$ao+sx*mI?~k~Wum}<TQ*h9BN>jT3
zc9+IyoX2)%jQvE}*SGODlWwYzI*UvT{sX=we}a^t`w47Ja0oEb)|#W+gbmekX`HBT
z%1|1uJZ7{K;m?}Rply@!+27LLvK>|KU}?p4Akx(E&FDz+ot&M_oje__EU7##PhZE`
z)5FEngTbKll-Ort(d~1xEzyQR8v<<zv?0)jKpO&W2(%&4hCmwvZ3wg>(1t)80)G(!
zG9EAX&3L?8U$SB%6wq|!Q*k>NDz7w#Oq5ai{PAQ$nhGe?Pa*RbRNi|UnNOz!4ztL-
z7#-)CL*~8dzyq0W`~M#yK*xdUXrnhB!1+ddsNf{SuirNWD7G|3gZj5DN84TWfR^Q2
zd+8tJl_a8cHxMD-4+zqKJ`OA;P3lbF85)6TI?<}>-;am1G^ubK9s(Rhgyh=rAaMBd
z-bCpq$K&ho*L6$@C5gjs9!bCTex=5mEN9qIlp^#nOZ%bUmogYVDaL&|;NM)JDPO<l
ztL$$x>mlf-^a}~H5ki7@1~g0X9zihF(0Re-1@8#xeA6)W)kf12U*0Ix<IDS?`XZtd
z;;4tzsA)#;2X`!{s74c;(Nr#N9xbd%=9~VT^5VU2ebAG;T{Q14Q)$M~n?>>b`_@yW
zYl>>tSB`d6;}tE}<%cX%Vil2c*Nfts%FBwpn)qj_d~s296F;452ZNDF=0Vb8c$4x-
zAkPO?AH`>Zx*k<tgZX58D7o7-ZT0OWkQLy1f}+30`!~?HCEE~aL!b?THU$1}5cvD&
z;cC)kOWWq*0vtqyH2&pykn7@pd|i#sH*g6Rz}9MyzrVj2NR!I`Tl<S#Q(QxenBxnk
z>w?)iIc*nX%$TZcIGx}b%{3N~rh;#;!#W_{BtxA|rUm~2Uy`jL^{wZx=Ia~3C1+<y
zDiL5zQab*oM(LJ5uw^B6sA4I(mgUT-V!075%h7ROHTjn1T6^jD`|aOlrxwyA7NnJb
z!KJm?>7QHg{NsMzlG0!7i|S3O|C(R7{{{}x&2VX*S=#ReQO(m)Xt|w<RB_ptu*+9Z
z$fvoAbmHrjexD>kzX15USgYrL-|zfr+}~B2RJwWpif;8-2pahJ<FP`TRJ;uj0S+QU
zYo!JB%ZhC>$6!R$iB?U&Zy$g4d1JFQsd#fcY2D)z?T72ByC$vto!duEJG@m|mXo05
z0>-v1Cr!zPj%!)2wU_=`yZZ6{T|1>o-6fhq(7NsF`*{3a^Whi@pPzj`{CD9I`%~aS
z8?(N-rL~8>7xPQr%F@Hyj-hXE&2)12{Jigy$gi+=U=n<Lr{A8TZ|-90=3y@=YC6a=
zb1O@yW#(3%Ha1K*l26w5TfVuQrH8#UeO;<tu)TAIxs4ms)syLDy@J7@`GUQkf`jS6
zwDxc#HXV7mSvtA9IJ*&hiv)*_r-!E-^V{{?x_Fvfy037uHuqTJLR2NlyL&jhS=y5I
zl6j&@Ok3irMEQ!Q<~b3qTw?F(PMyVIEF<){C$35FhjL<iHZ_N7Z_9M@vUKodnmal#
zVKQjGqo;$1y{((Gr_1m8&6>|>zK@FBH)ZZ*N$Kt4=IrrDtr-kACtGuOCUc3P&4R<-
z(u?`qTjK8Q>1N#mE>yXVy#w>h#fW+Yn`_OjoIRbESh@)S@pNZW7C~U_&Rk;Nbf-vu
zGpq#7Xp$$|Z0p8!cPGxEtbpYI=9-_QJ!(y7{rNw7+*G67PNHhd=W5XINxfw`9ZF7O
zLd$X%l$`X$mgQP|=^w`lsir2=`z(pDUeiHSdY|QwkE3zYWK9$yGsJ(u=TvJ*{1Fdw
zK1rL;U)TSXG^t<PJVJn_uyFHz&(|q@XrxKBem?v7(dYdel-|vsTcT+IS~aEpus3zv
z%tp4{+#pI$O`~PGXi83JYRhu;I%6}<mgQP|>G%E4-!(twOOseo0DlJN=k&ur9>1=l
z^l#qZe)i|Nzxw@>yA&S(wsGmN#*tqAl=eaU;j>h?w$u9DvA>#NQF6u8zgqm6d7qLi
zoAH_T)dao2qelCy#Xk-6pYy|w(j+K9yC43$#-&efcn}8hWj#A%aGQPn77wAH0uRy;
zw_`|C56tyC<X!yMzY-;<Gy7W-BwvS;Gt+HZ&VrJ&(`#9d-hbz<-?H5Qs`Pu`5oDN<
zEo(lD`(L%U?YjRg0`&L8%&5ysQ&fUlUq!cuh_hNRb6`5zdf1U6OLL{BB6`&Q(YWLH
z?6PJadF3{{r0>9gi#E4QGMY;40;1OY(P-GZ**C@PZ5%wk1w7gwHE4k!q~y5s`D2uf
zpl07+Y`*UB2dn?jh<~qtt1Tp<_))*{p8@b+u4wcAtwHG<K>2-IeMe6+Wx3M6DZzAN
zx(VJaHhr{~;52o5f)cEzV<kmtKS28pip?PVxA5n$%-$bFmeBIe_kYskPjNH3=KD$M
z`u{8RkI(lc(%*aE_2Y3y;lU)P&EskTtqN+=t~z|SE5Ai8H;=Z<pv5iA(e@{|pk+B}
zsyN!bWjQ&jINhRUIdXQ#U=&-nET={lmszze_h(E0IlJodQ`pr%hlflX9)yv8i5rdi
zSK!f$>UcDz-@E^~AJ(MsNNm=XTbCdHUVou+py^-7qwZgU2NmY}{Qc^`jz_QzX&9}0
zoc`W(Nr3IwwP<GaZP|Y*0?nUG=y{s7&F>X<WLi22IPOcBUIHp4_#+ODrvD%R`=j3?
zXKF&^@}AG%m+47W-JJgQep_ZoR_WK}_pSE+efA&D(a&jO)OEd^;orPDwAV54Y+mx;
z&KHzgEmND%#Xw2D`~yiWsLE7)@M3Y0RI<yj9Mw=dDW%u&kc&N>k&V&S9P}7}dz?f!
zJ~}xi#@N?ZA)%z3xU&)$bI%Mb9_`7;gqMDz_uo~b-^umYW0$jGWZ_uc&54EHR|m;@
zd9fhBJWNMLo{iCcGW)8!vvFKB*42+wg#lfI183RuP*tFGcx81Brtyj{WW46%;E`X|
zw2$zReOB`Fj87bNPL&q>V8}<|71`LCNBF2qZS>83P=#}vnaaIx@$h)Owb<Mo9waxi
z8oi>bFmF`Ckho=BaLTWp_ps)}S0SfQw+|KIEggB)duJuaUdf-a;}Q$6nVH??&hsHL
z*xfY2hK((?I|3`?iTdmvSvLasI5qCY2d~pyyiA^`dDVr314FehY5G;cY4Z~0&G)$I
zvS-}TEggB->?7H?aA*b2*S(y*XDS!h<McVJs(BEzDC^5w#K&2;M^m=nVIfto*x1mJ
zgBwfhN-po`Kq_ngs9v49m>D)>nILbbnLgKm!-4dOUW=zpti-(%ooX3V9{hX1N>O;h
zhQH9`dGRhhq+fDTWOBJ!X*|J=UBX4+%=7OKzhdKlkk3KiA|4jTIxT8{hmVzAJ$<Bz
z>nod{41cM@g$7Thx`xp6-C)Ozu*F;)dY&6n@{kX;&HL(g`g2kBynAk}83&h*PuQi%
zu@NvcT}SS51@^?uiC?J8LB7{)=Z6dVSR6A<J=lwjk9}1xE{!AZ>+O+SCI@rTYi+^8
z34w(Ey>6!^^<<&suu?-JhYjJCFiQ$z<EE=bhj4QqmWvKFnqa~~(Swl*Z<q6+$cjDd
zcbbon0VhQ7Y~(@o^0d?`E2^;NOxn67{i+b%kQ`>0$wkn@0U;@&Y&5)R7<a0OjdlE+
zBDL0RTwc^~j;J#S!sA8Wtli5+&-dwx`&FvpBCEQ5Wfv}<D0R{|yUGJ6`&oGRc^sV7
zIW<Yun~NE?3G*dab8w<z7C&V<7Y^onV#iGQ$p1w({u<G)kz#io=h$%(uNY@Gi&uq!
z1>#bt_f_Fe&ug1w_VA(faLR}+W?cAJ44itXfQ^B&rYYluxky(YE_LB78`n1KMywsk
zM)1Xf7<8BemFb&J&PcGpxRPNhCd<cIE!W+nA8`@s*I&()xbI7@YKza6<3pol*{d_j
zJlxls9{2tv8=<#qh8Y=f;3qS0^rf}cu<fO}t*$#CXXd$WdD4>$k*n;toAmg&f7nU?
z%u*g!dv>XKpvA@e#`unI^*rQ9B$gigg@b2|TkmFz@t{yzpB`&P^k2^T^_!k^p<=xx
z*JuI<>%CR>oGRyF<nmJmSBdlIazkZuUsNJMStw3sG!LD9_0npyxG)Q^l<20xL;1Ne
zOW%|ccC_Z4b0AT_e}UQKz1z8HbUth_i0~&{$9g{27Ue<LzEhY9gAb>M6_2h><REaI
z;_EkwH*mFI!kTf?T+H9F@k8>gDp=jSI%TdiA4?CVc^pgU!EcwdeD}*d$hiq6I1_ep
zcEMmo?&IOP;&`=?<6KNNoE%-~SdIEQ{cp;>t;VLp)$<*<R>41qY28(pg9>kLm(oLg
z$fPOwq$?7*$VSv%vFBr`sHQ>+;lHlj6}mUDV-+5cZ675hS%Vp`SNHGymW{@#1wH1y
zChT~}1H-!^9Nb-ZXtQY!AG_~Ie2hxu!`j}$VPh>Hmk(wSUekdO;!9@R_lwuSap$tY
z2j{DB*Xx<VjqqypO*Lx2^B^CCj~zdG`ZkwnhbvcW5D&}b7xxKP<>RW>db$47*pTaY
zt-4y4i*lor3OQSN$Xlei>x@q|v?b#|Ju>D&JC;@7YZDhAuD?(?8_z*pyJstd*Hz<Y
z=P`NZZ}|vH-g<;Z;JC{4jZ(OJHJ0}=nVmGg21|eWblbXs4|!$1&eI6LHmeseS~I8y
zRlAaA_RZpBZqJWePX|?_^7-6%+urf<Zr;6`K6QjYHyYmQCE;HN7wq;jQs==`PkDPe
zlaJDvk4~8$e8i!QF~+I}7Ipd-K5|t!8-GVduQLzM!-vQF5q=`rc1O>T-fXP&-;gow
zG~vG<7{!MluYualh7VE(JRD6{oGBb#jaTkf2E~m;eM4(s_bTMUp`%gI9Ctp}rio70
zWLKlcwAN?HfEtKLaJEh`;o*sQyV7$O1da~s`=su1;4bQS?{psm7n_K!w|B7+d!bX;
zeWzJC^JM<ZtGl@HN}a9d5>pLMFUNQ>aUNn1zMC?%oEUGCRva4Lg^wlsb`@OI<e^Ny
zvM5=bi#r`eA5SLi%<|&o2e~7;i0Zhnzfv$C^$RR>wF0Ujz0)z=)v^ld&#xqFb*jP2
zq#<{oeB$Hg+>a0Ew5vkjUuH|QSM%U8>r8p+Y!=SSBuic({9~xH=&Ws$94yW5foU>a
z7!SL0Z$60oy(>WESNUq>I=xvrLz9ED+>IB+9Qn|X+2<}p_@^`%zSFu`E;Nn&doJ!k
zv{Sk0+A>uZ`oA#=EF8wdtnzaUi))B+WcnU?6XO0!=QYHSXkUX@W((p{m#`sZxqLuL
z3J-1>=Bs(iT$lvPNAEeyLE6^Rr)9MqoLw@!NS+usoOw@IF3@7b^3Islg@nCKn|bl5
zpBpj0T|Mz6I-(kpGs;u8SaT3FI-{#Pq5oV?m}C|)E;4%M`DA2qAX6*yx@rv<^3Lw%
zQ;2>yp;qeITysKi{XN&qRucXFSM#CM3H$H;;Az<gPd3iFo0ybZ^04y3sQL|$Ik0%t
zAav+a74+&1y~9EYKc|0mrgsSok#c$iiq5lf>}7nQxB?4h=lf|DhH+tWr?FI_Co!%L
zzqweiBOgZX!WaW;_*nODq0)&iJPe4xJo~Xr74$pk?T+rlgHmpzU0pXGt`=M{S?5KJ
zJEp#sRtX%uGwN^O_?(SdwRN{GE^tv7c5VHO3>JpGAGGppQ6;R5udq9+abP?6qOw*1
z3uekE*Lv2nP!wgBHEksulK0mbD9W(VMJrL+_cDj@)8%&)xy1Zuq*l8z0bCTn&e{7(
zzXs7<Rmbjos<C*RZoWe>4}o=43ylT#7(f5gsa(Q--)NpWMvT+4;qFPN*Ax9+)kFVK
z3Nark>0BN&jE4?`?j25=$wA1+unqO+xOlbC;no~?7F4#l<(gL$<5p3C@!Q>ORK9%P
zg;#b1Ro+n(x4x=GZbW=<-NS5XNX?2nx`>56Q&bxF>vJG-a!+{vejaWFd9&NQRwGeJ
zdhpw4d;|(LsK?%5qsPL$ZT9zgh|qOb6C(Wim=VV5TdY~=HR1Yuw`&}z?q0SR4IIdh
z$q1d-$c9wt-r0TyEUeCcs~%axMy%5g{thNF->6Z(7_pX(G)Zs6$eKz7C+h83piZ=B
zaABcE7dGT~7@aXmB<#?2c3P@G56QL3TSNx2G3ZwLJ(cG*kkQhAW=r_*bldiaiimmO
z9BZ!zQ4Kcsi}vZ)`6?ILM}L{0-I;|5*<W0)++kzFI*oqu0c>JEt1)LBkB!WmcJcBB
z94KFXI>bqb@Z%-nk*l|{kbh~9&0=CaF4@0i$#ud{*tqIMY@Sc_D@B!ulf~IstuVl*
zb`u+35vdZIp<J}jjm!xNt3t3-s*U+l7Lwb2D1N9_g>9piCbc8REsX&88>6Cl;1(8k
zJsr-6`OaM1WyO53rXSHhK-jyV>6FH_+)C8eaUajS%EnQ%6_Z{Oe(=8YU?|EH^EuTq
zTNe^`*hkwXD#@S{al89-)@|kB)i~|Q9A_^2_e1}f{yf;N*4i_TP0a7QGB4%K=3txs
zew{AF{7AR<>HGrWDqP4Nm^DwI825SUUiF0EVlUm;Q=1rX{C>3=@;;RfS>xBy>wI}A
z^jy$a)P@VyOBFGvR#srgF9ugiFR~G)S>PR{TZxL!*K6juu^{*27n`glTs-1+DD8EG
z1O0AWoIQ#A#?ZZI;c=9Yhf*PhV<+-4TH5Ah@~3KOz9^idO3Vx7N6M{P*vQ4=YvU7-
z94Go+>W2&O?o?x2cO`jo;(q?`o@)(OkxzcD>$&#7{WRL9eH#L82(%&4hQNOXf#8_9
zEDH$>sL0l@R#e!Hd8o;lUJ->1?WKtl7SWhd6g0%uEE8Aiq27>v01eewJ3Us|gR@?5
zE?qd0f-s@+MwdM^Am95=_gtY6h-cd*J+4SXsN*=vqp@)?4cFN@aX})uD~yEB85f|F
za&e%gXb}3`81`sxPjA$0tz9wfbS7?fTvIByDHE#>+&OD*l!cPo=rQ%9W5F1)-^$86
z979hH9khK@D0)90lz!ypHk6#--YsDGR@@Jdxy2Y755HxGb?uz9u-9)$Y4N6g7$G!x
z^a5WmT$wbz_+VEnY)w0YH7BBByvZ{AdKX_vhJ*%3>iQ#Q#f%Bv7w$s*plylnBDuJA
zB>9noaWodBgtJB^WI*n+w&ZE;OpMuP7?l{Cgv-KDliSIp<9&mjk!VFUgzY4H);a}2
z_LukSC(J{U@oLQO1>+9FN2Z@;(ukdyS8F!2d}k)ix)!R9vx>q9@0Am+rY%G7Xw9^#
zvYT*z&%QTita$9`-zRm!m?R7;I%4l%8-$?ItKyDjgknYcl==%6@tFR~c$nC>TzHhE
z_wB+zg3t%!@0oYZgK3hdZq?xo?8@xU@o~+7gYUi>It7V%?5cV=WkC{73SWJ;d~-0^
zVxb|oe2Dw1l(swY#36jFY|rTwyd48)_>DiF9|vuDwY_`#WD)U=%}-BmNXO%ClTuZO
zZpYXj%!haR@yJ_fUdJ4}1fPz47_nr<L7b0{PVVcN3w8A~7dMSh#>0TR=-r2Nu+`w%
zg8bSH?0ue>qkc079fx`sU&%>>>=b|Zy`@QLU%9z>NV_C-Ik%}#mt%1l+S|E(Wnl*L
z$7nz8G~N$&*0l{EXGWu=;nio#*W(fPF;HpYs$6J?oIPw)d<qGw%fyTMYY~1*wXyGn
zaJc)%>s->z#Q6cMTvv~Y!qUgHddd_YK;ogCei|#X5SYXK)bBwW_E}$e?q8CDb8#Ev
zY$qn8YR8DhDpj#CH1PAvj|{<#_aUc@okB1yD{Q}T?<mM?@4V7r9stRgV@6jGJV?ad
zG?Zctvf+ENx1Nl1ESSA6m5zOzjsq^+yUW$(Ah@yAMSFP>j4MCfSDqM)IUN@n%n;8<
z&~`t=tk-!kFEUroch1J)y%W1I52t`PZc(p{-~=etkGpx|MG}r^OqpKsaxZea$=Y0F
zN8#a8<LEmD`!T=!cJB)t6VO5Nwr+X&VJI)Es>%(Hgwyui+Q&!3v1M)J-7cF_;JQP&
zn-4c1j2u(1trHT_aj$uoIO%lMi_M*S&?gp>N$Mx2_D@C6DVGzs_+-FUU(|2Vs9boh
z==18KK{{ljjGu4hW<#fM)hOMAF_7Q=dPL>+D0p-ziJcUXjc5M0LyNp)k$bN~qT_&6
z>=Jcem$Ez#TJ8!CH6ssU?X+-}Mb#ddo{}XeUSy5;uKiE?Jd4Jan@@iU=cGd9A#sCC
zBEj1qE8WlkFkYI9{xV5y7v>Let$4O_7s3m=jJ$k33n>vNOk|{rpjm2lep~N6yw!TB
zUUwz~lZ2%-R<BRPXbp`q+Nbis4v-qiu@8fFy344^G3gj;Gep&)BoSWcw7I=mh2Z!r
zuN8Y)i1z+Akw*=;VDhya_jG+M@#_@33o{K8aCr}Z`=w{z828@gmfP_#OnK<<CH5u?
z{d0G}loHQ@i<po4uEUWq>8LvK`fGnYy{9)?ZG1Y!b0$VlTy+?ePh4n!sCybVKe~`V
zDBuu=?bB6oEhp{+d<Ga-W#Jw7ZdheZ4lM0_pN+hp0?F_LK@O*Lz+S&;Vy1I0#&4NY
zl$9HfNzW$N)eH~AsFS}Y`7BPwd7YePPvxAjbgu4XvCc%l*Vo^mV(o?TOVy-AI&VR6
zK$WClKpflxb<B)+Bq7J9-}xtE3AkOi`u?I@`_Mbt+G4%NKE$TYJnLnV5AnB==c8C9
zXcyr<B4I-=L_Tcyo!r|G-QI_!>BQvYQQ6EB^CGfvZC6MSD~(*txR-sbw&z+rs@opX
zQz9L&8&5nsK6yXXIh~Boy$VKk`h*6PqwC?NcVJBN;sn^Gqz=i`j6|aqL+(v%I<)S_
zIn5fAfp;Piyy?Z!c-JB1+<IYO+>hwT_Kt{xvoFu6^MV6V5q0>`V3Lmb9exV0)8nxy
z%J1AU&jYA?J9~VX`+98LKEP?eTRLX$+}b~TcqTMY>3v*iXaR$xF##PO5dF{WLBiTI
z;jnsLS}8Fy70L<`*>`p<!l2{NOA<~+z~GhvFWO`Sp4Kwf9XRcenHM8O#;!?+FVBCP
zqUL(ssgPAH%2I>go4vl}@p};U=Jb^gmkIm6{(gJT!6ckYE}GM0OFV|_cbY$HR2q)I
z6^bcM%frak7Xpu-O@+(BY~`$38895%Wo=aX0X!9*BR00{9?bCD(phHxemGg&$@Ub@
zgyM;*I;&3;{e8ZuwQ67(W*$*pDmQaAW;a&y-Fj_6K|tCGxtqCoxavyJ4uqZiD5Nu&
zGSl(Uer|jqG2Tv<(MdShn2JFsQ(Tq>r6D$Roa^<q1vnhE@Lh#z0>VY!=f^~);T+5F
z_O<zGIIXz2<FVjmC>%veiBB4O^dAwQCK3Wy|4FMKuQ`ZC`{!krL~cdc<-@-oHF3rC
z%?}?M<?n!onzxI1!cGkJjC{#bi-S(PTeB}HW?)=l(G$y2`B-?o?_#GXd3bf#q`QdZ
z0VKbD)_a&rHgXgCJA~dl0@Do-9M_cO!8bC@Vs3dBq=&oPPV&q|knwT_>8Z)^WvnT>
zHYFK%64pw-dA=L%Me;OXiznl9Xu;#HW3%9T=SJU}EvZ<z_}yLE2a))A=9+G8N<7Sh
z21eOm4MY5uZrgN=-9cPq*ui%pxE=OLW{T25Y?NI#T7*>$Q&z|#annS+SZ6rMc1<ca
z+7;H`P|t<m8r0s`%17l?rDvx?4?^jIo>eL_er#1<*NIh-0IMF16H@}SaA2Z!W%RfN
z$WC)#<#BH<dar&Kv}?{GB3`F>E#N^4#%~(<s$hK_HpJcBC^~!_B0M*8;!O@A$>ZXF
zTct$EUcJ@s;oC^Oo1{2o=*vuOOurZy>=uJ@zj`hCI3XW9HfZ)3;Z}&Wl91p*$vF@&
zj8%ynQ2-UweN2`50tC4qU63P}hd`~P*C!2*Mc5}{pNp5Wu(5~t2re%PYtQexdvpuY
zU(Q!_G7}5Ni;F`O<>$)4I`0N!Rzf01L=SF%VL?1ri3U`AnC4;qz1UUV#$_O}Y{1>?
zhtiNRsP}5l_ycHn!LD|SnIm*njbx{5hrw>dmF3==$05^WqB7Pe!*j-)u}VEpLeIUQ
z+!6T<96H-jz0Elv%8`SWmOV<v1^2Yb*Jb<QxIryscu+2KhP<heY)D7m3GvE{4cV|*
z(m|=xJrb@_(U&gwJAi`9^3$)5x?+Lpu&6BYNUWN9*}!*57S{L<Gbsu<fFk#vvK<~p
zqal62&g5`g2<h%Wy<*00>>!4{v!am%E_bIsu}X(%g7|>z>of6IHrv>KULjN-#P-_M
zGYzf>HboY7Cj8pUTYgLDCSs|Su$rWLHk`Ufhw}*kd_U&2`iatL80{LNXBKoA!dIkb
zB|p!EzFJ1=*uAMRHnX-k|0)?4+MLbT2WBFr;wtCk@Po)5KG<>J%)M~-YV`1mTZoGl
z!c#3pgE77T@z=ezqM(=<x6avd2d+GPzoJXuBbayZM)7-oHlD7{EKMTjKk7*pzXoRK
zLe<89+PnN?Xx}xaUnk`>C|L)VJvYq5sD><$U6-;EQ}k@s0Q-DQd|7())Tmh4Bxc`q
zJDr8RZ2|dpL();+(Lz!Bawcl(vcjKdN8_^kb0w`F=}<KvvinBu0hGUhgzd^ur0HF%
zRh=3L2d44e0?9+T^J$Ht4l^G~>+gTiy;exfqr#r7CfcR&(Z8URc^=x|_1R=Gyck=%
z1uiux&Lrl!Z}@Trad@GnqrEaO6Qk8qp1zMwgKxLv44H@+T$8aZw`PZ6sl+<FNw3my
zv2u0y*)uod+)-`LlZOF#8lSrA(AafQo1<Pji;;mD)>0FC-AYGpUqc&(8++k>VHu+^
zG!6FCu6*p47LDGyJqD<KOd;Zrf**rHujA=Xt<Zk6nS271+D6jvtqY6ZPRh~iQOI1|
z|9=I6=Ig?AD8<z%`=N30a3((JMr2&vn8XMD1YhQBS|-MDw<o?}_a)09ze~qpG=1RD
z(b9?dex4)K$>Z}nGI~AJALHkwCDZ#X=yhi_{r^}8wv*JMxrGPNs<jGAt_y3&P@`H?
z6W(&O=zR=jVJ*wic;63fSx%ZNezT)xIU2>FX7<tYS^r$=KWA60pTe&GIXr593OxQf
zOZ>mRl9=b7_dV=}YAmw#dOD@+U5qJLk9~SfMx(fKe(8S06JRjzhdOl<CTQEyhCmwv
hZ3wg>(1t)80{<2S{<(SG{hwl9_s`)`|5M=c{{Z}(|IYvb

literal 0
HcmV?d00001

diff --git a/test/regression/basic_weight_windows/input.py b/test/regression/basic_weight_windows/input.py
new file mode 100644
index 000000000..742934594
--- /dev/null
+++ b/test/regression/basic_weight_windows/input.py
@@ -0,0 +1,63 @@
+import mcdc
+from mcdc.constant import INF
+import numpy as np
+import os
+
+os.environ["MCDC_LIB"] = "../mcdc-regression_test_data/"
+
+# Uses Pincell model as the basis
+
+# Material
+fuel = mcdc.Material(
+    nuclide_composition={
+        "U235": 0.0001654509603995036,
+        "U238": 0.022801089905717036,
+        "O16": 0.04593308173223308,
+    }
+)
+moderator = mcdc.Material(
+    nuclide_composition={
+        "H1": 0.05129627050184732,
+        "O16": 0.024622209840886707,
+        "B10": 4.103701640147785e-05,
+    }
+)
+
+# Geometry
+cylinder = mcdc.Surface.CylinderZ(radius=0.45720)
+pitch = 1.25984
+x0 = mcdc.Surface.PlaneX(x=-pitch / 2, boundary_condition="reflective")
+x1 = mcdc.Surface.PlaneX(x=pitch / 2, boundary_condition="reflective")
+y0 = mcdc.Surface.PlaneY(y=-pitch / 2, boundary_condition="reflective")
+y1 = mcdc.Surface.PlaneY(y=pitch / 2, boundary_condition="reflective")
+#
+mcdc.Cell(-cylinder, fill=fuel)
+mcdc.Cell(+x0 & -x1 & +y0 & -y1 & +cylinder, fill=moderator)
+
+# Source
+mcdc.Source(position=[0.0, 0.0, 0.0], isotropic=True, time=0.0, energy=14.1e6)
+
+# Setting
+mcdc.settings.N_particle = 20
+mcdc.settings.N_batch = 2
+mcdc.settings.time_boundary = 1.0
+mcdc.settings.active_bank_buffer = 1000
+
+# Mesh
+Nx, Ny = 20, 20
+x0, y0 = -pitch/2, -pitch/2
+dx, dy = pitch / Nx, pitch / Ny
+mesh = mcdc.MeshUniform(x=(x0, dx, Nx), y=(y0, dy, Ny))
+
+# Tally
+mcdc.Tally(mesh=mesh, scores=["flux"])
+
+# Weight windows
+ww_array = np.ones((1, 20, 20, 1, 3))
+# Actual bounds are set to arbitrary numbers
+ww_array[..., 0] = 0.4
+ww_array[..., 1] = 0.7
+ww_array[..., 2] = 0.9
+mcdc.simulation.weight_windows(ww_array, mesh=mesh)
+
+mcdc.run()

From a45e858cfa490216f25afb30d74c3adb624244a5 Mon Sep 17 00:00:00 2001
From: Nathan Glaser <glaserna@oregonstate.edu>
Date: Wed, 15 Apr 2026 20:28:36 -0700
Subject: [PATCH 35/53] back in black

---
 mcdc/code_factory/numba_objects_generator.py  | 8 ++++----
 mcdc/transport/simulation.py                  | 2 +-
 test/regression/basic_weight_windows/input.py | 2 +-
 3 files changed, 6 insertions(+), 6 deletions(-)

diff --git a/mcdc/code_factory/numba_objects_generator.py b/mcdc/code_factory/numba_objects_generator.py
index b03e05eb4..c466b92d6 100644
--- a/mcdc/code_factory/numba_objects_generator.py
+++ b/mcdc/code_factory/numba_objects_generator.py
@@ -1165,14 +1165,14 @@ def _accessor_3d_element(object_name, attribute_name, stride_2, stride_3, setter
     return text
 
 
-def _accessor_4d_element(object_name, attribute_name, stride_2, stride_3, stride_4, setter=False):
+def _accessor_4d_element(
+    object_name, attribute_name, stride_2, stride_3, stride_4, setter=False
+):
     text = f"@njit\n"
     if setter:
         text += f"def {attribute_name}(index_1, index_2, index_3, index_4, {object_name}, data, value):\n"
     else:
-        text += (
-            f"def {attribute_name}(index_1, index_2, index_3, index_4, {object_name}, data):\n"
-        )
+        text += f"def {attribute_name}(index_1, index_2, index_3, index_4, {object_name}, data):\n"
     text += f'    offset = {object_name}["{attribute_name}_offset"]\n'
     text += f'    stride_2 = {object_name}["{stride_2}"]\n'
     text += f'    stride_3 = {object_name}["{stride_3}"]\n'
diff --git a/mcdc/transport/simulation.py b/mcdc/transport/simulation.py
index 8807977d1..0ddb2bd7e 100644
--- a/mcdc/transport/simulation.py
+++ b/mcdc/transport/simulation.py
@@ -352,7 +352,7 @@ def step_particle(particle_container, program, data):
         # Weight windows
         if simulation["weight_windows"]["active"]:
             technique.weight_windows(particle_container, program, data)
-        
+
         # Weight roulette
         else:
             technique.weight_windows(particle_container, simulation, data)
diff --git a/test/regression/basic_weight_windows/input.py b/test/regression/basic_weight_windows/input.py
index 742934594..c730df1bc 100644
--- a/test/regression/basic_weight_windows/input.py
+++ b/test/regression/basic_weight_windows/input.py
@@ -45,7 +45,7 @@
 
 # Mesh
 Nx, Ny = 20, 20
-x0, y0 = -pitch/2, -pitch/2
+x0, y0 = -pitch / 2, -pitch / 2
 dx, dy = pitch / Nx, pitch / Ny
 mesh = mcdc.MeshUniform(x=(x0, dx, Nx), y=(y0, dy, Ny))
 

From c3b70c7b690da16eda51a094aa8542607b1070f4 Mon Sep 17 00:00:00 2001
From: Nathan Glaser <glaserna@oregonstate.edu>
Date: Wed, 15 Apr 2026 20:36:25 -0700
Subject: [PATCH 36/53] change lower bound and add comment explanations for
 choices

---
 .../regression/basic_weight_windows/answer.h5 | Bin 39920 -> 39920 bytes
 test/regression/basic_weight_windows/input.py |   6 +++---
 2 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/test/regression/basic_weight_windows/answer.h5 b/test/regression/basic_weight_windows/answer.h5
index 93163bbc349c8468b636d0816d868c413af83ed4..32f66b1dde1dfef118382d2539d8c4972572b301 100644
GIT binary patch
delta 6520
zcmWNVd0Y)`9K|J)WJ?I8ENMlOB;{0+5K?HNP+C<)d+BOlR9bJhec$&L=}n80D3z9b
z=PpIjB4jBfdgsrX&wOT{d4A`7&zYHixOn>E;__~5iZE|(z}-s>yc$0GacF!H8zL{S
z(RXA+KwL&MPm(~y<(YG5;+arm=<MP#7=lWMO*r)^6W^;wv;tKJQN7Esi)PM1dA+Hi
zkURs~Av{}>q8PZ+wL4mtdl>R!K^j5aLy-S-%xl$A4jP|-Zy9HDps!cDD?^G+LGUq+
zqLWT6Y;|CEd{JkEX20b?Q63As8&7xKOJt#C>j~R?`b@}(@|@OD9L9mS|30_}v2p8f
zvhvjd7T%V*-gbGxLU!AZa~a<VXd0ib-KWpMv-1x6GP8paXH#5zhUsXaO%{B#ArN>@
ze4%qA6Amh~l@jt40(E<C^KfQZ&<t8vCvNc;meiqYvkn%%h1Q9u?B-y*is|r|{2_Qt
zT{H=fXF*AC+v=U<0%LoO#ab^CSYoHWl>Uu{`ZW8e489@QHLq5FuQY_HY}~tZl8u^f
zeVN8(19+AbEA`Bsfu+x^RL)EfLf{jZu@YrB9V~XXyLLJO8&RIKj$tgg^l=+}C}hD-
zy0>sYc>}(P@Z#k^hv1b+k6gIPL554zlxjN@TOE%US-7$A;YgnJ7kegLM%r9xs!aTd
zv~$<cB4GDNZ>j2G0^QdX-RTA#Oyn)m^XDCcQrokwDStTl+1<5+vilGR;+5Ywx;*Eg
zRzsF*IL(Co?#}DN1QXND!!xED1ZKx&x?eUjVY0hc)l`QA<>ww3774OoyKeWe!Cof#
z_}*PITr-T}@mlW-@0oaK+*Gfb!@`aV&BVE_9Qc@xt^UeOKtcM_{+kjU?6|s*P5(WN
zwHASt0NP^~)IJBQX})8D!8-oSrJsZI;$N0Y*f4PCzR)v*pFpD4^Hir#1d@Y}FHb+i
zLfTNclCB_u)bld}Ei(jocoVWuo*<z8<Lht#n*=UboeN05Fo1FAE3TW*u`nzwKH>0&
ziM~ze-VOG!;K99D_I(%u-pdqKHdB%bmGG+7&#DN-pS-1fmC8Y=m&5M4H!SeJq*}{{
zGw|70t6({Ug=K#pZ2fqB2=(pRZ7)*@w3h}da*+$`TJCo3=rsmHc}jV%q>$npkr$}E
zN{(PtwquPI14;2B9hJs(?AAIt_q3Ub9Tst&@+wT!QTXz23EySH)_QeBloA8U+sFIr
zHnBjVo6G$RVB%Hv%O0yI1hO*g&2)I#uym(XpM1bV<H}!<cOOL7b3e01e+H3#&GlgO
zM>dWe+3Wqcj})~=iRIPC0hm4iF8c3NA6%7VnvAb75w0KbgnFNj-pb133G11lL|jYB
zJ-&~D{4$3!ZyP#nugcWxeJ60yY(L+>%>?|Z1(>qzhgzIxY1IN9$|CCfDx??){IIPA
zAn+r8(uZwDz)>l>Ws@@t)9%8}U9xNhM~>d`BWE3WL$OPb+~L~}Gtb<>24OmuB1<D+
zy}DIeWh0A%t0vW(eDVmmTjmEl@iI_p*+V_CnSrr^y$c#ENf9<ETK+gZh>)OPQHK%<
zc*lp|J}JY%_1<~0h^an!SqZ<^T~A=@*r?gWJvv5%xSU-=82FUxbE=;?fZooV3j8F%
zmKY4M=Vh3<)L*yHZJI#!YZvxEJIWAf3=8Q?6{MXuzu*fYSDaxtb8NtU7}`p1-qrD>
z+AU2fy@^cxY|)=x%EiI+;^C#sE)ei@<*REVM-(_Hb8~e%6AC*HTKfnSkT3njy>SN<
zT{cDbS%GxW>&~0L+B^vVB`wrZHUqO2r`ApT(Loz=sGO(lV&mt3FPCS(XF>3qFlBfn
z3szjmeze#VSaeHMc0i4Z4lUo_i4Fp1IAvO)H`p-94LC1Y&&J~e?-zYPPRHAeb^5)e
zTegOt*erXTzz54xl~oM{sACTDVI-)+N<)q&hY&DdWPi2HlL?8(3SU+EInex0eL>kb
zMBul;7FKo)6Uina<I@66OtybJdW@fewzp=!qTvLBw-mk;;bubn1>fySR|2dj;iAC;
z1oDRtY!m39LnTZ;elaQ1H|l5VO?cRlP_((>)WbxWTZxUwN)Do)zA*16(fp40(bv3O
zCbUd%`ssgU;?ce;%ApmtEST&I$$C9PUhqbvkRro@)7$1B?=uJtt*FfIvtglw*1D<n
zEfZHVt-2C-aiFkyXN%t$1Gq`QB6o%f?pKe-4v4X_GpF}l_I4&>BMgg<g)wpRb^fZH
zuPo@iGK{5NVd2tdAr0X-L)goF8M=HA2aG?Il{bPFI9OJ^Tf)hj1KYA&pH{dJflqn%
z&Yw>lRIRemuzA2nL-e3c50k)i>gFAz0}M#^Bq-0eGcaABFFmMEN77^Q8(wM*c-2Pk
z7rjrQv@A_V#E6ND%l+c2i<xj1D-Jn&bqEm-)W3r>Of1lPM&{cX(5Cp?sHU_qL0!d5
zv5sLw+hC*6jesF!uD6X6zs12HOTlEjFT<E#@qNt$uVLJ2G6+#r<RGHuu5m~`3%$9w
zgiZQci1u5V+)%?pe%_zBzEk9mhW&^-899uwGq=w*hjC!UYMMwlV}rX`O04x76T0Ft
zE{+^FC<B*yUlo#sGt~J{MZt@K?X!BK<8%TgZvWO7k>_>WV=TnUc&=Uan(q3EiSDW!
z59STou;nk8+SxRO+po`0Y2RbPu>8YVejPaqbCt?h6KuHW{L=fZK8!-OM>WsZad6Vl
z<lGJ^4o)5Lx1CzX!Y4T``^pn!{?9q{H5{BFQ&8-<xsMSG2cJfDHhDAAAo;hWID)_`
zZ?VEa6E<Z0#&~-T7}(?3Jzv4kf_3veXT2GL*KcB4+ewkCyUPrejL?Cs^|39aOHN#C
zvuaReV^Q6LqR|!>s!iwKD~hpDs_mrXBuw^&TXH#8Q4}WZl)m(U$-=1rW4kdb9qF89
zv|l9Yg6;h7>w6R6N)=4pZp}nFcluje&=7b%uB#4Nv7y{*KU2Aubb}hRO4XBrONw7r
z&XQnhl$$lM>L*~>oAF(<lmqjnIc*<Ek(bMaItK?2AX-h`u%ep-IleQhi4;o~l<T(c
zPz4K|<}kN95jKi7C4IW0nW!9GbHtQJM`@p@<DGu;%@5d!Gf0w&DUVGAZXj^sx`Dj7
z3j?;9PC92Q2y{e7du<scwcu94U))K$r}#^E8iNTplUxU5ve3-Ed7;B1w_cQeLvqHB
zM7*I~6XhqFM(SJ6A6_Ctz`jUXsg+c>pv%FRmbxrlPky4bb`62@i2Xhrt_+}WU9(G)
z858Wa#;)ZPbfjL)jk>jn1==U#$V)P2-_=TAGM!{#vDurp@yml4mfC+{tviA5mpm*C
zJD4aKdH61O4cQmMr_7aClT~bEDTS9qzWJ?xGbRRJEG#DOhn*ECaQm*pqM|(n(tiw@
zMUZIL<B<uxN`m9Aev!{_1=40GUnl>%OrY7<FZ2%&11Ie^V<`_=sT2%$USjoQApcZt
z?MWt#&y=3=kz-@ztI2W2WdyQm&T=UvM|2Dnc?7cv2;@?R*P8e-5V>u@ER!q<b2YaQ
z-rr6hC$_%ZPTGuJv`z0lnQqY+t?3+VQVXuEwH{yT$fWKaJU&2R>wd?-X<>9&r;L6)
zqD@DY=RoDycLsjGVN*nO2Vq*IFm2^RM|*((=K?kfxP&A26iqs)o}o$W9SQ8DWE@!$
zL!$A+_aBu;txP=3U8lWxB@26ouBh6OA(L>#FkDZUiMHKVOGmGhJ)!x@x@k`q&V*DL
z?3y4emGPNV*0SX6PelG~+($0VvSQt%kAqlh`1Sp~5J|uZVz42a-0I|?dOxY;14HKn
z?V^|%-7s^-ErtN)n}CS>V{r!J+7~YhC1bR>&evK#X9(45a{9YtSlGUDG`o*v{j7=F
zx0ZJVnqH`$y+Lw9bXKOjfHD9br^^#6A#}Ve%sF93=6v9i?}62%&E|H}4$jFju`=Pu
zytoe?VXmEi#Tsnr1(#^Zvls~Y(HXfQN(iB@>*1=N0s`mSn|xfxNw!DxlrHrgt~%??
z#Y@><$6k5DF$21H9;mjh^T5SBG9r2@DS+c(>=oR@+Ew3*D~pQpqIO=<Wmg;q2JT*4
z77~rN?drR>DLG=U!R)rz@pQOyNv&{`sD|lr1N9c=7`UCRPTQayhd-T+(aTl^5bn6$
zb83Grxb7PM@#m=kTP85k$CE~Z?L<pJh;R<JZhvTFvpk*TvsmV%&@7bGw*Sszmq09g
z=dyg&T7(Qm_dRro#YfSh)k%rDs63p$t8jG+{F5I_OlxPrd?sJ{zwvN%_lQwLSA?LL
zn#Omy!wrQixRk3^qT!$!#OpDcj%$pRZpmMH7+TW%ds3VdivWhB?~w}y;8-Yrsjf~U
zKYx-cTOZKyz5e%qb?UbvC8C`Ciq981GqnM2gZAKWgnUOu0MbL)F<woXkn6brM>59?
zRCNm%S7tIa{hRuJ4izBe+m}Vsl!sVdS7_I27zXJfi#x&hO;E|0ek-6$h5pnfxp~UL
zT9TJFm(BdnBXqTW!D)kByqWe~X8Sh{DW`J3K7RfHvPZ@ak2W1aCmhdmmLy}%%-#GI
z2Xdjo^=N@f%SO@8$RyuGw&=gQK?A1}Atipn<Jz=62^CIDcX|v?s_*#Ru9kvn#sw$!
zi3E(VG&hj<%)v;5ZwF<alM~pwUlvkqa?r`^LR~bNiQoelh9l>*AkZG3&#xMVh#<Yf
zjV3wx`ugg0c6cy^7bldBn8m?1BHC4q;(`|Q<-fH={bADCyRG129PZ|iu6n<s9E-WS
zKm0R~Lo5IF0PeS$xTEQFQFv7q_AfhqFt?tD{e5dFk9yyip?M&@wqa`#qUwl>x4}7B
zL36!%IW-3ghs!SM?n{DX{yDGoG71*|J3L>n90?D`jo;GoPPo(X+~>jWXv}wPN(}s$
z1&OhtAk$83oEp#%+UuT$6`q5hW$$dUr;@(&gwtW{m3x@9>7GAYLl(X)X{AB=3*~EY
zfL<l0ex8mA5KIR{=H%2>=NN=&>@a!DNkIQ(-^ip84I3RwSL$k%!RA@t-M6*zU>^^x
z`=?h1E&s%=6wO-*&s=Zl<9QEVf_JyvYjVTA0k5KQ(Ikut2;Yz@rs3}VxSM3Q9~Ldz
z(kM)cMi0N~wpZEFU{e<Q17nJ!;C^P;<pH`K^uCYmt)ZpCXI4n?r-BFC_YX*FG*m*r
zkmY$kE*qaIhdMr<j7R)neB+&WR6L$BO<p(Xfv_j^=7|<h%>0d4?7EqPeN956V&ZwY
z*n4K=o9HF9$hH=i+4|zMDE+VQyNBrPpssMXC<i6>D<#`2F&F_NHw?AI&9GOm-MLQC
z6ZjIg?XY|jHtt^?`Eh$H=9HWYrcMOlpuMr%o~%5`s(v2I`*R-@`zrhY;-kRKyW;DU
zS_5wBXv^04a@3dAtyB{vN8{AQdC_(Z%^PQZ`s~YbS#;3Xj2?+-_4ntDjLS*?uNLz-
zyC(^!0uKEY;7>q2Z&_hLrWAOmg2iOKVt{DhwhV#$h?h0C4C$vK{ZQV9w+m574d!OJ
z@`b`uFu!M&bPBA_{BDsr6^o^Zc_x=jRU@$G*Obvh9+I|y@QGMlfb*|K+uq->1E)iJ
zlfWfkeC!C)IscZD4z0V_@0r+=BBZ~V{1@zq$@;g$7B$5P^Sqz=@$wZIeu-Rl?^7<m
zjBM*Zz9bL(IP?USA4zBv$p2y0m5xQl`Sx<+8Q_016dT_k52}Wm)$EKvY-jwrqL2At
z>Gd^Td3uq!WA4Yp=j8+Wm-+XO?2Sf$N`=%g#V8$R6OA87TMJ>)xGnm8S}yXPVh#Pa
zhU2+*eZS(7a+Jjg`An^fhtS#@jdcp4IApGDx1}Zu-vZjkzJJp|mjCzf`K(<0$yU*Q
zaXu6KGrPmK>`Dj!;LU==JG1fO93<@B3-Q!{Lz3YMU%ZHYQNWptz^>o<iIl)Ze~iZl
zMkF_tLFvJcuQxXqL#NoU;Dn4jf`zK7+|iY2d%E2EZFM%fL!%~xH4|{@U&XIN_gwN%
zRl@l4wlwfJeBhW(rJ>QosMhvo4h*GAA7uEHBDIBE{VXR3_6J|*THL<}n=4<^xOc=N
z^R$!T`r}D(2;oq66`};(OF1EqIu0mUHeFzPD;%@;Ldv~2m&4<4&3^?_IS{trF;dU;
z!IuIDxvx?Y@ZKv+d)`upP~$|Ni75}nmdq8e`BefQr!(Gulc6a3!)u?DRt87*<+ElY
z<v6J9Wx3--Ce(@yURF$`VV0V`a!#NUxfYbO0{p(!FgH$(Q{gLsh@we=kWDyFz3<%E
zXGz2A6p5D`#&VGvb6E64vp?3VWmy%jbw<Ue?4qqP(I|}C`|#+tLWl^+vCTGRBk#13
zuiMWe{J8t!@yN?aTwH8xy*)?`BHn+rpO2=)Zt$s(z_V0X>O7?7_h*q`Cl41kU4LQ+
z`zW5a4U=(rP5B(THlr5)zxIhPoNB<e<?F>|mK393JA>NJa>co)lU%JkLUCHSqK8``
z6~<i8uR0u20X0U`BUZBnRPTgpnHBj6ZB$VAr{qIBWbkeY-4k^ZN{g-Xg5jF8SdFEf
zkHMH>sclWApq$;Lt4gb)p@+S<V1r{Miu8-*4O`2RE>NAL(US`23u9a#-D*(ftj_uo
z^a#&y1jveUGGT4@BB}amA*2+}W=%GvqMci^xW&)|21l<=_nSn4_j$74I*(lJ7t!z9
zq#KR$9Lb{E(qcrUWoO<w6auZ48gZw#-=jc0EhLVsC<Bs3%(m(D2blCow$I*{iCaYn
zCR^e&F{|3BB}*@XxKWpUPnJ0rD59qRF*y(`xHHntT@LOK*G)7wKE`Cjt8Y*3ic!+8
z@2z802@MTz$zJs+sA=r<@g-d*#JFUC@qGr$_jI|GFVDh%&UXWaPEkVeV|CfW<M13b
zo7>FCa~H$6&Ts6OU?w(t$97LFsezD<YfXM_1v*lBYgO0ZhX1&C__8i4`lv<S{Ac5k
zwclf#V@x^Z6E~TjaSBC{ZF+Dy$rTYhfptaiOK{cPg{zoai?fnJx^Eh*v2t|a0AftA
zG<A>mD$4pua0x0(uc#r}-c~ecq~HmD5qIq;lee(uZo=kHEgEL^CdI1U9q{w%)zL}!
zWGF=bU3Ag50C%P07HX<}v1V`$FTZOh%tEck-#qk#!p8L*ciUaV^0@iN(GnkQeL&ef
zmXHC(Qqi&-Je3I0^33{ns2GvIc3V(X#>>%QLDgr==0W~bY`Jk$EcV{Y8%_J?PX5@T
z{fiLG0AImcxfXJEXAdg|=K4lsY?FeEq<aO7I#N#Y|L{YZ-@&o!!em6!UwX<n`(m!O
z*F3;E8aI#J5!q>!0oh&*?0Xmm;Wp=^e)=>7cl>)M$d?2O{ezT}{s}kieW%<UIh=x5
zTMSJeo%h6vu0tE7rTtNV+FC12G!MU$!lvXpGck8#cHKV5M3^NsMKs0aBPvH-Pb4xL
z`?(G$YRQK}Lh`hJJ6{wgrTsWX%E9=UJNZp(cN{)PJT5VQlY+-N)>ChSYLIy>V{KAU
z7Um96c(xTLI>2y5MqT$%DtOrUtxN^8QTr%D==$b(_{;O$EUb>j?HSRh8h%=kPHVTk
zFc*tM-WpH#i~2%qRl{Bjk#MXNMUmf&WW322U7PU81sZ?tBxs)^ARINET`X6K)TCbJ
zC+g*p_<MtuQ<j3WGb&242HqH<%r!;z+ZaL6Nz~d*i7X;n(aOIVwaBe2K5<}aEa(!!
zLQS4A7@3@qtzH>{KJRVSaeY2`_VZQVvxy*VZ8-Aj@q<$6bQW?cX(Yn$Q@CNs%{&N=
zo@Ab0mjROLY?;L*kOr<lu`r;Z#klORq*xxd9M~Gv<CzM|vf1iWddI1Fq-o{NR&fAh
zdB*Y|<>fd#`7+!&=`jp<3&uqqEyZPC1^yQOLcFP6V>c6)j`a<d?+Bl4e9XGvzC|Vw
z`Hr)@{@Q1wKO$mh#hyG+5A07IT2cjeXQM&2c{;ilbjq~Hiy^J&tZdp<g+beXqYZ!Z
zg)qN#vbfDN9aqj9QH*Bt@T9Z*ireDqs@EZ0ypqqa{P*u2znI!%$?}kR;ytW;W0#rN
g3#$!_-)yc?e@6Bx+4uX1rKoD+r7x@MP7quCKRpDOI{*Lx

delta 6520
zcmWlbcOcd8AH}UoB@H7*DM?0UMKl~yMj>CjC`zeRMj|1XWMoHV?<@1V)-}s05=AKP
zG%6_)N=wo2^ZWDj&+|N=_xqgJIWsjpGc`QaR~Je3?Po7*oMd6Mka2BpAsy;c1Ao2f
zOcb`n2lQeXN>9X$=O$T5WbRzd8fIeM-uG&2R&%hzr6Bv#=n(Q-Rtvf+v0(OWd-HBj
z4)%Qy;QKf~2=OODuXAoNVBv7Jd9@22qJ5GwUS4#_-$*u4m1ki0hKh~q9t;v5^X0e)
zFo&>dRb2GZOKc3)E0y^SkKi!7sjK`a2M_E1)6lPDqpCxwTkjteD@(WV|25~Ju}3Dy
zu#SU?(nUYory)GosZbVq$HvrNM}D(vHiUxdi(Z*SII$~#Ti#_Bn7#d-*Bv?VQ>fms
z`tJbP&O7_O?~w+vr>9mg<|Q3J&sVIK?c`w1RuAiZCkDbtW1<Ji#DSNt(%*z}@Ko*F
zU#~V6z8CG+>2qTuS6ct2PQVac?w(V=JI;bYqMCI0ayIVz3T<qZ9zf^B_hX5NSa_Xh
z!n`%i2EW6rjr3C-bX@za6*Wdjsd2NVImw)fH_j6+-DymSRi4}>vXTYEWWCG(|F+Sg
znwc`0*dix#T5JCx##;=ArL5Tq5&coD@Qr~G?x_=GH#Vp*-4xHWSn#pjca71)LbPG$
zeAy2MKF0b!^lM_nF3086(lHKvR(bl05f4zdev<NCl?83K>M&`9_|y5VSId)6vrzJt
zme4ZEfyUj`83PFxhQ6+)<=8OMZP|LUSeAh>L#l!7)BqB*kCW|;n5gwS<~nJ|!RhSn
znsHt%ENoPL;hab8_WjQAgK;c~1lHT_izX%@@}5#COh-$Z(p&+P0UjUN6vi^})_o0W
zS&BUyH~2PN?6YE`>C?{q-#6G$q~~-5v~jQy*2*^)%m!cg;Zm)eLkNFHxqVK22w8JQ
z$rmeFh_&04P@KfT+_yQkr%epp=Dg({b!4FXl=yKzS0;G&^8N~Zz=H4}YC)RnFx+I+
zZ}<qX@I`5biA^6H%qr4s%Gwi5JTZ8BK;4@Ky$ks#h5VUlojb}YzQMv3dt?4bRvgr>
zS10!q58KH<cJ=s07RZWuHb>b*2(uLwYfBx%m~j8y>_iTfCbi_kZCD5y*nFs@o`KCW
z*2QW(EKrqqh;?-^&>w7+9=Mr-xEGrtS;mAa>BwEHXKU!->M1|NFT=qeUHALDKeLb#
zAfaJR?EYMzzMz>b2ih%{e>^K<<D>48ygyGENP0K2-NKZK0I3tZUj`22f{0G!#99uX
zowyeMMVJNNKF04*V-7x+xtKh2X2ai8VBnK33x5{L%dgF_QJY@S{%AcDq*<<a^T+tv
zP-vf_=A0#b%<K#b{mO!><2jneJ|=>^RTH1~GO_c<)A~N*05evS6z$s}!j!r5q;|8h
z($APOTFHV<%HW#S+HCYbS8@LJiqI12xob3WK}fyL)PpD%7G29sB?+!X?D3q^=VQa@
z(u!m&E)EGUb2mTt?`I-fP4VZif;Z?B&nKyEVd3PR;J-yjhj4bhPs_}e1LqRT^+!}T
z0^(ie*LJfZdyPBamC#0q?N(%@vhh`MuSUXS77m#o%xb(kjG5ySZ)JZEBec=~<kg5F
zgjAn*TqVQAfVaL|dkF_p6a`<ZB8eD{O!`F6B@U$dbQFpSlJ$Jx9^brt2vd8OW^xOS
zK<}r&#LC|cEFP*~cj70Z=9o|B9}<}OaJ%HLbu|a~Kc+8a7I5Hr$>B=yC<omStG1Gs
zaiFzjBq~jC1Xu4}j{ej+gb!Y`rf*V)v9Z)*>Ai;>Y<=|jNgL@si+G7UOIMPO%krl;
z#Hn-8ryC?IahL&F@&4gq85Vjio+wm@vr%(OG5(qFF!Y7U|2|u?p`Sya5ea4C@9S?0
z9b_i*mdyIZ-5$o<l`1v8zd494im0O#le}g9ODRQj7&kUp9V^^B0_XMr-aFQFAg^q^
zl606L+))vBmQL&lhT;niH&$|BCcL0KEjf(AuV(X+^Bl~d7$4a%L9pFo$BOR+Xtvhh
z_p;Dr!`fIms`oqx?b!=171ue)!z(V8vm<brFmdpe9YP0rOx1WL8?HNckOK%_#9fFH
zUhrnXC*)4~;Wh$fpDak^l*c2`csuu3%#@9WB1J=<tYQ4{7&2{MBo34w{VCGO#+Bt3
zvBy0)2&C{G)L{%`#Cp_s+oll+rZXe<S+ViOdrA9q2V#;}G*iXKnegBX7;oD^OvWia
z;(a_FIbADOrMA)W?90jTeeo=Kl^)Y@%N_<(gmjh6FUUsD!+9;~ULt%7Z<g#9;NV<p
zeEka@HeSgOHWlf!Ft&_u>L4MuGcOK)qU~fMb9t(SQXB^}wr6O%VMEw*?`n$snITZW
z_7v%^7=cgWwhv$aaq!k`Ve<HrA#7ZKYzxDmjq69B^|l|QqeH4lsEZ(F5=oiwXrvGm
z&Q<GhSc(P9?LFfsLF_y}jQ2nJVbEND`RM5|@roAQ#ebCplk8LvDFRd!H;&8g92Rse
zLWEB*BOa>U)PGr>4vAk@(T&^bINJN%u6cyWlp~4qR>TUo)Xb6PmX6?ujcs1(IR?1T
z+}PAo%m(RNxxGJInFXt8`K-haCMXf@)2~LE=s35dNuCH1SN61ztu6y+##Bn}2+<ri
ze9;hajYzw`)-PG<!^qI<Ee>~NB71lFDotVnX3S)vN+Kn>L~4A?E18fQ<^4HCVnN>3
zqgRWt-@Z|?Su=a$&rK3vzw%*{5dELM^bx`ZqMxQ;-SK3g!^6s|{R|sEpLWgMnPS4>
z^Bi}{=OGwRn0qHD5yUfTF!XMrBSY4BQ&T4ckG_+m1r_Lc)hVvqn9PF1*kZeaFcGsm
z-kvsI&Vj{}WUfsk9NeC_Q)(4pV-vai*p%B4OqLnn&)PuvUx~JOabh(aef3>dx4noE
zvi2K1o6p3&g~X-BuM8X=op|ri#ll2#f6%vbI=20h^yz3C#97N8#&QiNE^K|FtQ$s$
zjq;N~&rv#>GHogk`!FE%kz}eUMTdZHfwEsWlR$RwhXNMSh<55OQ3+$A`DgWme<YI;
zWU<t*u1y@q=}4p6D{*W@PaJBr{J$v4Ctp6L5i0(r^Xw6k+cGI0g>6BE@71rHloS&M
zrDbJr_I5UwNsgBl8ZwcvkbGz6ISW5hue>|%L5FJiHJbe}5n@eYmcQ>aF!=qe0Q=P&
z40&hnkN7bNS~^*jWK_n0w%F0UhEsGTYN;-!nJ~foBr&Bnjg2?4-i)R0!zke1vi0{Y
z2hrSfnmKP6SZ7xgd1;)DbR$;{Zi4M9a+cI^M><6Iz5a8pp9%H*mmgq`2^p30q!Wt_
zh$THZ7En)zf7Nfzj1~rRTw*ve=ZOL`qD*>`9>@Si$lE+)WDs!$#xb^<#DlgrHaZ9}
zARlA#%&LG8qx&&RX$TudqebDok_<?`>m679Isz$Oli3Rd%c&QZmNXIF@VKMb9G^A=
zX?z>RSN5?`)v*3#)k-?jW!Ag(j4`n9wzfDqj6sxK?c-`}1}ff?E|TTznNaSV-sWOJ
zFuf%u!#|Ra+Lwt=r-^)SNjrD$H31AKcZ2l1CkcBgs!kphWWZlxlhbG@177K+Yjl!W
zSW3&NPDma?oJ*;by)zv}Oa3-b>JA}tx6*+nM2KmJdA!+`$p)*jaaCIi2ln@97cMt*
zKqnok)6XSTA7HJuNTCg4bb>W?qK|<Fo0|uI5Jddwx)qA@M8{KCiLfKYxIy17v(R)9
zdG{rlw<DPNp{Acv?aG3LI3%(q*tqDgn|PQ(biP&RUsfMuBJxt2fdJ8%j7FzV*7FRZ
zi?+G)gbC3D*i^3>0$?PDbFi>J5k3L`Ic@t>%7BdJ&#c>iY&3e>Zsc=fLH*@G_EVn$
z#H=^%X@9{$vQE8stkEC_R=yrNevJ;<Z|j{Z&#~~Cy{uiNjtP_1;jW&<>bQ)?9j-TU
zFe#SMxMx2HySF$!Df%}Ioo|iD)rp=Uzf+cEx5yI1>HfV1b&m=AmHzFT9~(yGS|xeG
z5>DlDPi}rvTy|ci!x{&u%FOsHD%{5jjFcZ4$V9ombHN&iEa)}GZgaP(K+g;`=c;ls
zH{7>kN+A&)UcX*;wH70pd#^>eXF23W$JWxg6ChaSR5&$Ih@`7(LJc{2uud_!x8Jq^
zteY0x&n@e*Lb*Bm3|}lZyxIQQOxT-*k%-Zo+uJJeZaJx4Hnak_a>qLCEh^D6nx!(c
zI|p2HX=l%Rr$G9tv}9Ch5=5sYsde8Y(b5^UI&4P-KBi>9<B}#L;IjF|64y#R2-w!%
z9GVI_?yb9R{k+g~;7IesRc8@FsROAs3zng0s$L8DK`0?9F2jfvg6x}m`_|gUgDe?Y
z;K577ySk#!3YJ;0El#2D$}fj(x4uxDeg#w_%`*#f3en9oU9?1sia&D~E%*kqz;khp
z@Tf~HWY+)DY_(59`45%*wrUUID<yuWQ0^X1jM^CX-m8Gkszwd9vzd_d_St{-@MVZ*
z=}-=lWJ1xInEK0xPDYHxhEiLVLP$2%T?!eEMeOcdd5_AIaI;ryrptkhBR?#+^GDKf
zy@k3_fK!L0PkYDhm)F3$(9>wBtQ_$bYni_8<+$RPs%KDNfGKzNvSQmpJmKk^y>T}V
z4F05qcfQ1Gl_>Y4TT8GoxRgm+5f_EcdI5VM*XBWAUgJUHhDxIL+?{@Mhl;7l1EuQH
zQP{KY{Nx7?88vqH6X*Aw!@s(}a_4S7L}yl3(Z;JZXlg!t5xTbslVKBC_sgmgVLEGD
zJ6ethUu&v0-&SL}w0Cn)H3c$SAs!Ff3$b+YZu7P!g%EfiO4=asC=b%2u1g0S%TcSM
zKfPja045wq=N1gJu-v?FR{1p<$qUg+cDHEIPv|IfYJQ4*^~-|IoIs>JRbSk=F9jZc
zWP_JF73kb_%iUik6V6jdg{2yEQBYDXu6?r-(bebwiGQLX)v@bqNJ}}M=LP3o*k6R9
z7?Rv+)u9}in+AB*W+XuGPr_46mjrCDOitqw&4j%Ey`DMyFbI8D**(1ZA<@aTm9kB%
z;P*n*SjsgA=S5z&@A*wdu3OYv*@<ezEw;Pq-)Ms6;NOqR`*U!7`6*L9!CJ&d1(;X<
ztbu)#y>hK<70Mp$7dT&53^qyalt_78J``ru-nM=#M4h(Qk%8|IP`z5lsh^RFNoUKf
zvHCQeTpQ)x6`YS{itmkjQ_7%xYG{ZSmjRb3+UQh63c>>yKL~^t!###)wJ)m{T-DZI
z5&H_T{DHkd-WDom_{|PI^v!`#p=RqLiBbq_br*#DmV@MO!WSU9iw3Wo8-7fhQX!RT
z`8AkT1%r)4yNn)YL;n6xxxuJRTwm6bb0Dk=vmqCxo4j*C8y{G+d{Zgn`CM-o-^hcm
zhr*SSoh1l7oT7Sa_&SahSIP=DIpUAI#1r4yEcCpcUZ27&1@9!W<CYAt({i?mhm?`<
z-I{Ox0seTL+~htm>l2TZdV!tYuPafU-fAT!)&!mQvz?KmHTbPNsX6hi90z#BwEcr9
z*sZOtqW`o8j4-jy%uC5|q`K`om`#PW(>C=hEd}s;uFn#sH-Z_W9LWE@5lcgyG8)Xo
zaj^f*xRLJ}{HJxX%g~gRkM2ZH)XQ0KsQqz!ckOX9v?fEm_<v<Wf_DGAm|!*B_<c3w
z%Q9fKTz&uRpCOnYH{PwWmkPn^{aO2OmEmA(*V2-;6x{vXRVx`*g6*kB3huqcD&V`x
za;OsXtPjb9+0{66(QkI=>tYC{<i=iUs|F(|bbkfOm4?0HT1}O-6dag6I5Dy#8M~hR
zSLl1X6rBdum#1Z2;B01ekbfoNe-o2Cs*YaR>#QNhyD}VcVM9UzVR^U~ZD3;=Q;2FO
z@y;*&`FKC!|MAqjREQQiIs|E_B8Otw;pI>Z!QUC3ne-McN%xk^ze5A>-zY!QLD2xL
z{*yp4$gaibSB9-8(ksy)pRn$%HVt~?RgXr61MztxDqVOD6+ahSKR-U022JJ)i|0S$
zFihPyXVnk{FXLR5qSN`fSX{cTMkfP{XSrm5<xrvfA<yNgN;&3v)7eLwvoOCb;du~`
zA3mmwGrZF?;p)dGS**0pg(}~bzjIbpkYfTAep1Oel^O8-k!LO@ejnSL>=A_Es7)?u
z*QhYO7a>u#qXIfljTh|99bnp!9k%QfVMCiw`GL<;aQ0{W;F_JKP*zB<8jCpv$;V$?
z@}H)|^qna?%jyoMN4ak2wuQj(MLO>u5*5V#A%_)pf-p89qu5lb0pni}{CdfW$o$pT
zv#gs?`RhMX)ej5tw5aL$x^Oafn5;OtYZnENe{*NIQ);l&zbm?-qZDortCTB`mcx9H
zKwxHXE~fd8^Y2-e2)%&tl~O@zaB&!`^5m<4V(TFT|2D$+C;1%JqmyAss#ABCHT1`^
z#X-(BkvpglqqNGtrD5_`kMJ@=?Y;`s^Umj~n7m|0jwa&wkd#6G^Tkp~J}Gv)97{n?
zg_`^8z<QL$+RYDG=OcypPi=N41<&aL@B2?u(585L`J=cZC^Vp@#g~G05^`h;Zvxyy
z4){-z9)ff#?L<XOMg)?(%l>Pya>tRolam&;G0@iVb`#9M2WiiY?{tkk7%X{rtV^*R
zYK=`_&g`m%-Q$g?UB1-d#|Nvmyh6Ds`aLVUU9}3de2FVb<8`pU^XV$7r3QW($qr_{
zmDsYw<H7;Y3dCC8P}p*)2!33orhcs=62|fag?@d#kEOgdIzI)A(4AC26|tuh?qhE@
zj)a%O?)3ZznNJy5c-C(;T1<vb?B>i%eaRsAtd2Bl_5g9o_J{Kcc%S@PN=xY>f@Lo6
z=A}2onx1e<(7FKMZktP9AeAEcV&mu=O&S7782zYQi@`%mvrm(V5h#5!K3htJNrdw4
z74-UioL$G&suf*{-2IM&S!(%^IqY%k`gkBj{byt2kCzZ#Ua>#yQ!)02ZvIgpl!rTc
zZ-e=EL?YcYm`S!OLE-fmX&00VAk+75$>i@0%pXwPCjGqv!PFPgao4h;_Mg|eg?+V%
zxudgAj&!XNl$L}z$)aiqHs+}2$<;&EI`zEjOg&;f8f>d&YY?s5@cO{k93=na@qO{K
z62a@d<yh=O1a`)MXb2}<(>buhhCdG9UPu?nn@Pd3<_*`;`~t{jZC%=BOU5m}u;J_0
zH3%BdxwTrY90jj7eRy3$LB6D@KZ&_F7fZS>j%wLlg^{|2%n|)$T$Jm%;jQx+QtS3B
zBd7?TdcXE42|t0chq!E=d^t)w=7u9(YoVO6Rq68QQgnGxGJd{F#nn3+2|HqGsNVK#
zKxmGNnpVNV7jso`IJZn`&?5uxnOQHp#dA?V*xUA_!5y~N+ew*~f*H7F*lp^!tr8@^
z?N&`;xoGkbmRa^W3v<*ogM%p-z-^S)c2n;@Vu%pz;L9K;^Wo5!vsCcq3vPNHRDs_z
zRhE}dG(z=Lj!3951@5__8BHq*!uh-la5gJ|vlx$tkY*KJR%NBI2~>Z~ZqsaS&w@of
zNzT|NwhTNyVn>U<R=`A~ymZflQdrtJI&}Uhf`dNu?(5AJC?4ozF6?*++K#PPQw<-$
z)oby(SDqbS4DcK}!xx7m5|4k1>SjW*An&&8)fn_l{<$fzu?{E3-!%W>RAD->qP>u4
zM4E*I|3z2PpzaiMc)s=#mXcOwi?2|oK*=%s)mQTx?3$~*9{;im*-f)YH(jd5{_pK?
zpYF<mQ$f|+Yi*UNi43cq*hWR~atFmN-4z&_s7(1<m4$B2uS&Y>s8F}xcK^+2E_%OV
z%>|z%P>f%WsvnBRmGhP#>V--$_K#$4aK09WK_CAb^*4f)ypa6GpLm(VLP-4z`x-3$
z;2Y|&qZtvaqn%BgD~R^{iz8c~hi|$D`aU%k*sW1K{U?V4ztxYqq|&p|FLkEZk&yuB
zHMcJw_(8#oLI1VK41@8!L7(|$G7Qt?(px2aZbRd^X8TdDa_BjV?Gt%N1#P3blfs(^
zV({*|%+;7g!KK4J3#%zv5T&i#q_I$3*>jhhUzXo9@z1)8(%Ppkcui}q`T&()&7A2+
gQrgXnC)?A^TfxQkG09~G4=P*u>B}l#C2!&R9|OUw&Hw-a

diff --git a/test/regression/basic_weight_windows/input.py b/test/regression/basic_weight_windows/input.py
index c730df1bc..38d3e9805 100644
--- a/test/regression/basic_weight_windows/input.py
+++ b/test/regression/basic_weight_windows/input.py
@@ -55,9 +55,9 @@
 # Weight windows
 ww_array = np.ones((1, 20, 20, 1, 3))
 # Actual bounds are set to arbitrary numbers
-ww_array[..., 0] = 0.4
-ww_array[..., 1] = 0.7
-ww_array[..., 2] = 0.9
+ww_array[..., 0] = 0.55  # Forces roulette on split particles from 1.0
+ww_array[..., 1] = 0.7  # arbitrary in the middle
+ww_array[..., 2] = 0.9  # forces splitting on all particles born with w=1.0
 mcdc.simulation.weight_windows(ww_array, mesh=mesh)
 
 mcdc.run()

From 344a6fd15102dd38c517d767c55a9c9111cbf2aa Mon Sep 17 00:00:00 2001
From: Nathan Glaser <glaserna@oregonstate.edu>
Date: Wed, 15 Apr 2026 20:49:48 -0700
Subject: [PATCH 37/53] fix typo on weight_roulette method call in simulation

---
 mcdc/transport/simulation.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/mcdc/transport/simulation.py b/mcdc/transport/simulation.py
index 0ddb2bd7e..15a71acf0 100644
--- a/mcdc/transport/simulation.py
+++ b/mcdc/transport/simulation.py
@@ -355,7 +355,7 @@ def step_particle(particle_container, program, data):
 
         # Weight roulette
         else:
-            technique.weight_windows(particle_container, simulation, data)
+            technique.weight_roulette(particle_container, simulation)
 
 
 @njit

From c4598b11e3376a5e884c0c21433a8c358d6475e0 Mon Sep 17 00:00:00 2001
From: Nathan Glaser <glaserna@oregonstate.edu>
Date: Thu, 16 Apr 2026 17:33:29 -0700
Subject: [PATCH 38/53] create forced collision object

---
 mcdc/mcdc_get/__init__.py          |  2 ++
 mcdc/mcdc_get/forced_collisions.py | 32 ++++++++++++++++++++++++++++++
 mcdc/mcdc_set/__init__.py          |  2 ++
 mcdc/mcdc_set/forced_collisions.py | 32 ++++++++++++++++++++++++++++++
 mcdc/numba_types.py                |  7 +++++++
 mcdc/object_/simulation.py         |  3 +++
 mcdc/object_/technique.py          | 31 +++++++++++++++++++++++++++++
 7 files changed, 109 insertions(+)
 create mode 100644 mcdc/mcdc_get/forced_collisions.py
 create mode 100644 mcdc/mcdc_set/forced_collisions.py

diff --git a/mcdc/mcdc_get/__init__.py b/mcdc/mcdc_get/__init__.py
index d3b533388..a82ba3d95 100644
--- a/mcdc/mcdc_get/__init__.py
+++ b/mcdc/mcdc_get/__init__.py
@@ -86,6 +86,8 @@
 
 import mcdc.mcdc_get.settings as settings
 
+import mcdc.mcdc_get.forced_collisions as forced_collisions
+
 import mcdc.mcdc_get.implicit_capture as implicit_capture
 
 import mcdc.mcdc_get.population_control as population_control
diff --git a/mcdc/mcdc_get/forced_collisions.py b/mcdc/mcdc_get/forced_collisions.py
new file mode 100644
index 000000000..f482ecc55
--- /dev/null
+++ b/mcdc/mcdc_get/forced_collisions.py
@@ -0,0 +1,32 @@
+# The following is automatically generated by code_factory.py
+
+from numba import njit
+
+
+@njit
+def cell_IDs(index, forced_collisions, data):
+    offset = forced_collisions["cell_IDs_offset"]
+    return data[offset + index]
+
+
+@njit
+def cell_IDs_all(forced_collisions, data):
+    start = forced_collisions["cell_IDs_offset"]
+    size = forced_collisions["cell_IDs_length"]
+    end = start + size
+    return data[start:end]
+
+
+@njit
+def cell_IDs_last(forced_collisions, data):
+    start = forced_collisions["cell_IDs_offset"]
+    size = forced_collisions["cell_IDs_length"]
+    end = start + size
+    return data[end - 1]
+
+
+@njit
+def cell_IDs_chunk(start, length, forced_collisions, data):
+    start += forced_collisions["cell_IDs_offset"]
+    end = start + length
+    return data[start:end]
diff --git a/mcdc/mcdc_set/__init__.py b/mcdc/mcdc_set/__init__.py
index 38ce0dddd..737ee2a57 100644
--- a/mcdc/mcdc_set/__init__.py
+++ b/mcdc/mcdc_set/__init__.py
@@ -86,6 +86,8 @@
 
 import mcdc.mcdc_set.settings as settings
 
+import mcdc.mcdc_set.forced_collisions as forced_collisions
+
 import mcdc.mcdc_set.implicit_capture as implicit_capture
 
 import mcdc.mcdc_set.population_control as population_control
diff --git a/mcdc/mcdc_set/forced_collisions.py b/mcdc/mcdc_set/forced_collisions.py
new file mode 100644
index 000000000..04657889a
--- /dev/null
+++ b/mcdc/mcdc_set/forced_collisions.py
@@ -0,0 +1,32 @@
+# The following is automatically generated by code_factory.py
+
+from numba import njit
+
+
+@njit
+def cell_IDs(index, forced_collisions, data, value):
+    offset = forced_collisions["cell_IDs_offset"]
+    data[offset + index] = value
+
+
+@njit
+def cell_IDs_all(forced_collisions, data, value):
+    start = forced_collisions["cell_IDs_offset"]
+    size = forced_collisions["cell_IDs_length"]
+    end = start + size
+    data[start:end] = value
+
+
+@njit
+def cell_IDs_last(forced_collisions, data, value):
+    start = forced_collisions["cell_IDs_offset"]
+    size = forced_collisions["cell_IDs_length"]
+    end = start + size
+    data[end - 1] = value
+
+
+@njit
+def cell_IDs_chunk(start, length, forced_collisions, data, value):
+    start += forced_collisions["cell_IDs_offset"]
+    end = start + length
+    data[start:end] = value
diff --git a/mcdc/numba_types.py b/mcdc/numba_types.py
index 8ca7ba4f4..e9271af56 100644
--- a/mcdc/numba_types.py
+++ b/mcdc/numba_types.py
@@ -561,6 +561,12 @@
     ('gpu_storage', int64),
 ])
 
+forced_collisions = into_dtype([
+    ('active', bool),
+    ('cell_IDs_offset', int64),
+    ('cell_IDs_length', int64),
+])
+
 implicit_capture = into_dtype([
     ('active', bool),
 ])
@@ -813,6 +819,7 @@ def set_simulation(N: dict):
         ('N_tracklength_tally', int64),
         ('settings', settings),
         ('implicit_capture', implicit_capture),
+        ('forced_collisions', forced_collisions),
         ('weighted_emission', weighted_emission),
         ('weight_roulette', weight_roulette),
         ('population_control', population_control),
diff --git a/mcdc/object_/simulation.py b/mcdc/object_/simulation.py
index fa5b5f9bd..3c13a481e 100644
--- a/mcdc/object_/simulation.py
+++ b/mcdc/object_/simulation.py
@@ -3,6 +3,7 @@
 
 from mcdc.object_.technique import (
     ImplicitCapture,
+    ForcedCollisions,
     PopulationControl,
     WeightRoulette,
     WeightedEmission,
@@ -80,6 +81,7 @@ class Simulation(ObjectSingleton):
 
     # Techniques
     implicit_capture: ImplicitCapture
+    forced_collisions: ForcedCollisions
     weighted_emission: WeightedEmission
     weight_roulette: WeightRoulette
     population_control: PopulationControl
@@ -165,6 +167,7 @@ def __init__(self):
 
         # Techniques
         self.implicit_capture = ImplicitCapture()
+        self.forced_collisions = ForcedCollisions()
         self.weighted_emission = WeightedEmission()
         self.weight_roulette = WeightRoulette()
         self.population_control = PopulationControl()
diff --git a/mcdc/object_/technique.py b/mcdc/object_/technique.py
index 1c3bacc1b..d0fa2d13a 100644
--- a/mcdc/object_/technique.py
+++ b/mcdc/object_/technique.py
@@ -1,4 +1,10 @@
+from typing import TYPE_CHECKING, Annotated
+from numpy.typing import NDArray
+import numpy as np
 from mcdc.object_.base import ObjectSingleton
+if TYPE_CHECKING:
+    from mcdc.object_.cell import Cell
+from mcdc.constant import FILL_MATERIAL
 from mcdc.print_ import print_error
 
 # ======================================================================================
@@ -18,6 +24,31 @@ def __call__(self, active: bool = True):
         self.active = active
 
 
+# ======================================================================================
+# ForcedCollisions
+# ======================================================================================
+
+
+class ForcedCollisions(ObjectSingleton):
+    #Annotations for Numba mode
+    label: str = "forced_collisions"
+    active: bool
+
+    cell_IDs: list[np.int64]
+
+    def __init__(self):
+        self.active = False
+        self.cell_IDs = []
+    
+    def __call__(self, cells):
+        for cell in cells:
+            if cell.fill_type != FILL_MATERIAL:
+                print_error(
+                    f"Invalid cell fill on cell: \n{cell}\nForced collision technique is only valid on cells with material fill"
+                )
+            self.cell_IDs.append(cell.ID)
+
+
 # ======================================================================================
 # Weighted emission
 # ======================================================================================

From 2b924d3c5a652a96b0ab5f0a3af66374ac77d477 Mon Sep 17 00:00:00 2001
From: Nathan Glaser <glaserna@oregonstate.edu>
Date: Thu, 16 Apr 2026 18:13:47 -0700
Subject: [PATCH 39/53] add notion for forced_collision_distance

---
 mcdc/transport/physics/interface.py | 38 +++++++++++++++++++++++------
 1 file changed, 31 insertions(+), 7 deletions(-)

diff --git a/mcdc/transport/physics/interface.py b/mcdc/transport/physics/interface.py
index ef3cef349..453e22c89 100644
--- a/mcdc/transport/physics/interface.py
+++ b/mcdc/transport/physics/interface.py
@@ -30,6 +30,20 @@ def particle_speed(particle_container, simulation, data):
 # ======================================================================================
 
 
+@njit
+def total_xs(particle_container, simulation, data):
+    particle = particle_container[0]
+    if particle["particle_type"] == PARTICLE_NEUTRON:
+        module = neutron
+        type_total = NEUTRON_REACTION_TOTAL
+    elif particle["particle_type"] == PARTICLE_ELECTRON:
+        module = electron
+        type_total = ELECTRON_REACTION_TOTAL
+    else:
+        return 0.0
+    return module.macro_xs(type_total, particle_container, simulation, data)
+
+
 @njit
 def macro_xs(reaction_type, particle_container, simulation, data):
     particle = particle_container[0]
@@ -57,14 +71,8 @@ def neutron_production_xs(reaction_type, particle_container, simulation, data):
 
 @njit
 def collision_distance(particle_container, simulation, data):
-    particle = particle_container[0]
-
     # Get total cross-section
-    SigmaT = 0.0
-    if particle["particle_type"] == PARTICLE_NEUTRON:
-        SigmaT = macro_xs(NEUTRON_REACTION_TOTAL, particle_container, simulation, data)
-    elif particle["particle_type"] == PARTICLE_ELECTRON:
-        SigmaT = macro_xs(ELECTRON_REACTION_TOTAL, particle_container, simulation, data)
+    SigmaT = total_xs(particle_container, simulation, data)
 
     # Vacuum material?
     if SigmaT == 0.0:
@@ -76,6 +84,22 @@ def collision_distance(particle_container, simulation, data):
     return distance
 
 
+@njit
+def forced_collision_distance(particle_container, surface_distance, simulation, data):
+    # Get total cross-section
+    SigmaT = total_xs(particle_container, simulation, data)
+
+    # Vacuum material?
+    if SigmaT == 0.0:
+        return INF
+
+    # Sample collision distance
+    xi = rng.lcg(particle_container)
+    
+    distance = - math.log(1 - xi*(1-math.exp(-surface_distance * SigmaT))) / SigmaT
+    return distance
+
+
 @njit
 def collision(particle_container, collision_data_container, program, data):
     particle = particle_container[0]

From 93672174736be763f6f72ad45b76e2eef8b4b53e Mon Sep 17 00:00:00 2001
From: Nathan Glaser <glaserna@oregonstate.edu>
Date: Thu, 16 Apr 2026 18:31:35 -0700
Subject: [PATCH 40/53] first pass at implementation of forced-collisions

---
 mcdc/transport/simulation.py |  5 +++-
 mcdc/transport/technique.py  | 47 ++++++++++++++++++++++++++++++++++++
 2 files changed, 51 insertions(+), 1 deletion(-)

diff --git a/mcdc/transport/simulation.py b/mcdc/transport/simulation.py
index 6f2e3f06f..6eefc25e0 100644
--- a/mcdc/transport/simulation.py
+++ b/mcdc/transport/simulation.py
@@ -401,7 +401,10 @@ def move_to_event(particle_container, simulation, data):
     )
 
     # Distance to next collision
-    d_collision = physics.collision_distance(particle_container, simulation, data)
+    if technique.in_forced_collision_cell(particle_container, simulation, data):
+        d_collision = technique.forced_collisions(particle_container, d_boundary, simulation, data)
+    else:
+        d_collision = physics.collision_distance(particle_container, simulation, data)
 
     # ==================================================================================
     # Determine event(s)
diff --git a/mcdc/transport/technique.py b/mcdc/transport/technique.py
index fe629af33..eb4d3bc5d 100644
--- a/mcdc/transport/technique.py
+++ b/mcdc/transport/technique.py
@@ -9,7 +9,54 @@
 import mcdc.transport.particle as particle_module
 import mcdc.transport.particle_bank as particle_bank_module
 import mcdc.transport.rng as rng
+from mcdc.transport.physics import interface as physics
 import mcdc.transport.util as util
+import mcdc.mcdc_get as mcdc_get
+
+# ======================================================================================
+# Forced Collisions
+# ======================================================================================
+
+
+@njit
+def forced_collisions(particle_container, surface_distance, program, data):
+    simulation = util.access_simulation(program)
+    fc_object = simulation["forced_collisions"]
+    SigmaT = physics.total_xs(particle_container, simulation, data)
+
+    # create collided and transmitted particles
+    collided_container = particle_container
+    collided = collided_container[0]
+
+    transmitted_container = util.local_array(1, type_.particle_data)
+    transmitted = transmitted_container[0]
+    particle_module.copy_as_child(transmitted_container, collided_container)
+    
+    # update transmitted particle
+    weight_multiplier = math.exp(-surface_distance * SigmaT)
+    transmitted["w"] *= weight_multiplier
+    particle_module.move(transmitted_container, surface_distance, simulation, data)
+    particle_bank_module.bank_active_particle(transmitted_container, program)
+
+    # update collided particle
+    collided["w"] *= (1 - weight_multiplier)
+    # return distance to forced collision, let simulation handle the rest (tallies)
+    return physics.forced_collision_distance(collided_container, surface_distance, simulation, data) 
+
+
+@njit
+def in_forced_collision_cell(particle_container, simulation, data):
+    fc_object = simulation["forced_collisions"]
+    # not active, dont need to query cells
+    if not fc_object["active"]:
+        return False
+    # active, need to check if in active cell
+    cell_ids = mcdc_get.forced_collisions.cell_IDs_all(fc_object, data)
+    if particle_container[0]["cell_ID"] not in cell_ids:
+        return True
+    # not in active cell
+    return False
+
 
 # ======================================================================================
 # Weight Roulette

From f2e7817a955757ef4b90341cf39791250fb9bf15 Mon Sep 17 00:00:00 2001
From: Nathan Glaser <glaserna@oregonstate.edu>
Date: Thu, 16 Apr 2026 18:58:09 -0700
Subject: [PATCH 41/53] rely on roulette to kill particles instead of tracking
 only surface crossings

---
 mcdc/mcdc_get/forced_collisions.py | 58 ++++++++++++++++++++++++++++++
 mcdc/mcdc_set/forced_collisions.py | 58 ++++++++++++++++++++++++++++++
 mcdc/numba_types.py                |  4 +++
 mcdc/object_/technique.py          | 16 ++++++++-
 mcdc/transport/simulation.py       |  2 ++
 mcdc/transport/technique.py        | 14 ++++++++
 6 files changed, 151 insertions(+), 1 deletion(-)

diff --git a/mcdc/mcdc_get/forced_collisions.py b/mcdc/mcdc_get/forced_collisions.py
index f482ecc55..24e9b1bbd 100644
--- a/mcdc/mcdc_get/forced_collisions.py
+++ b/mcdc/mcdc_get/forced_collisions.py
@@ -30,3 +30,61 @@ def cell_IDs_chunk(start, length, forced_collisions, data):
     start += forced_collisions["cell_IDs_offset"]
     end = start + length
     return data[start:end]
+
+
+@njit
+def threshold_weights(index, forced_collisions, data):
+    offset = forced_collisions["threshold_weights_offset"]
+    return data[offset + index]
+
+
+@njit
+def threshold_weights_all(forced_collisions, data):
+    start = forced_collisions["threshold_weights_offset"]
+    size = forced_collisions["threshold_weights_length"]
+    end = start + size
+    return data[start:end]
+
+
+@njit
+def threshold_weights_last(forced_collisions, data):
+    start = forced_collisions["threshold_weights_offset"]
+    size = forced_collisions["threshold_weights_length"]
+    end = start + size
+    return data[end - 1]
+
+
+@njit
+def threshold_weights_chunk(start, length, forced_collisions, data):
+    start += forced_collisions["threshold_weights_offset"]
+    end = start + length
+    return data[start:end]
+
+
+@njit
+def target_weights(index, forced_collisions, data):
+    offset = forced_collisions["target_weights_offset"]
+    return data[offset + index]
+
+
+@njit
+def target_weights_all(forced_collisions, data):
+    start = forced_collisions["target_weights_offset"]
+    size = forced_collisions["target_weights_length"]
+    end = start + size
+    return data[start:end]
+
+
+@njit
+def target_weights_last(forced_collisions, data):
+    start = forced_collisions["target_weights_offset"]
+    size = forced_collisions["target_weights_length"]
+    end = start + size
+    return data[end - 1]
+
+
+@njit
+def target_weights_chunk(start, length, forced_collisions, data):
+    start += forced_collisions["target_weights_offset"]
+    end = start + length
+    return data[start:end]
diff --git a/mcdc/mcdc_set/forced_collisions.py b/mcdc/mcdc_set/forced_collisions.py
index 04657889a..eda92e70e 100644
--- a/mcdc/mcdc_set/forced_collisions.py
+++ b/mcdc/mcdc_set/forced_collisions.py
@@ -30,3 +30,61 @@ def cell_IDs_chunk(start, length, forced_collisions, data, value):
     start += forced_collisions["cell_IDs_offset"]
     end = start + length
     data[start:end] = value
+
+
+@njit
+def threshold_weights(index, forced_collisions, data, value):
+    offset = forced_collisions["threshold_weights_offset"]
+    data[offset + index] = value
+
+
+@njit
+def threshold_weights_all(forced_collisions, data, value):
+    start = forced_collisions["threshold_weights_offset"]
+    size = forced_collisions["threshold_weights_length"]
+    end = start + size
+    data[start:end] = value
+
+
+@njit
+def threshold_weights_last(forced_collisions, data, value):
+    start = forced_collisions["threshold_weights_offset"]
+    size = forced_collisions["threshold_weights_length"]
+    end = start + size
+    data[end - 1] = value
+
+
+@njit
+def threshold_weights_chunk(start, length, forced_collisions, data, value):
+    start += forced_collisions["threshold_weights_offset"]
+    end = start + length
+    data[start:end] = value
+
+
+@njit
+def target_weights(index, forced_collisions, data, value):
+    offset = forced_collisions["target_weights_offset"]
+    data[offset + index] = value
+
+
+@njit
+def target_weights_all(forced_collisions, data, value):
+    start = forced_collisions["target_weights_offset"]
+    size = forced_collisions["target_weights_length"]
+    end = start + size
+    data[start:end] = value
+
+
+@njit
+def target_weights_last(forced_collisions, data, value):
+    start = forced_collisions["target_weights_offset"]
+    size = forced_collisions["target_weights_length"]
+    end = start + size
+    data[end - 1] = value
+
+
+@njit
+def target_weights_chunk(start, length, forced_collisions, data, value):
+    start += forced_collisions["target_weights_offset"]
+    end = start + length
+    data[start:end] = value
diff --git a/mcdc/numba_types.py b/mcdc/numba_types.py
index b2ec721af..46c5c4d01 100644
--- a/mcdc/numba_types.py
+++ b/mcdc/numba_types.py
@@ -565,6 +565,10 @@
     ('active', bool),
     ('cell_IDs_offset', int64),
     ('cell_IDs_length', int64),
+    ('threshold_weights_offset', int64),
+    ('threshold_weights_length', int64),
+    ('target_weights_offset', int64),
+    ('target_weights_length', int64),
 ])
 
 implicit_capture = into_dtype([
diff --git a/mcdc/object_/technique.py b/mcdc/object_/technique.py
index 3b85b5f5a..75a751faf 100644
--- a/mcdc/object_/technique.py
+++ b/mcdc/object_/technique.py
@@ -36,18 +36,32 @@ class ForcedCollisions(ObjectSingleton):
     active: bool
 
     cell_IDs: list[np.int64]
+    threshold_weights: list[float]
+    target_weights: list[float]
 
     def __init__(self):
         self.active = False
         self.cell_IDs = []
+        self.threshold_weights = []
+        self.target_weights = []
     
-    def __call__(self, cells):
+    def __call__(self, cells, threshold_weights=None, target_weights=None):
+        if threshold_weights is None:
+            threshold_weights = [0.5] * len(cells)
+        if target_weights is None:
+            target_weights = [1.0] * len(cells)
+        if len(cells) != len(threshold_weights) or len(cells) != len(target_weights):
+            print_error(
+                f"Expected cells, threshold_weights, and target_weights to be the same size, but got {len(cells)}, {len(threshold_weights)}, and {len(target_weights)} instead"
+            )
         for cell in cells:
             if cell.fill_type != FILL_MATERIAL:
                 print_error(
                     f"Invalid cell fill on cell: \n{cell}\nForced collision technique is only valid on cells with material fill"
                 )
             self.cell_IDs.append(cell.ID)
+        self.threshold_weights = threshold_weights
+        self.target_weights = target_weights
 
 
 # ======================================================================================
diff --git a/mcdc/transport/simulation.py b/mcdc/transport/simulation.py
index c1a2ea7a4..807229505 100644
--- a/mcdc/transport/simulation.py
+++ b/mcdc/transport/simulation.py
@@ -353,6 +353,8 @@ def step_particle(particle_container, program, data):
         if simulation["weight_windows"]["active"]:
             technique.weight_windows(particle_container, program, data)
 
+        elif simulation["forced_collisions"]["active"]:
+            technique.forced_collision_roulette(particle_container, program, data)
         # Weight roulette
         else:
             technique.weight_roulette(particle_container, simulation)
diff --git a/mcdc/transport/technique.py b/mcdc/transport/technique.py
index 4b76ef303..da27f18a6 100644
--- a/mcdc/transport/technique.py
+++ b/mcdc/transport/technique.py
@@ -46,6 +46,20 @@ def forced_collisions(particle_container, surface_distance, program, data):
     return physics.forced_collision_distance(collided_container, surface_distance, simulation, data) 
 
 
+@njit
+def forced_collision_roulette(particle_container, program, data):
+    simulation = util.access_simulation(program)
+    if not in_forced_collision_cell(particle_container, simulation, data):
+        return
+    particle = particle_container[0]
+    fc_object = simulation["forced_collisions"]
+    cell_ids = mcdc_get.forced_collisions.cell_IDs_all(fc_object, data)
+    index = cell_ids.index(particle["cell_ID"])
+    threshold = mcdc_get.forced_collisions.threshold_weights(index ,fc_object, data)
+    target = mcdc_get.forced_collisions.target_weights(index, fc_object, data)
+    roulette_from_weight_bounds(particle_container, threshold, target)
+
+
 @njit
 def in_forced_collision_cell(particle_container, simulation, data):
     fc_object = simulation["forced_collisions"]

From 05cc1b197c0938e7d8cd16c2a67d7768145aed3d Mon Sep 17 00:00:00 2001
From: Nathan Glaser <glaserna@oregonstate.edu>
Date: Thu, 16 Apr 2026 23:29:10 -0700
Subject: [PATCH 42/53] error throwing on the object for forced collisions

---
 mcdc/object_/technique.py                     |   5 +-
 .../transport/technique/forced_collisions.py  | 104 ++++++++++++++++++
 2 files changed, 108 insertions(+), 1 deletion(-)
 create mode 100644 test/unit/transport/technique/forced_collisions.py

diff --git a/mcdc/object_/technique.py b/mcdc/object_/technique.py
index 75a751faf..f8fb8c219 100644
--- a/mcdc/object_/technique.py
+++ b/mcdc/object_/technique.py
@@ -54,12 +54,15 @@ def __call__(self, cells, threshold_weights=None, target_weights=None):
             print_error(
                 f"Expected cells, threshold_weights, and target_weights to be the same size, but got {len(cells)}, {len(threshold_weights)}, and {len(target_weights)} instead"
             )
+
+        cell_ids = []
         for cell in cells:
             if cell.fill_type != FILL_MATERIAL:
                 print_error(
                     f"Invalid cell fill on cell: \n{cell}\nForced collision technique is only valid on cells with material fill"
                 )
-            self.cell_IDs.append(cell.ID)
+            cell_ids.append(cell.ID)
+        self.cell_IDs = cell_ids
         self.threshold_weights = threshold_weights
         self.target_weights = target_weights
 
diff --git a/test/unit/transport/technique/forced_collisions.py b/test/unit/transport/technique/forced_collisions.py
new file mode 100644
index 000000000..7e9b81654
--- /dev/null
+++ b/test/unit/transport/technique/forced_collisions.py
@@ -0,0 +1,104 @@
+import numpy as np
+import os
+import sys
+import pytest
+
+os.environ["MCDC_LIB"] = "../../../regression/mcdc-regression_test_data/"
+
+# =============================================================================
+# State helpers
+# =============================================================================
+
+
+def _clear_mcdc_modules():
+    for name in list(sys.modules):
+        if name == "mcdc" or name.startswith("mcdc."):
+            del sys.modules[name]
+
+@pytest.fixture
+def fresh_mcdc():
+    _clear_mcdc_modules()
+    import mcdc
+    yield mcdc 
+    _clear_mcdc_modules()
+
+
+# =============================================================================
+# Model base fixture
+# =============================================================================
+
+
+@pytest.fixture
+def pin_cell_model(fresh_mcdc):
+  mcdc = fresh_mcdc
+  # Material
+  fuel = mcdc.Material(
+      nuclide_composition={
+          "U235": 1.0 
+      }
+  )
+  moderator = mcdc.Material(
+      nuclide_composition={
+          "H1": 1.0
+      }
+  )
+
+  # Geometry
+  cylinder = mcdc.Surface.CylinderZ(radius=0.5)
+  pitch = 1.5
+  x0 = mcdc.Surface.PlaneX(x=-pitch / 2, boundary_condition="reflective")
+  x1 = mcdc.Surface.PlaneX(x=pitch / 2, boundary_condition="reflective")
+  y0 = mcdc.Surface.PlaneY(y=-pitch / 2, boundary_condition="reflective")
+  y1 = mcdc.Surface.PlaneY(y=pitch / 2, boundary_condition="reflective")
+  #
+  fuel_cell = mcdc.Cell(-cylinder, fill=fuel)
+  mod_cell = mcdc.Cell(+x0 & -x1 & +y0 & -y1 & +cylinder, fill=moderator)
+
+  yield mcdc, fuel_cell, mod_cell
+
+
+# =============================================================================
+# Error throwing in object creation
+# =============================================================================
+
+
+@pytest.mark.parametrize(
+    "cells_builder, thresholds, targets, expected_msg",
+    [
+        (
+            lambda fuel_cell, mcdc: [fuel_cell],
+            [0.5, 0.5],
+            [1.0],
+            "Expected cells, threshold_weights, and target_weights to be the same size",
+        ),
+        (
+            lambda fuel_cell, mcdc: [fuel_cell],
+            [0.5],
+            [1.0, 1.0],
+            "Expected cells, threshold_weights, and target_weights to be the same size",
+        ),
+        (
+            lambda fuel_cell, mcdc: [mcdc.Cell(fill=mcdc.Universe(cells=[fuel_cell]))],
+            None,
+            None,
+            "Invalid cell fill on cell",
+        ),
+    ],
+)
+def test_forced_collisions_error_throw(
+    pin_cell_model, capsys, cells_builder, thresholds, targets, expected_msg
+):
+    mcdc, fuel_cell, mod_cell = pin_cell_model
+
+    cells = cells_builder(fuel_cell, mcdc)
+
+    with pytest.raises(SystemExit):
+        mcdc.simulation.forced_collisions(
+            cells,
+            threshold_weights=thresholds,
+            target_weights=targets,
+        )
+
+    captured = capsys.readouterr()
+    assert expected_msg in captured.out
+

From 36ae825ba68100da35ba08eea282a1062ab9a75a Mon Sep 17 00:00:00 2001
From: Nathan Glaser <glaserna@oregonstate.edu>
Date: Thu, 16 Apr 2026 23:56:08 -0700
Subject: [PATCH 43/53] attempt at forced collision cell test

---
 mcdc/object_/technique.py                     |  4 +-
 .../transport/technique/forced_collisions.py  | 59 +++++++++++--------
 2 files changed, 34 insertions(+), 29 deletions(-)

diff --git a/mcdc/object_/technique.py b/mcdc/object_/technique.py
index f8fb8c219..a8744778a 100644
--- a/mcdc/object_/technique.py
+++ b/mcdc/object_/technique.py
@@ -55,14 +55,12 @@ def __call__(self, cells, threshold_weights=None, target_weights=None):
                 f"Expected cells, threshold_weights, and target_weights to be the same size, but got {len(cells)}, {len(threshold_weights)}, and {len(target_weights)} instead"
             )
 
-        cell_ids = []
         for cell in cells:
             if cell.fill_type != FILL_MATERIAL:
                 print_error(
                     f"Invalid cell fill on cell: \n{cell}\nForced collision technique is only valid on cells with material fill"
                 )
-            cell_ids.append(cell.ID)
-        self.cell_IDs = cell_ids
+            self.cell_IDs.append(cell.ID)
         self.threshold_weights = threshold_weights
         self.target_weights = target_weights
 
diff --git a/test/unit/transport/technique/forced_collisions.py b/test/unit/transport/technique/forced_collisions.py
index 7e9b81654..82f99153c 100644
--- a/test/unit/transport/technique/forced_collisions.py
+++ b/test/unit/transport/technique/forced_collisions.py
@@ -1,27 +1,15 @@
 import numpy as np
 import os
-import sys
 import pytest
+import mcdc
+from mcdc.main import preparation
+from mcdc.transport import technique
+import mcdc.numba_types as types_
+from mcdc.transport.util import access_simulation
+from mcdc.transport import geometry
 
 os.environ["MCDC_LIB"] = "../../../regression/mcdc-regression_test_data/"
 
-# =============================================================================
-# State helpers
-# =============================================================================
-
-
-def _clear_mcdc_modules():
-    for name in list(sys.modules):
-        if name == "mcdc" or name.startswith("mcdc."):
-            del sys.modules[name]
-
-@pytest.fixture
-def fresh_mcdc():
-    _clear_mcdc_modules()
-    import mcdc
-    yield mcdc 
-    _clear_mcdc_modules()
-
 
 # =============================================================================
 # Model base fixture
@@ -29,8 +17,7 @@ def fresh_mcdc():
 
 
 @pytest.fixture
-def pin_cell_model(fresh_mcdc):
-  mcdc = fresh_mcdc
+def pin_cell_model():
   # Material
   fuel = mcdc.Material(
       nuclide_composition={
@@ -54,7 +41,7 @@ def pin_cell_model(fresh_mcdc):
   fuel_cell = mcdc.Cell(-cylinder, fill=fuel)
   mod_cell = mcdc.Cell(+x0 & -x1 & +y0 & -y1 & +cylinder, fill=moderator)
 
-  yield mcdc, fuel_cell, mod_cell
+  yield fuel_cell, mod_cell
 
 
 # =============================================================================
@@ -66,19 +53,19 @@ def pin_cell_model(fresh_mcdc):
     "cells_builder, thresholds, targets, expected_msg",
     [
         (
-            lambda fuel_cell, mcdc: [fuel_cell],
+            lambda fuel_cell: [fuel_cell],
             [0.5, 0.5],
             [1.0],
             "Expected cells, threshold_weights, and target_weights to be the same size",
         ),
         (
-            lambda fuel_cell, mcdc: [fuel_cell],
+            lambda fuel_cell: [fuel_cell],
             [0.5],
             [1.0, 1.0],
             "Expected cells, threshold_weights, and target_weights to be the same size",
         ),
         (
-            lambda fuel_cell, mcdc: [mcdc.Cell(fill=mcdc.Universe(cells=[fuel_cell]))],
+            lambda fuel_cell: [mcdc.Cell(fill=mcdc.Universe(cells=[fuel_cell]))],
             None,
             None,
             "Invalid cell fill on cell",
@@ -88,9 +75,9 @@ def pin_cell_model(fresh_mcdc):
 def test_forced_collisions_error_throw(
     pin_cell_model, capsys, cells_builder, thresholds, targets, expected_msg
 ):
-    mcdc, fuel_cell, mod_cell = pin_cell_model
+    fuel_cell, mod_cell = pin_cell_model
 
-    cells = cells_builder(fuel_cell, mcdc)
+    cells = cells_builder(fuel_cell)
 
     with pytest.raises(SystemExit):
         mcdc.simulation.forced_collisions(
@@ -102,3 +89,23 @@ def test_forced_collisions_error_throw(
     captured = capsys.readouterr()
     assert expected_msg in captured.out
 
+
+# =============================================================================
+# Method tests
+# =============================================================================
+
+
+def test_in_forced_collision_cell(pin_cell_model):
+    # instantiate
+    fuel_cell, _ = pin_cell_model
+    mcdc.simulation.forced_collisions([fuel_cell])
+    # preparation
+    program, data = preparation()
+    simulation = access_simulation(program)
+    # make particle
+    particle_container = np.zeros(1, types_.particle)
+    particle = particle_container[0]
+    particle["cell_ID"] = -1
+    # inspect geometry
+    geometry.inspect_geometry(particle_container, simulation, data)
+    assert technique.in_forced_collision_cell(particle_container, simulation, data)
\ No newline at end of file

From 4517cf960a25d2b357685b26eb3deee871013c18 Mon Sep 17 00:00:00 2001
From: Nathan Glaser <glaserna@oregonstate.edu>
Date: Fri, 17 Apr 2026 19:26:43 -0700
Subject: [PATCH 44/53] set active to be true when called...

---
 mcdc/object_/technique.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/mcdc/object_/technique.py b/mcdc/object_/technique.py
index a8744778a..6aaf1e800 100644
--- a/mcdc/object_/technique.py
+++ b/mcdc/object_/technique.py
@@ -46,6 +46,7 @@ def __init__(self):
         self.target_weights = []
     
     def __call__(self, cells, threshold_weights=None, target_weights=None):
+        self.active = True
         if threshold_weights is None:
             threshold_weights = [0.5] * len(cells)
         if target_weights is None:

From 83dd21939018fa8673c708c220ea90227d35adb9 Mon Sep 17 00:00:00 2001
From: Nathan Glaser <glaserna@oregonstate.edu>
Date: Fri, 17 Apr 2026 19:27:15 -0700
Subject: [PATCH 45/53] add copy particle runtime state to particle module

---
 mcdc/transport/particle.py | 11 +++++++++++
 1 file changed, 11 insertions(+)

diff --git a/mcdc/transport/particle.py b/mcdc/transport/particle.py
index e524e1c21..5d6954c01 100644
--- a/mcdc/transport/particle.py
+++ b/mcdc/transport/particle.py
@@ -49,3 +49,14 @@ def copy_as_child(child_particle_container, parent_particle_container):
 
     # Evolve parent seed
     rng.lcg(parent_particle_container)
+
+
+@njit
+def copy_run_state(target_particle_container, source_particle_container):
+    target_particle = target_particle_container[0]
+    source_particle = source_particle_container[0]
+    target_particle["alive"] = source_particle["alive"]
+    target_particle["material_ID"] = source_particle["material_ID"]
+    target_particle["cell_ID"] = source_particle["cell_ID"]
+    target_particle["surface_ID"] = source_particle["surface_ID"]
+    target_particle["event"] = source_particle["event"]
\ No newline at end of file

From d92ce79012f0175df4a374a7e85a77aeade680f4 Mon Sep 17 00:00:00 2001
From: Nathan Glaser <glaserna@oregonstate.edu>
Date: Fri, 17 Apr 2026 19:28:10 -0700
Subject: [PATCH 46/53] update forced collision to generate collided and
 transmitted properly

---
 mcdc/transport/technique.py | 86 +++++++++++++++++++++++++++++--------
 1 file changed, 69 insertions(+), 17 deletions(-)

diff --git a/mcdc/transport/technique.py b/mcdc/transport/technique.py
index da27f18a6..ad7e2edce 100644
--- a/mcdc/transport/technique.py
+++ b/mcdc/transport/technique.py
@@ -8,12 +8,14 @@
 import mcdc.mcdc_get.weight_windows as ww_get
 import mcdc.numba_types as type_
 from mcdc.transport.mesh import get_indices as get_mesh_indices
+import mcdc.transport.geometry as geometry_module
 import mcdc.transport.particle as particle_module
 import mcdc.transport.particle_bank as particle_bank_module
 import mcdc.transport.rng as rng
 from mcdc.transport.physics import interface as physics
 import mcdc.transport.util as util
 import mcdc.mcdc_get as mcdc_get
+from mcdc.print_ import print_error
 
 # ======================================================================================
 # Forced Collisions
@@ -23,53 +25,103 @@
 @njit
 def forced_collisions(particle_container, surface_distance, program, data):
     simulation = util.access_simulation(program)
-    fc_object = simulation["forced_collisions"]
+
+    # find weight multiplier
     SigmaT = physics.total_xs(particle_container, simulation, data)
+    weight_multiplier = math.exp(-surface_distance * SigmaT)
 
-    # create collided and transmitted particles
+    # alias input particle as collided particle
     collided_container = particle_container
+ 
+    # transmitted particle
+    bank_transmitted_particle(
+        collided_container, weight_multiplier, surface_distance, program, data
+    )
+
+    # update collided particle
     collided = collided_container[0]
+    collided["w"] *= (1 - weight_multiplier)
 
-    transmitted_container = util.local_array(1, type_.particle_data)
-    transmitted = transmitted_container[0]
+    # return distance to forced collision, let simulation handle the rest (tallies)
+    collision_distance =  physics.forced_collision_distance(
+        collided_container, surface_distance, simulation, data
+    ) 
+    return collision_distance
+
+
+@njit
+def bank_transmitted_particle(collided_container, weight_multiplier, surface_distance, program, data):
+    simulation = util.access_simulation(program)
+
+    # create child copy of collided particle history
+    transmitted_container = util.local_array(1, type_.particle)
     particle_module.copy_as_child(transmitted_container, collided_container)
-    
-    # update transmitted particle
-    weight_multiplier = math.exp(-surface_distance * SigmaT)
+    particle_module.copy_run_state(transmitted_container, collided_container)
+
+    # assign weight
+    transmitted = transmitted_container[0]
     transmitted["w"] *= weight_multiplier
+
+    # update position and perform surface crossing
     particle_module.move(transmitted_container, surface_distance, simulation, data)
-    particle_bank_module.bank_active_particle(transmitted_container, program)
+    geometry_module.surface_crossing(transmitted_container, simulation, data)
 
-    # update collided particle
-    collided["w"] *= (1 - weight_multiplier)
-    # return distance to forced collision, let simulation handle the rest (tallies)
-    return physics.forced_collision_distance(collided_container, surface_distance, simulation, data) 
+    # particle could leave through BC, so check if alive before banking
+    if transmitted["alive"]:
+        particle_bank_module.bank_active_particle(transmitted_container, program)
 
 
 @njit
 def forced_collision_roulette(particle_container, program, data):
     simulation = util.access_simulation(program)
+    fc_object = simulation["forced_collisions"]
+
+    # check if particle is in a cell with forced collisions
     if not in_forced_collision_cell(particle_container, simulation, data):
         return
-    particle = particle_container[0]
-    fc_object = simulation["forced_collisions"]
-    cell_ids = mcdc_get.forced_collisions.cell_IDs_all(fc_object, data)
-    index = cell_ids.index(particle["cell_ID"])
+    
+    # get index into arrays
+    index = get_forced_collision_cell_index(particle_container, fc_object, data) 
+
+    # get weights
     threshold = mcdc_get.forced_collisions.threshold_weights(index ,fc_object, data)
     target = mcdc_get.forced_collisions.target_weights(index, fc_object, data)
+
+    # roulette
     roulette_from_weight_bounds(particle_container, threshold, target)
 
 
+@njit
+def get_forced_collision_cell_index(particle_container, fc_object, data):
+    particle = particle_container[0]
+    
+    # grab all cell ids
+    cell_ids = mcdc_get.forced_collisions.cell_IDs_all(fc_object, data)
+
+    # find the cell index
+    for index in range(len(cell_ids)):
+        if cell_ids[index] == particle["cell_ID"]:
+            return index
+    
+    # should never hit this, but just to be safe
+    print_error(
+        f"Failed to find {particle['cell_ID']} in list of cell ids for forced collisions"
+    )
+
+
 @njit
 def in_forced_collision_cell(particle_container, simulation, data):
     fc_object = simulation["forced_collisions"]
+
     # not active, dont need to query cells
     if not fc_object["active"]:
         return False
+
     # active, need to check if in active cell
     cell_ids = mcdc_get.forced_collisions.cell_IDs_all(fc_object, data)
-    if particle_container[0]["cell_ID"] not in cell_ids:
+    if particle_container[0]["cell_ID"] in cell_ids:
         return True
+
     # not in active cell
     return False
 

From d4cc65f9dd8b65e4800c674a0270578012860e17 Mon Sep 17 00:00:00 2001
From: Nathan Glaser <glaserna@oregonstate.edu>
Date: Fri, 17 Apr 2026 20:33:48 -0700
Subject: [PATCH 47/53] factor out tracklength score procedure in simulation,
 score transmitted particle

---
 mcdc/transport/simulation.py  | 36 ++-----------------------------
 mcdc/transport/tally/score.py | 40 +++++++++++++++++++++++++++++++++++
 mcdc/transport/technique.py   |  6 ++++++
 3 files changed, 48 insertions(+), 34 deletions(-)

diff --git a/mcdc/transport/simulation.py b/mcdc/transport/simulation.py
index 807229505..d2517cf3d 100644
--- a/mcdc/transport/simulation.py
+++ b/mcdc/transport/simulation.py
@@ -447,40 +447,8 @@ def move_to_event(particle_container, simulation, data):
     # ==================================================================================
     # Move particle
     # ==================================================================================
-
     # Score tracklength tallies
-    if simulation["cycle_active"]:
-        # Cell tallies
-        cell = simulation["cells"][particle["cell_ID"]]
-        for i in range(cell["N_tally"]):
-            tally_base_ID = int(mcdc_get.cell.tally_IDs(i, cell, data))
-            tally_base = simulation["tallies"][tally_base_ID]
-
-            # Skip non-tracklength tallies
-            if tally_base["child_type"] != TALLY_TRACKLENGTH:
-                continue
-
-            tally = simulation["tracklength_tallies"][tally_base["child_ID"]]
-            tally_module.score.tracklength_tally(
-                particle_container, distance, tally, simulation, data
-            )
-
-        # Other tracklength tallies
-        for i in range(simulation["N_tracklength_tally"]):
-            tally = simulation["tracklength_tallies"][i]
-
-            # Skip cell tallies
-            if tally["spatial_filter_type"] == SPATIAL_FILTER_CELL:
-                continue
-
-            tally_module.score.tracklength_tally(
-                particle_container, distance, tally, simulation, data
-            )
-
-    if settings["neutron_eigenvalue_mode"]:
-        tally_module.score.eigenvalue_tally(
-            particle_container, distance, simulation, data
-        )
+    tally_module.score.score_tracklength_tallies(particle_container, distance, simulation, data)
 
     # Move particle
-    particle_module.move(particle_container, distance, simulation, data)
+    particle_module.move(particle_container, distance, simulation, data)
\ No newline at end of file
diff --git a/mcdc/transport/tally/score.py b/mcdc/transport/tally/score.py
index 924c03386..4c1fea6f6 100644
--- a/mcdc/transport/tally/score.py
+++ b/mcdc/transport/tally/score.py
@@ -26,6 +26,8 @@
     SCORE_NET_CURRENT,
     SCORE_ENERGY_DEPOSITION,
     SPATIAL_FILTER_MESH,
+    SPATIAL_FILTER_CELL,
+    TALLY_TRACKLENGTH,
 )
 from mcdc.transport.geometry.surface import get_normal_component
 from mcdc.transport.tally.filter import get_filter_indices
@@ -140,6 +142,44 @@ def collision_tally(
 # ======================================================================================
 
 
+@njit
+def score_tracklength_tallies(particle_container, distance, simulation, data):
+    particle = particle_container[0]
+
+    if simulation["cycle_active"]:
+        # Cell tallies
+        cell = simulation["cells"][particle["cell_ID"]]
+        for i in range(cell["N_tally"]):
+            tally_base_ID = int(mcdc_get.cell.tally_IDs(i, cell, data))
+            tally_base = simulation["tallies"][tally_base_ID]
+
+            # Skip non-tracklength tallies
+            if tally_base["child_type"] != TALLY_TRACKLENGTH:
+                continue
+
+            tally = simulation["tracklength_tallies"][tally_base["child_ID"]]
+            tracklength_tally(
+                particle_container, distance, tally, simulation, data
+            )
+
+        # Other tracklength tallies
+        for i in range(simulation["N_tracklength_tally"]):
+            tally = simulation["tracklength_tallies"][i]
+
+            # Skip cell tallies
+            if tally["spatial_filter_type"] == SPATIAL_FILTER_CELL:
+                continue
+
+            tracklength_tally(
+                particle_container, distance, tally, simulation, data
+            )
+
+    if simulation["settings"]["neutron_eigenvalue_mode"]:
+        eigenvalue_tally(
+            particle_container, distance, simulation, data
+        )
+
+
 @njit
 def tracklength_tally(particle_container, distance, tally, simulation, data):
     particle = particle_container[0]
diff --git a/mcdc/transport/technique.py b/mcdc/transport/technique.py
index ad7e2edce..c2a1517a3 100644
--- a/mcdc/transport/technique.py
+++ b/mcdc/transport/technique.py
@@ -11,6 +11,7 @@
 import mcdc.transport.geometry as geometry_module
 import mcdc.transport.particle as particle_module
 import mcdc.transport.particle_bank as particle_bank_module
+import mcdc.transport.tally as tally_module
 import mcdc.transport.rng as rng
 from mcdc.transport.physics import interface as physics
 import mcdc.transport.util as util
@@ -62,6 +63,11 @@ def bank_transmitted_particle(collided_container, weight_multiplier, surface_dis
     transmitted = transmitted_container[0]
     transmitted["w"] *= weight_multiplier
 
+    # score tracklength tallies
+    tally_module.score.score_tracklength_tallies(
+        transmitted_container, surface_distance, simulation, data
+    )
+
     # update position and perform surface crossing
     particle_module.move(transmitted_container, surface_distance, simulation, data)
     geometry_module.surface_crossing(transmitted_container, simulation, data)

From 319095306870c8f8a48fea373a35b751151029b5 Mon Sep 17 00:00:00 2001
From: Nathan Glaser <glaserna@oregonstate.edu>
Date: Fri, 17 Apr 2026 20:40:16 -0700
Subject: [PATCH 48/53] add regression test for forced collisions

---
 .../slab_beam_forced_collision/answer.h5      | Bin 0 -> 41200 bytes
 .../slab_beam_forced_collision/input.py       |  63 ++++++++++++++++++
 2 files changed, 63 insertions(+)
 create mode 100644 test/regression/slab_beam_forced_collision/answer.h5
 create mode 100644 test/regression/slab_beam_forced_collision/input.py

diff --git a/test/regression/slab_beam_forced_collision/answer.h5 b/test/regression/slab_beam_forced_collision/answer.h5
new file mode 100644
index 0000000000000000000000000000000000000000..6d43b155f7dd17b3ae577e0671b2136e4093014d
GIT binary patch
literal 41200
zcmeHQZ;Ts96(9R@;aXI1!WTFE+w_wIq_`n~3l-=5t!YD2B`I74!RPbs`JDLH&e|LK
z0-`B@iW{K1v<TfNa0`h0fJ9UYDH0zzRR~4!ffFiHg#exSfFdMX!-psmQp3!?_w4i7
zyL)FGXV;!R$!})f+u3=)eQ##=&CcHQsl!JfxPRCEOdJyvY#WQ%N&2~iN23lUlqdc)
zCPGMSNLRv;2r==;SQN{5VEwY(>61r}FfK{;6UEb5d>f=Ds6f$?LLhbIz*9ouDbzZS
z^ak1)X4ysEEa!^`E9;d?-pG}yypkyt@_Lz>mQl$s=(bG?UW^1BXo!u2EE;%-fd!Ji
zk8#Tu=?6n0+p>lC01xEGxxC1S!%#edo-jrgN&d#zvO`P!QGDw7iNhh@OAu~rFF?Bp
zkK34!UWIb%O={cQ<)YXjT*Gn=*ld1pyWE<6*BNV^gmOTYF>MS_2q^KktB7qd5w0-k
zoU}D<*;Cktyk$=jM3fLBVCt0cHFX*%WILgLU={tohBW>dsMm=rz!KvADGstzqU66j
zw8<vU``^j--$DCw#(PS$E5x%PEq-gQWB^ja{e+E)Ltygwp%ce!Y|wEDCv;O5fB2Do
z565|bW}iW}$(sEwx#c+W3z-?cz|%N<v-_~zC>m+QTA0aTp5*B|MXOS>DvVJc#hPUL
z-OZ#*C<GJ&3IT<HLO>y)5Kssx1QY@a0fm4<Kq1g00uYa%aE!+%c0t8F6ks3MF>Y7F
ze0(1y)-ZqllaSCbpwRdX<fk!TJq-D46mWO~@)H#2IRW`91s)))e(w<i6bGhgW0e9p
z4(-9<B)j}x2LUYBkS@K?U5@Oo{(g5kzjHg|m4GO@!5Cj35b4dwfukDeOm~J5s8hn%
zwEuW$8Yryr5IBs5VhRtALu=h69&kLmudk=Rz?J03{28E5>nrP7;IlcTNNAX)eyE+x
zS2^ATC9^W4E#G*iE$f|WObR!p1qo<GNURx9kyu9%Lk;D{dV^Rapu9Z{b!{|B;EIBl
zu7?k{)`J?q=9TcHF&W)^`1F5w-?xf1X%{n_kLN%M|E0X3uWdFTTKMJR-M=`O{C)M-
zq1wB)C!x09ufP27=Pv#9&Rdhe+P?h2pI_fL`Ty2;?$P|7*xPSSes*$x$3NdZi}N(7
zRsAXi6aoqXg}|^Q(EoXOLIcdyJY3)~7D^5|9`Ih=*59j1IpDG^VDs(Ke}8#O17(NR
zUt)G#Bg*IaV!AG7=P+%bU`G!eJ9UJYoj&=rAU4DKJ**#VwhW#PY4KZYC4G=`dj68%
zZ*-EIKLtwEo(4+sw*=~T?E~&grjQ+txywx<JGRGNj^ezD54p?vo$Kwl{j$?(4X^+!
zZ{gxMJKfy-&Rg%-GpPSWOLaRnYW;fN0f)MTi$Al}?<`^SQXh72=L)jdoL=Na{&lM6
zqm<}QX+0@I3jp24d|&F_?`$>he?bGKW&f)C9t+XH{^PNzf#M1efx}qnSxwN^8{3d$
ztWF7EQ*Zm|^?BoY4HTE{#NXo*^}`LkYuY0ow-2Ws)sMN$MNqD`-(5~ax#f?$%lVzV
zS-aZ${;sPU=pK<E@V8xg<I%7A@XLr#(C5Sb!sC@d;6cWGJe|o_@)z}1ekN1N&a>m`
ztZtO8=DtgwU(6PCUY<8Lo^d=~%9xeBVA%)o(lZ(3e0s*3o6}7whuSulr_D?yU!?0I
zUu^GOOwXD61xq)wi;PjZ*y||{y`X0+Cf{^aF*8QFR5bbCB5}-F70c9}>(7;}bf&y$
zWYd+!60b_+%ax*;$w9r4=S|Xc{HmP4XgAN`t(?tUWju?q^IZQtzb5U6GIYysj-Joy
z#>Gs*($fpYS)Eb&f>o&Gb7s*hwU<lHXYxKO*f*6nGN^aSELJ+TX3R8l>9Vdjmv76N
zi~7b}QZ8C%)&&>j&*ckx>tbAgvAH%qQ?!iPj41%I$~syE$GEJ|rtLcg<q}o`K3kr*
zIcMtSGCu>Y0LnL9vl;D4*k^72hsMnWb~`Jda6dOm-P5D)aw(LHe9B#J8s)T4yUY2V
z+c{1|?IzPcOCHv<52Ccsa_i&hPc*0r31s+ht(o#eq7x68Pm=lcuK#Z|&`-@H1eW2j
zyzlvS#D|auw1Q_JTYcXD6Y4ELx6}y$UsLLbt9aX{KH}coC6r4f-R0^impb4sNADR=
z9dwuTJJ-A4>DT=9uNq*105^dNW<R|3`1RjVf7#!H{yf*~_e=hPcnsUP)N35+Rv_(z
z`r#Ys*7Jv&j@>o!2Fk4->B?>v-bA^zqs^kO3EJPW{&-h*U<wD<55KJeC_(Rs`!z1z
zQFw5JXuZ$Qb}ReXh({<8ctAfq&NO^rZlvHYI{n5`F7>371eB*xZt9r3+%(F~f6iTw
z_TN>HyUX1xuGe=2+2^5U+P!L^y6!zifW99#g{>eri{SfJ)VOSse=kscNe}{-={PcY
z{55;Npc}c$d>f~1nR6K)Foi&LoYCk^vtuv2Tt5RBke@E7J&F!D0<vAkP~8=L-;cD<
znRz>AuN%6VTTIXDrD8eH*Fy{xAiE=brSqI9j5}k^uh=a~-(Qs9JEVB~)K{TKy542X
z*S`kk8eVWc-!3<C4)QVaQMe&Gy!<=yB~X^|8zR+}Dw{w+_U|N~Q^S6rs#vffmAM3t
z1AY$cYJyWmL9w<bj!QtA<4qLvM8Lnp&9BVf{2*}1Z{_`;H2$nga`Jvss((~y=kq-v
z?X}*u^*CdAH(*LL2ST_;313jKt9LZ*>hf9lW|3Vkec4@(>@W63cR3B&^|ZSjjicA5
z-Q{S!T+O)4C6K)~<1V*Z?%?cdXCUlqaCkhV@Zd(;iW}`yc(mbhP~p*r$4P}p8y+)(
zz{Bggq?Y4K@k5?%K&kp2b_C?-5}K!JbM2Zh=ovf2P@dH<3Qh>wNu*mPo#Q;r{HR`f
zUux^~XA}LBW|GvhWlH#x(sLTMx+>u>H=6u)ZC>V#cK+_!)ulk%mD7G=IS%JFPE6w)
z*JPigZDpgUz0)|6lH#}&pPnmNRe=)q=gA%+bF$wXc-@;A`1Z2j>ldH_T955erv)1U
z*?wcF&O~!32^4t6U&)>+UgU)ejq5aC$L-lu%P*!dU$woZwL7zp_D$E5?cM9T!^ZC%
z*I}mG|M=_Iah#_+Mo88u;R_0Owc)=tHg{ZgZyun<pYw3vQ+KJB{jHRoB%oYFwub8!
z?$Z9{4$iJ#41`?`4v$w99^5Tj@z+-af`^!qb^bqBc^vbakI}-V+{^m6r%)~HCj#<(
zjK-&BoCndhf)jUo4gFVsLS80Sd`W43O08~4_{)tZe@%Tw=8Sf}bKC?(MPDh7-s73S
z3#45+J-5s6iFu6^Z#)b5k=;Q~{f-I(@;ITP;;Z<cA>lxACqu97&040Z8}Oh*{+A%Z
zO;RWOO+Et^Oru?T*>C)tAg^umx=L)g5RmORhU#5L`yp-bR@{lk^(fjc`Gai7)W5FO
zphYCNhWW-7YUyav&T&0z1JM7t(;s&OvUDd1XS@NFy~k9JJ5|5s-XK7We~sZecR8|C
z?HlfL8nTza?Jnnc?%?d|?}4zZ!Qt^Ag$H+zR@^E45crC}`yHqGm^?mt&Bv<W?XIog
z`wP{+{vaUF$7p<te5bws;KZGFqW{teN_s3)!k1KzJIzb@%Z;YqY4v+DXSDO3<0jYv
z&EI>EJM9XjT{(>tw_gA|iCzWjHBQuD1UcEOM%&s(M|*jkNTKQpv>!sEQsg_nd8TCQ
z{DY?A%cSr^QaR7R;MvMD_%0`7_6JTEG6w&C=z@q*k)Klk*6Gjf&zmGqw0Hkr`vK@6
zTX>CB^96-*Vw@$gHS0feZ&n=HYd>_CBfM|E<SwTn`}X(U<p{+k$v)iA8XR|UcJ)Re
z>}qg$tOo**!CB(nz0%&F|K-Y)XLlsuy!GG{cVa)AJoK+K<1b%|B-vv-xBcPDZzuVe
op<gZT33I9rg@8gpA)pXY2q*-G1%biM>uv^OUN<;A8iBy$e@5sbTmS$7

literal 0
HcmV?d00001

diff --git a/test/regression/slab_beam_forced_collision/input.py b/test/regression/slab_beam_forced_collision/input.py
new file mode 100644
index 000000000..ca346fbbb
--- /dev/null
+++ b/test/regression/slab_beam_forced_collision/input.py
@@ -0,0 +1,63 @@
+import numpy as np
+import os
+import mcdc
+
+os.environ["MCDC_LIB"] = "../mcdc-regression_test_data/"
+
+# ======================================================================================
+# Set model
+# ======================================================================================
+# Finite homogeneous pure-absorbing slab
+
+# Set materials
+# Set materials
+generate_material = lambda atomdensity: mcdc.Material(nuclide_composition={"H1": atomdensity})
+m1 = generate_material(0.0001)
+
+# Set surfaces
+s1 = mcdc.Surface.PlaneX(x=0.0, boundary_condition="vacuum")
+s2 = mcdc.Surface.PlaneX(x=1.0, boundary_condition="vacuum")
+
+# Set cells
+low_xs_cell = mcdc.Cell(region=+s1 & -s2, fill=m1)
+
+# ======================================================================================
+# Set source
+# ======================================================================================
+# Isotropic beam from left-end
+
+mcdc.Source(
+    position=(0.0, 0.0, 0.0),
+    direction=(1.0, 0.0, 0.0)
+)
+
+# ======================================================================================
+# Set tallies, settings, and run MC/DC
+# ======================================================================================
+
+# Tallies
+# energy deposition for actual VR against analog
+mesh = mcdc.MeshUniform()
+mcdc.Tally(
+    mesh=mesh,
+    scores=["energy_deposition"],
+)
+# flux to make sure tracklength is unbiased
+mcdc.Tally(
+  cell=low_xs_cell,
+  scores=["flux"]
+)
+# net-current to make sure surface is unbiased
+mcdc.Tally(
+  surface=s2,
+  scores=["net-current"]
+)
+
+# Settings
+mcdc.settings.N_particle = 5000
+mcdc.settings.N_batch = 2
+
+mcdc.simulation.forced_collisions(cells=[low_xs_cell])
+
+# Run
+mcdc.run()

From 388d5642276f700466a033ab25f87ba40ca3b492 Mon Sep 17 00:00:00 2001
From: Nathan Glaser <glaserna@oregonstate.edu>
Date: Fri, 17 Apr 2026 21:29:06 -0700
Subject: [PATCH 49/53] remove unit test for in cell check for forced
 collisions

---
 .../transport/technique/forced_collisions.py  | 58 ++++++-------------
 1 file changed, 18 insertions(+), 40 deletions(-)

diff --git a/test/unit/transport/technique/forced_collisions.py b/test/unit/transport/technique/forced_collisions.py
index 82f99153c..498b3c83f 100644
--- a/test/unit/transport/technique/forced_collisions.py
+++ b/test/unit/transport/technique/forced_collisions.py
@@ -1,12 +1,14 @@
 import numpy as np
 import os
 import pytest
+from numba import njit
 import mcdc
 from mcdc.main import preparation
 from mcdc.transport import technique
 import mcdc.numba_types as types_
 from mcdc.transport.util import access_simulation
 from mcdc.transport import geometry
+from mcdc.mcdc_get import forced_collisions as get_fc
 
 os.environ["MCDC_LIB"] = "../../../regression/mcdc-regression_test_data/"
 
@@ -18,30 +20,22 @@
 
 @pytest.fixture
 def pin_cell_model():
-  # Material
-  fuel = mcdc.Material(
-      nuclide_composition={
-          "U235": 1.0 
-      }
-  )
-  moderator = mcdc.Material(
-      nuclide_composition={
-          "H1": 1.0
-      }
-  )
-
-  # Geometry
-  cylinder = mcdc.Surface.CylinderZ(radius=0.5)
-  pitch = 1.5
-  x0 = mcdc.Surface.PlaneX(x=-pitch / 2, boundary_condition="reflective")
-  x1 = mcdc.Surface.PlaneX(x=pitch / 2, boundary_condition="reflective")
-  y0 = mcdc.Surface.PlaneY(y=-pitch / 2, boundary_condition="reflective")
-  y1 = mcdc.Surface.PlaneY(y=pitch / 2, boundary_condition="reflective")
-  #
-  fuel_cell = mcdc.Cell(-cylinder, fill=fuel)
-  mod_cell = mcdc.Cell(+x0 & -x1 & +y0 & -y1 & +cylinder, fill=moderator)
-
-  yield fuel_cell, mod_cell
+    # Material
+    fuel = mcdc.Material(nuclide_composition={"U235": 1.0})
+    moderator = mcdc.Material(nuclide_composition={"H1": 1.0})
+
+    # Geometry
+    cylinder = mcdc.Surface.CylinderZ(radius=0.5)
+    pitch = 1.5
+    x0 = mcdc.Surface.PlaneX(x=-pitch / 2, boundary_condition="reflective")
+    x1 = mcdc.Surface.PlaneX(x=pitch / 2, boundary_condition="reflective")
+    y0 = mcdc.Surface.PlaneY(y=-pitch / 2, boundary_condition="reflective")
+    y1 = mcdc.Surface.PlaneY(y=pitch / 2, boundary_condition="reflective")
+    #
+    fuel_cell = mcdc.Cell(-cylinder, fill=fuel)
+    mod_cell = mcdc.Cell(+x0 & -x1 & +y0 & -y1 & +cylinder, fill=moderator)
+
+    yield fuel_cell, mod_cell
 
 
 # =============================================================================
@@ -93,19 +87,3 @@ def test_forced_collisions_error_throw(
 # =============================================================================
 # Method tests
 # =============================================================================
-
-
-def test_in_forced_collision_cell(pin_cell_model):
-    # instantiate
-    fuel_cell, _ = pin_cell_model
-    mcdc.simulation.forced_collisions([fuel_cell])
-    # preparation
-    program, data = preparation()
-    simulation = access_simulation(program)
-    # make particle
-    particle_container = np.zeros(1, types_.particle)
-    particle = particle_container[0]
-    particle["cell_ID"] = -1
-    # inspect geometry
-    geometry.inspect_geometry(particle_container, simulation, data)
-    assert technique.in_forced_collision_cell(particle_container, simulation, data)
\ No newline at end of file

From f64c85d7feba08b2820071ab062760efb0cbdd93 Mon Sep 17 00:00:00 2001
From: Nathan Glaser <glaserna@oregonstate.edu>
Date: Fri, 17 Apr 2026 21:59:00 -0700
Subject: [PATCH 50/53] add docstrings to all added methods

---
 mcdc/transport/particle.py          | 10 +++
 mcdc/transport/physics/interface.py | 36 +++++++++++
 mcdc/transport/tally/score.py       | 14 +++++
 mcdc/transport/technique.py         | 94 ++++++++++++++++++++++++++---
 4 files changed, 147 insertions(+), 7 deletions(-)

diff --git a/mcdc/transport/particle.py b/mcdc/transport/particle.py
index 5d6954c01..771289487 100644
--- a/mcdc/transport/particle.py
+++ b/mcdc/transport/particle.py
@@ -53,6 +53,16 @@ def copy_as_child(child_particle_container, parent_particle_container):
 
 @njit
 def copy_run_state(target_particle_container, source_particle_container):
+    """
+    Helper for copying runtime particle state into new particle.
+
+    Parameters
+    ----------
+    target_particle_container : ndarray
+        Container holding the particle to copy data into.
+    source_particle_container : ndarray
+        Container holding the particle to copy data from.
+    """
     target_particle = target_particle_container[0]
     source_particle = source_particle_container[0]
     target_particle["alive"] = source_particle["alive"]
diff --git a/mcdc/transport/physics/interface.py b/mcdc/transport/physics/interface.py
index 453e22c89..1e892254d 100644
--- a/mcdc/transport/physics/interface.py
+++ b/mcdc/transport/physics/interface.py
@@ -32,6 +32,23 @@ def particle_speed(particle_container, simulation, data):
 
 @njit
 def total_xs(particle_container, simulation, data):
+    """
+    Convenience helper for getting specifically the total cross section. 
+
+    Parameters
+    ----------
+    particle_container : ndarray
+        Container holding the particle.
+    simulation : object
+        Simulation object.
+    data : object
+        Simulation data for array access.
+
+    Returns
+    -------
+    float
+        Total macroscopic cross section.
+    """
     particle = particle_container[0]
     if particle["particle_type"] == PARTICLE_NEUTRON:
         module = neutron
@@ -86,6 +103,25 @@ def collision_distance(particle_container, simulation, data):
 
 @njit
 def forced_collision_distance(particle_container, surface_distance, simulation, data):
+    """
+    Method for finding the distance for a forced collision particle to travel.
+
+    Parameters
+    ----------
+    particle_container : ndarray
+        Container holding the particle.
+    surface_distance:
+        The distance to the next surface along the particles direction.
+    simulation : object
+        Simulation object.
+    data : object
+        Simulation data for array access.
+    
+    Returns
+    -------
+    distance : float
+        Distance for particle to travel.
+    """
     # Get total cross-section
     SigmaT = total_xs(particle_container, simulation, data)
 
diff --git a/mcdc/transport/tally/score.py b/mcdc/transport/tally/score.py
index 4c1fea6f6..0b4cae0e7 100644
--- a/mcdc/transport/tally/score.py
+++ b/mcdc/transport/tally/score.py
@@ -144,6 +144,20 @@ def collision_tally(
 
 @njit
 def score_tracklength_tallies(particle_container, distance, simulation, data):
+    """
+    Helper for scoring traveled distance on all track length tallies in the simulation.
+
+    Parameters
+    ----------
+    particle_container : ndarray
+        Container holding the particle.
+    distance:
+        The distance traveled to score.
+    simulation : object
+        Simulation object.
+    data : object
+        Simulation data for array access.
+    """
     particle = particle_container[0]
 
     if simulation["cycle_active"]:
diff --git a/mcdc/transport/technique.py b/mcdc/transport/technique.py
index c2a1517a3..44798e2c8 100644
--- a/mcdc/transport/technique.py
+++ b/mcdc/transport/technique.py
@@ -25,20 +25,38 @@
 
 @njit
 def forced_collisions(particle_container, surface_distance, program, data):
+    """
+    Applies the method of forced collisions, splitting the source particle into a collided and transmitted particle. This method returns the distance for the collided particle to travel, letting the `simulation.move_to_event` handle the actual transport for the collided particle.
+
+    Parameters
+    ----------
+    particle_container : ndarray
+        Container holding the original particle to copy over all data from.
+    surface_distance : float
+        The distance to the surface the transmitted particle will be moved to.
+    program : object
+        Program object containing simulation state with forced collision data.
+    data : object
+        Simulation data for array access.
+
+    Returns
+    -------
+    collision_distance : float
+        Distance for the collided component to travel.
+    """
     simulation = util.access_simulation(program)
 
     # find weight multiplier
     SigmaT = physics.total_xs(particle_container, simulation, data)
     weight_multiplier = math.exp(-surface_distance * SigmaT)
-
-    # alias input particle as collided particle
-    collided_container = particle_container
  
     # transmitted particle
     bank_transmitted_particle(
-        collided_container, weight_multiplier, surface_distance, program, data
+        particle_container, weight_multiplier, surface_distance, program, data
     )
 
+    # alias input particle as collided particle
+    collided_container = particle_container
     # update collided particle
     collided = collided_container[0]
     collided["w"] *= (1 - weight_multiplier)
@@ -51,13 +69,29 @@ def forced_collisions(particle_container, surface_distance, program, data):
 
 
 @njit
-def bank_transmitted_particle(collided_container, weight_multiplier, surface_distance, program, data):
+def bank_transmitted_particle(particle_container, weight_multiplier, surface_distance, program, data):
+    """
+    Helper for creating and banking the transmitted component. If the transmitted particle leaves the simulation through a boundary surface, the particle is not banked. Additionally, the particle is scored over all tracklength tallies via the helper in `tally.score`.
+
+    Parameters
+    ----------
+    particle_container : ndarray
+        Container holding the original particle to copy over all data from.
+    weight_multiplier : float
+        The multiplier to adjust the particle weight by.
+    surface_distance : float
+        The distance to the surface the transmitted particle will be moved to.
+    program : object
+        Program object containing simulation state with forced collision data.
+    data : object
+        Simulation data for array access.
+    """
     simulation = util.access_simulation(program)
 
     # create child copy of collided particle history
     transmitted_container = util.local_array(1, type_.particle)
-    particle_module.copy_as_child(transmitted_container, collided_container)
-    particle_module.copy_run_state(transmitted_container, collided_container)
+    particle_module.copy_as_child(transmitted_container, particle_container)
+    particle_module.copy_run_state(transmitted_container, particle_container)
 
     # assign weight
     transmitted = transmitted_container[0]
@@ -79,6 +113,18 @@ def bank_transmitted_particle(collided_container, weight_multiplier, surface_dis
 
 @njit
 def forced_collision_roulette(particle_container, program, data):
+    """
+    Roulette procedure for forced collision. Particle is only rouletted if in a cell marked for forced collisions
+
+    Parameters
+    ----------
+    particle_container : ndarray
+        Container holding the particle.
+    program : object
+        Program object containing simulation state with forced collision data.
+    data : object
+        Simulation data for array access.
+    """
     simulation = util.access_simulation(program)
     fc_object = simulation["forced_collisions"]
 
@@ -99,6 +145,23 @@ def forced_collision_roulette(particle_container, program, data):
 
 @njit
 def get_forced_collision_cell_index(particle_container, fc_object, data):
+    """
+    Helper for getting the getter index for weight roulette parameters
+
+    Parameters
+    ----------
+    particle_container : ndarray
+        Container holding the particle.
+    fc_object : object
+        Forced collision object for use in mcdc_get methods.
+    data : object
+        Simulation data for array access.
+
+    Returns
+    -------
+    index : int
+        The flattened index for getting weight roulette parameters.
+    """
     particle = particle_container[0]
     
     # grab all cell ids
@@ -117,6 +180,23 @@ def get_forced_collision_cell_index(particle_container, fc_object, data):
 
 @njit
 def in_forced_collision_cell(particle_container, simulation, data):
+    """
+    Check if particle is in a cell marked for forced collision
+
+    Parameters
+    ----------
+    particle_container : ndarray
+        Container holding the particle.
+    program : object
+        Program object containing simulation state with forced collision data.
+    data : object
+        Simulation data for array access.
+
+    Returns
+    -------
+    bool
+        True if particle in cell marked for forced collision.
+    """
     fc_object = simulation["forced_collisions"]
 
     # not active, dont need to query cells

From 41ff45afdc11513ad3b914b376c811bbd689eae8 Mon Sep 17 00:00:00 2001
From: Nathan Glaser <glaserna@oregonstate.edu>
Date: Fri, 17 Apr 2026 22:01:27 -0700
Subject: [PATCH 51/53] back in black

---
 mcdc/object_/technique.py                     |  5 ++--
 mcdc/transport/particle.py                    |  2 +-
 mcdc/transport/physics/interface.py           |  8 +++---
 mcdc/transport/simulation.py                  | 10 ++++---
 mcdc/transport/tally/score.py                 | 12 +++------
 mcdc/transport/technique.py                   | 27 +++++++++++--------
 .../slab_beam_forced_collision/input.py       | 19 +++++--------
 7 files changed, 40 insertions(+), 43 deletions(-)

diff --git a/mcdc/object_/technique.py b/mcdc/object_/technique.py
index 6aaf1e800..af719f6a2 100644
--- a/mcdc/object_/technique.py
+++ b/mcdc/object_/technique.py
@@ -3,6 +3,7 @@
 import numpy as np
 from mcdc.constant import FILL_MATERIAL, INF
 from mcdc.object_.base import ObjectSingleton
+
 if TYPE_CHECKING:
     from mcdc.object_.cell import Cell
 from mcdc.object_.mesh import MeshBase, MeshUniform
@@ -31,7 +32,7 @@ def __call__(self, active: bool = True):
 
 
 class ForcedCollisions(ObjectSingleton):
-    #Annotations for Numba mode
+    # Annotations for Numba mode
     label: str = "forced_collisions"
     active: bool
 
@@ -44,7 +45,7 @@ def __init__(self):
         self.cell_IDs = []
         self.threshold_weights = []
         self.target_weights = []
-    
+
     def __call__(self, cells, threshold_weights=None, target_weights=None):
         self.active = True
         if threshold_weights is None:
diff --git a/mcdc/transport/particle.py b/mcdc/transport/particle.py
index 771289487..edeb0bd67 100644
--- a/mcdc/transport/particle.py
+++ b/mcdc/transport/particle.py
@@ -69,4 +69,4 @@ def copy_run_state(target_particle_container, source_particle_container):
     target_particle["material_ID"] = source_particle["material_ID"]
     target_particle["cell_ID"] = source_particle["cell_ID"]
     target_particle["surface_ID"] = source_particle["surface_ID"]
-    target_particle["event"] = source_particle["event"]
\ No newline at end of file
+    target_particle["event"] = source_particle["event"]
diff --git a/mcdc/transport/physics/interface.py b/mcdc/transport/physics/interface.py
index 1e892254d..8fbf878ca 100644
--- a/mcdc/transport/physics/interface.py
+++ b/mcdc/transport/physics/interface.py
@@ -33,7 +33,7 @@ def particle_speed(particle_container, simulation, data):
 @njit
 def total_xs(particle_container, simulation, data):
     """
-    Convenience helper for getting specifically the total cross section. 
+    Convenience helper for getting specifically the total cross section.
 
     Parameters
     ----------
@@ -116,7 +116,7 @@ def forced_collision_distance(particle_container, surface_distance, simulation,
         Simulation object.
     data : object
         Simulation data for array access.
-    
+
     Returns
     -------
     distance : float
@@ -131,8 +131,8 @@ def forced_collision_distance(particle_container, surface_distance, simulation,
 
     # Sample collision distance
     xi = rng.lcg(particle_container)
-    
-    distance = - math.log(1 - xi*(1-math.exp(-surface_distance * SigmaT))) / SigmaT
+
+    distance = -math.log(1 - xi * (1 - math.exp(-surface_distance * SigmaT))) / SigmaT
     return distance
 
 
diff --git a/mcdc/transport/simulation.py b/mcdc/transport/simulation.py
index d2517cf3d..cd72f7f44 100644
--- a/mcdc/transport/simulation.py
+++ b/mcdc/transport/simulation.py
@@ -411,7 +411,9 @@ def move_to_event(particle_container, simulation, data):
 
     # Distance to next collision
     if technique.in_forced_collision_cell(particle_container, simulation, data):
-        d_collision = technique.forced_collisions(particle_container, d_boundary, simulation, data)
+        d_collision = technique.forced_collisions(
+            particle_container, d_boundary, simulation, data
+        )
     else:
         d_collision = physics.collision_distance(particle_container, simulation, data)
 
@@ -448,7 +450,9 @@ def move_to_event(particle_container, simulation, data):
     # Move particle
     # ==================================================================================
     # Score tracklength tallies
-    tally_module.score.score_tracklength_tallies(particle_container, distance, simulation, data)
+    tally_module.score.score_tracklength_tallies(
+        particle_container, distance, simulation, data
+    )
 
     # Move particle
-    particle_module.move(particle_container, distance, simulation, data)
\ No newline at end of file
+    particle_module.move(particle_container, distance, simulation, data)
diff --git a/mcdc/transport/tally/score.py b/mcdc/transport/tally/score.py
index 0b4cae0e7..6c3633d51 100644
--- a/mcdc/transport/tally/score.py
+++ b/mcdc/transport/tally/score.py
@@ -172,9 +172,7 @@ def score_tracklength_tallies(particle_container, distance, simulation, data):
                 continue
 
             tally = simulation["tracklength_tallies"][tally_base["child_ID"]]
-            tracklength_tally(
-                particle_container, distance, tally, simulation, data
-            )
+            tracklength_tally(particle_container, distance, tally, simulation, data)
 
         # Other tracklength tallies
         for i in range(simulation["N_tracklength_tally"]):
@@ -184,14 +182,10 @@ def score_tracklength_tallies(particle_container, distance, simulation, data):
             if tally["spatial_filter_type"] == SPATIAL_FILTER_CELL:
                 continue
 
-            tracklength_tally(
-                particle_container, distance, tally, simulation, data
-            )
+            tracklength_tally(particle_container, distance, tally, simulation, data)
 
     if simulation["settings"]["neutron_eigenvalue_mode"]:
-        eigenvalue_tally(
-            particle_container, distance, simulation, data
-        )
+        eigenvalue_tally(particle_container, distance, simulation, data)
 
 
 @njit
diff --git a/mcdc/transport/technique.py b/mcdc/transport/technique.py
index 44798e2c8..4d81ccd0b 100644
--- a/mcdc/transport/technique.py
+++ b/mcdc/transport/technique.py
@@ -49,7 +49,7 @@ def forced_collisions(particle_container, surface_distance, program, data):
     # find weight multiplier
     SigmaT = physics.total_xs(particle_container, simulation, data)
     weight_multiplier = math.exp(-surface_distance * SigmaT)
- 
+
     # transmitted particle
     bank_transmitted_particle(
         particle_container, weight_multiplier, surface_distance, program, data
@@ -59,19 +59,24 @@ def forced_collisions(particle_container, surface_distance, program, data):
     collided_container = particle_container
     # update collided particle
     collided = collided_container[0]
-    collided["w"] *= (1 - weight_multiplier)
+    collided["w"] *= 1 - weight_multiplier
 
     # return distance to forced collision, let simulation handle the rest (tallies)
-    collision_distance =  physics.forced_collision_distance(
+    collision_distance = physics.forced_collision_distance(
         collided_container, surface_distance, simulation, data
-    ) 
+    )
     return collision_distance
 
 
 @njit
-def bank_transmitted_particle(particle_container, weight_multiplier, surface_distance, program, data):
+def bank_transmitted_particle(
+    particle_container, weight_multiplier, surface_distance, program, data
+):
     """
-    Helper for creating and banking the transmitted component. If the transmitted particle leaves the simulation through a boundary surface, the particle is not banked. Additionally, the particle is scored over all tracklength tallies via the helper in `tally.score`.
+    Helper for creating and banking the transmitted component. If the
+    transmitted particle leaves the simulation through a boundary surface, the
+    particle is not banked. Additionally, the particle is scored over all
+    tracklength tallies via the helper in `tally.score`.
 
     Parameters
     ----------
@@ -131,12 +136,12 @@ def forced_collision_roulette(particle_container, program, data):
     # check if particle is in a cell with forced collisions
     if not in_forced_collision_cell(particle_container, simulation, data):
         return
-    
+
     # get index into arrays
-    index = get_forced_collision_cell_index(particle_container, fc_object, data) 
+    index = get_forced_collision_cell_index(particle_container, fc_object, data)
 
     # get weights
-    threshold = mcdc_get.forced_collisions.threshold_weights(index ,fc_object, data)
+    threshold = mcdc_get.forced_collisions.threshold_weights(index, fc_object, data)
     target = mcdc_get.forced_collisions.target_weights(index, fc_object, data)
 
     # roulette
@@ -163,7 +168,7 @@ def get_forced_collision_cell_index(particle_container, fc_object, data):
         The flattened index for getting weight roulette parameters.
     """
     particle = particle_container[0]
-    
+
     # grab all cell ids
     cell_ids = mcdc_get.forced_collisions.cell_IDs_all(fc_object, data)
 
@@ -171,7 +176,7 @@ def get_forced_collision_cell_index(particle_container, fc_object, data):
     for index in range(len(cell_ids)):
         if cell_ids[index] == particle["cell_ID"]:
             return index
-    
+
     # should never hit this, but just to be safe
     print_error(
         f"Failed to find {particle['cell_ID']} in list of cell ids for forced collisions"
diff --git a/test/regression/slab_beam_forced_collision/input.py b/test/regression/slab_beam_forced_collision/input.py
index ca346fbbb..7b391604e 100644
--- a/test/regression/slab_beam_forced_collision/input.py
+++ b/test/regression/slab_beam_forced_collision/input.py
@@ -11,7 +11,9 @@
 
 # Set materials
 # Set materials
-generate_material = lambda atomdensity: mcdc.Material(nuclide_composition={"H1": atomdensity})
+generate_material = lambda atomdensity: mcdc.Material(
+    nuclide_composition={"H1": atomdensity}
+)
 m1 = generate_material(0.0001)
 
 # Set surfaces
@@ -26,10 +28,7 @@
 # ======================================================================================
 # Isotropic beam from left-end
 
-mcdc.Source(
-    position=(0.0, 0.0, 0.0),
-    direction=(1.0, 0.0, 0.0)
-)
+mcdc.Source(position=(0.0, 0.0, 0.0), direction=(1.0, 0.0, 0.0))
 
 # ======================================================================================
 # Set tallies, settings, and run MC/DC
@@ -43,15 +42,9 @@
     scores=["energy_deposition"],
 )
 # flux to make sure tracklength is unbiased
-mcdc.Tally(
-  cell=low_xs_cell,
-  scores=["flux"]
-)
+mcdc.Tally(cell=low_xs_cell, scores=["flux"])
 # net-current to make sure surface is unbiased
-mcdc.Tally(
-  surface=s2,
-  scores=["net-current"]
-)
+mcdc.Tally(surface=s2, scores=["net-current"])
 
 # Settings
 mcdc.settings.N_particle = 5000

From 4c4ccdc27ad6403b18a8898dc0c29830fb892458 Mon Sep 17 00:00:00 2001
From: Nathan Glaser <glaserna@oregonstate.edu>
Date: Fri, 17 Apr 2026 22:04:10 -0700
Subject: [PATCH 52/53] add python api line for forced collisions

---
 docs/source/pythonapi/index.rst | 1 +
 1 file changed, 1 insertion(+)

diff --git a/docs/source/pythonapi/index.rst b/docs/source/pythonapi/index.rst
index a0f9a6445..86a687a62 100644
--- a/docs/source/pythonapi/index.rst
+++ b/docs/source/pythonapi/index.rst
@@ -87,6 +87,7 @@ Defining techniques
 Techniques are enabled by calling methods on the ``mcdc.simulation`` singleton:
 
 - ``mcdc.simulation.implicit_capture(active=True)``
+- ``mcdc.simulation.forced_collisions(cells=[], weight_thresholds=[], weight_targets=[])``
 - ``mcdc.simulation.weighted_emission(active=True, weight_target=1.0)``
 - ``mcdc.simulation.weight_roulette(weight_threshold=0.0, weight_target=1.0)``
 - ``mcdc.simulation.weight_windows(weight_windows, mesh=None, energy=None)``

From 575f1e205bafdc4bdd0c9fc844920892cc9ab608 Mon Sep 17 00:00:00 2001
From: Nathan Glaser <glaserna@oregonstate.edu>
Date: Sat, 18 Apr 2026 14:46:29 -0700
Subject: [PATCH 53/53] fix various typos and bugs causing tests to fail

---
 mcdc/transport/physics/interface.py | 12 +++++-------
 mcdc/transport/technique.py         |  7 ++++---
 2 files changed, 9 insertions(+), 10 deletions(-)

diff --git a/mcdc/transport/physics/interface.py b/mcdc/transport/physics/interface.py
index 8fbf878ca..f1e2c8a6a 100644
--- a/mcdc/transport/physics/interface.py
+++ b/mcdc/transport/physics/interface.py
@@ -51,14 +51,12 @@ def total_xs(particle_container, simulation, data):
     """
     particle = particle_container[0]
     if particle["particle_type"] == PARTICLE_NEUTRON:
-        module = neutron
-        type_total = NEUTRON_REACTION_TOTAL
+        total = NEUTRON_REACTION_TOTAL
+        return neutron.macro_xs(total, particle_container, simulation, data)
     elif particle["particle_type"] == PARTICLE_ELECTRON:
-        module = electron
-        type_total = ELECTRON_REACTION_TOTAL
-    else:
-        return 0.0
-    return module.macro_xs(type_total, particle_container, simulation, data)
+        total = ELECTRON_REACTION_TOTAL
+        return electron.macro_xs(total, particle_container, simulation, data)
+    return 0.0
 
 
 @njit
diff --git a/mcdc/transport/technique.py b/mcdc/transport/technique.py
index 4d81ccd0b..a82d0c796 100644
--- a/mcdc/transport/technique.py
+++ b/mcdc/transport/technique.py
@@ -2,6 +2,7 @@
 import math
 
 from numba import njit
+import numba
 
 ####
 
@@ -139,6 +140,8 @@ def forced_collision_roulette(particle_container, program, data):
 
     # get index into arrays
     index = get_forced_collision_cell_index(particle_container, fc_object, data)
+    if index < 0:
+        return
 
     # get weights
     threshold = mcdc_get.forced_collisions.threshold_weights(index, fc_object, data)
@@ -178,9 +181,7 @@ def get_forced_collision_cell_index(particle_container, fc_object, data):
             return index
 
     # should never hit this, but just to be safe
-    print_error(
-        f"Failed to find {particle['cell_ID']} in list of cell ids for forced collisions"
-    )
+    return -1
 
 
 @njit