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 01/20] 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 02/20] 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 03/20] 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 04/20] 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 05/20] 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 06/20] 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 07/20] 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 08/20] 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 09/20] 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 10/20] 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 11/20] 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 12/20] 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 13/20] 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 14/20] 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 15/20] 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 16/20] 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

From 26f9bacf13d8506ed02e82e8d7dda8f8c6eacbc8 Mon Sep 17 00:00:00 2001
From: Nathan Glaser <glaserna@oregonstate.edu>
Date: Wed, 29 Apr 2026 18:50:51 -0700
Subject: [PATCH 17/20] move unit test to object for forced collisions

---
 .../technique/forced_collisions.py            | 21 +++++++------------
 1 file changed, 8 insertions(+), 13 deletions(-)
 rename test/unit/{transport => object_}/technique/forced_collisions.py (83%)

diff --git a/test/unit/transport/technique/forced_collisions.py b/test/unit/object_/technique/forced_collisions.py
similarity index 83%
rename from test/unit/transport/technique/forced_collisions.py
rename to test/unit/object_/technique/forced_collisions.py
index 498b3c83f..6276d5943 100644
--- a/test/unit/transport/technique/forced_collisions.py
+++ b/test/unit/object_/technique/forced_collisions.py
@@ -1,17 +1,6 @@
 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/"
-
 
 # =============================================================================
 # Model base fixture
@@ -21,8 +10,14 @@
 @pytest.fixture
 def pin_cell_model():
     # Material
-    fuel = mcdc.Material(nuclide_composition={"U235": 1.0})
-    moderator = mcdc.Material(nuclide_composition={"H1": 1.0})
+    fuel = mcdc.MaterialMG(
+        capture=np.array([1.0 / 3.0]),
+        scatter=np.array([[1.0 / 3.0]]),
+    )
+    moderator = mcdc.MaterialMG(
+        capture=np.array([1.0 / 3.0]),
+        scatter=np.array([[1.0 / 3.0]]),
+    )
 
     # Geometry
     cylinder = mcdc.Surface.CylinderZ(radius=0.5)

From 91d317b221567d6431ebbb64960090ffed39bcc3 Mon Sep 17 00:00:00 2001
From: Nathan Glaser <glaserna@oregonstate.edu>
Date: Wed, 29 Apr 2026 19:08:35 -0700
Subject: [PATCH 18/20] update theory for VR for forced collisions

---
 docs/source/theory/variance_reduction.rst | 53 +++++++++++++++++++++++
 1 file changed, 53 insertions(+)

diff --git a/docs/source/theory/variance_reduction.rst b/docs/source/theory/variance_reduction.rst
index 84a33183c..0c037aa2a 100644
--- a/docs/source/theory/variance_reduction.rst
+++ b/docs/source/theory/variance_reduction.rst
@@ -42,6 +42,59 @@ effective in highly absorbing media.
    to eliminate very-low-weight particles, a memory overhead can build
    up over time.
 
+
+Forced Collisions
+------------------
+
+.. warning::
+
+   Forced collisions are currently valid only for neutral particles.
+
+Forced collisions force a neutral particle to undergo a collision within selected material cells before it reaches the next surface.
+This is done by splitting the incident particle into two components:
+
+- a transmitted component that travels to the next surface without collision, with weight
+
+  .. math::
+
+     w_{\text{trans}} = w_0 e^{-\Sigma_t d}
+
+- a collided component that undergoes a forced collision in the cell, with weight
+
+  .. math::
+
+     w_{\text{coll}} = w_0 \left(1 - e^{-\Sigma_t d}\right)
+
+where :math:`d` is the distance to the next surface and :math:`\Sigma_t` is the total macroscopic cross section.
+The collision distance is sampled from the following distribution:
+
+.. math::
+
+   s = -\frac{1}{\Sigma_t} \ln \left(1 - \xi \left(1 - e^{-\Sigma_t d}\right)\right)
+
+Additionally, collided particles are continually forced to undergo collisions in the cell of interest.
+To prevent tracking particles with extremely low weight, weight roulette is used.
+
+**Usage:**
+
+.. code-block:: python3
+
+  mcdc.simulation.forced_collisions(cells=[cell])
+
+Optional roulette parameters can be supplied for each forced-collision cell:
+
+.. code-block:: python3
+
+  mcdc.simulation.forced_collisions(
+      cells=[cell_1, cell_2],
+      threshold_weights=[0.25, 0.25],
+      target_weights=[1.0, 1.0]
+  )
+
+If ``threshold_weights`` or ``target_weights`` are omitted, default values of ``0.5`` and ``1.0`` are assumed for each cell.
+Forced collisions may only be enabled on cells with material fills.
+
+
 Weight Roulette
 ----------------
 

