Annihilate positrons below the positron energy cutoff

docs/source/io_formats/settings.rst

@@ -201,7 +201,9 @@ time.
     *Default*: 0.0
 
   :energy_positron:
-    The energy under which positrons will be killed.
+    The energy under which positrons will be killed. A killed positron
+    deposits its kinetic energy locally but still produces a pair of
+    annihilation photons.
 
     *Default*: 0.0
 

src/physics.cpp

@@ -471,13 +471,17 @@ void process_charged_secondary(
   Particle& p, Direction u, double E, ParticleType type)
 {
   int idx = type.transport_index();
-  if (idx == C_NONE || E < settings::energy_cutoff[idx])
+  if (idx == C_NONE)
     return;
 
-  if (settings::electron_treatment == ElectronTreatment::TTB) {
+  // Below the energy cutoff, the kinetic energy is deposited locally without
+  // producing bremsstrahlung
+  if (E >= settings::energy_cutoff[idx] &&
+      settings::electron_treatment == ElectronTreatment::TTB) {
     thick_target_bremsstrahlung(p, type, u, E);
   }
 
+  // A positron annihilates regardless of its kinetic energy
   if (type == ParticleType::positron()) {
     Direction photon_u = isotropic_direction(p.current_seed());
     p.create_secondary(

tests/unit_tests/test_energy_cutoff.py

@@ -61,3 +61,46 @@ def test_energy_cutoff(run_in_tmpdir):
     # Despite killing particles below the cutoff, the total heating should be
     # equal to the source energy
     assert heating[0] == pytest.approx(source_energy)
+
+
+def test_positron_cutoff_annihilation(run_in_tmpdir):
+    """Positrons below the energy cutoff must still annihilate"""
+
+    def run(positron_cutoff):
+        model = openmc.Model()
+        m = openmc.Material()
+        m.add_nuclide('Pb208', 1.0)
+        m.set_density('g/cm3', 11.35)
+        sph = openmc.Sphere(r=100.0, boundary_type='reflective')
+        model.geometry = openmc.Geometry([openmc.Cell(fill=m, region=-sph)])
+
+        model.settings.run_mode = 'fixed source'
+        model.settings.source = openmc.IndependentSource(
+            particle='photon',
+            energy=openmc.stats.Discrete([10.0e6], [1.0]),
+        )
+        model.settings.particles = 100
+        model.settings.batches = 2
+        # With local energy deposition, charged particles produce no
+        # bremsstrahlung, so the positron cutoff should have no effect
+        model.settings.electron_treatment = 'led'
+        if positron_cutoff is not None:
+            model.settings.cutoff = {'energy_positron': positron_cutoff}
+
+        tally = openmc.Tally()
+        tally.filters = [
+            openmc.EnergyFilter([0.0, 0.5e6, 0.52e6, 10.0e6]),
+            openmc.ParticleFilter(['photon'])
+        ]
+        tally.scores = ['flux']
+        model.tallies = [tally]
+        model.run(apply_tally_results=True)
+        return tally.mean.ravel()
+
+    # Every positron is below a cutoff above the source energy
+    flux_ref = run(None)
+    flux_cutoff = run(20.0e6)
+
+    # Annihilation photons contribute to the bin containing 511 keV
+    assert flux_ref[1] > 0.0
+    assert flux_cutoff == pytest.approx(flux_ref)