From df62de4f8e8c8a93b79840820a6c3840d22f1d48 Mon Sep 17 00:00:00 2001
From: Nathan Glaser <glaserna@oregonstate.edu>
Date: Wed, 29 Apr 2026 19:16:15 -0700
Subject: [PATCH 19/20] add check to bypass forced collisions if not neutron

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

diff --git a/mcdc/transport/technique.py b/mcdc/transport/technique.py
index 6fbd99ced..f25a65b72 100644
--- a/mcdc/transport/technique.py
+++ b/mcdc/transport/technique.py
@@ -17,7 +17,7 @@
 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
+from mcdc.constant import PARTICLE_NEUTRON
 
 # ======================================================================================
 # Forced Collisions
@@ -131,6 +131,10 @@ def forced_collision_roulette(particle_container, program, data):
     data : object
         Simulation data for array access.
     """
+    # skip if not a neutron
+    if particle_container[0]["particle_type"] != PARTICLE_NEUTRON:
+        return
+
     simulation = util.access_simulation(program)
     fc_object = simulation["forced_collisions"]
 
@@ -203,6 +207,10 @@ def in_forced_collision_cell(particle_container, simulation, data):
     bool
         True if particle in cell marked for forced collision.
     """
+    # skip if not a neutron
+    if particle_container[0]["particle_type"] != PARTICLE_NEUTRON:
+        return False
+    
     fc_object = simulation["forced_collisions"]
 
     # not active, dont need to query cells

From 6d9e99736bef4837b63601149c3119eb11d723dc Mon Sep 17 00:00:00 2001
From: Nathan Glaser <glaserna@oregonstate.edu>
Date: Wed, 29 Apr 2026 19:16:50 -0700
Subject: [PATCH 20/20] back in black

---
 mcdc/transport/simulation.py | 2 +-
 mcdc/transport/technique.py  | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/mcdc/transport/simulation.py b/mcdc/transport/simulation.py
index cefc8c95a..5d3050573 100644
--- a/mcdc/transport/simulation.py
+++ b/mcdc/transport/simulation.py
@@ -350,7 +350,6 @@ def step_particle(particle_container, program, data):
     # Apply techniques
     # ==================================================================================
 
-
     # Skip if not alive
     if not particle["alive"]:
         return
@@ -366,6 +365,7 @@ def step_particle(particle_container, program, data):
     if simulation["forced_collisions"]["active"]:
         technique.forced_collision_roulette(particle_container, program, data)
 
+
 @njit
 def move_to_event(particle_container, simulation, data):
     settings = simulation["settings"]
diff --git a/mcdc/transport/technique.py b/mcdc/transport/technique.py
index f25a65b72..77d5bcaa4 100644
--- a/mcdc/transport/technique.py
+++ b/mcdc/transport/technique.py
@@ -210,7 +210,7 @@ def in_forced_collision_cell(particle_container, simulation, data):
     # skip if not a neutron
     if particle_container[0]["particle_type"] != PARTICLE_NEUTRON:
         return False
-    
+
     fc_object = simulation["forced_collisions"]
 
     # not active, dont need to query cells