From 80bef50aabf56e8af968e1958cf42bf5d00fbe7a Mon Sep 17 00:00:00 2001
From: GuySten <guyste@post.bgu.ac.il>
Date: Tue, 9 Dec 2025 21:12:08 +0200
Subject: [PATCH 01/14] extend level scattering to support incident photon

---
 include/openmc/distribution_energy.h |  5 +-
 openmc/data/energy_distribution.py   | 88 ++++++++++++++++++----------
 openmc/data/reaction.py              |  7 +--
 src/distribution_energy.cpp          | 29 ++++++++-
 4 files changed, 87 insertions(+), 42 deletions(-)

diff --git a/include/openmc/distribution_energy.h b/include/openmc/distribution_energy.h
index 9b08ed039dc..b6c8ad2a446 100644
--- a/include/openmc/distribution_energy.h
+++ b/include/openmc/distribution_energy.h
@@ -61,8 +61,9 @@ class LevelInelastic : public EnergyDistribution {
   double sample(double E, uint64_t* seed) const override;
 
 private:
-  double threshold_;  //!< Energy threshold in lab, (A + 1)/A * |Q|
-  double mass_ratio_; //!< (A/(A+1))^2
+  double a_;
+  double b_;
+  double c_;
 };
 
 //===============================================================================
diff --git a/openmc/data/energy_distribution.py b/openmc/data/energy_distribution.py
index 069ab1b9beb..463cd85daac 100644
--- a/openmc/data/energy_distribution.py
+++ b/openmc/data/energy_distribution.py
@@ -1,5 +1,6 @@
 from abc import ABC, abstractmethod
 from collections.abc import Iterable
+from math import sqrt
 from numbers import Integral, Real
 from warnings import warn
 
@@ -897,42 +898,54 @@ class LevelInelastic(EnergyDistribution):
 
     Parameters
     ----------
-    threshold : float
-        Energy threshold in the laboratory system, :math:`(A + 1)/A * |Q|`
-    mass_ratio : float
-        :math:`(A/(A + 1))^2`
+    q_value : float
+        Q-value of the reaction.
+    mass : float
+        mass of the nucleus in units of neutron mass.
+    particle : ParticleType
+        incident particle type, defaults to neutron
 
     Attributes
     ----------
-    threshold : float
-        Energy threshold in the laboratory system, :math:`(A + 1)/A * |Q|`
-    mass_ratio : float
-        :math:`(A/(A + 1))^2`
-
+    q_value : float
+        Q-value of the reaction.
+    mass : float
+        mass of the nucleus in units of neutron mass.
+    particle : ParticleType
+        incident particle type.
     """
 
-    def __init__(self, threshold, mass_ratio):
+    def __init__(self, q_value, mass, particle = 'neutron'):
         super().__init__()
-        self.threshold = threshold
-        self.mass_ratio = mass_ratio
+        self.q_value = q_value
+        self.mass = mass
+        self.particle = particle
 
     @property
-    def threshold(self):
-        return self._threshold
+    def q_value(self):
+        return self._q_value
 
-    @threshold.setter
-    def threshold(self, threshold):
-        cv.check_type('level inelastic threhsold', threshold, Real)
-        self._threshold = threshold
+    @q_value.setter
+    def q_value(self, q_value):
+        cv.check_type('level inelastic q_value', q_value, Real)
+        self._q_value = q_value
 
     @property
-    def mass_ratio(self):
-        return self._mass_ratio
-
-    @mass_ratio.setter
-    def mass_ratio(self, mass_ratio):
-        cv.check_type('level inelastic mass ratio', mass_ratio, Real)
-        self._mass_ratio = mass_ratio
+    def mass(self):
+        return self._mass
+
+    @mass.setter
+    def mass(self, mass):
+        cv.check_type('level inelastic mass', mass, Real)
+        self._mass = mass
+        
+    def particle(self):
+        return self._particle
+    
+    @particle.setter
+    def particle(self, particle):
+        cv.check_type('product particle type', particle, str)
+        self._particle = particle
 
     def to_hdf5(self, group):
         """Write distribution to an HDF5 group
@@ -945,8 +958,9 @@ def to_hdf5(self, group):
         """
 
         group.attrs['type'] = np.bytes_('level')
-        group.attrs['threshold'] = self.threshold
-        group.attrs['mass_ratio'] = self.mass_ratio
+        group.attrs['q_value'] = self.q_value
+        group.attrs['mass'] = self.mass
+        group.attrs['particle'] = np.bytes_(self.particle)
 
     @classmethod
     def from_hdf5(cls, group):
@@ -963,9 +977,18 @@ def from_hdf5(cls, group):
             Level inelastic scattering distribution
 
         """
-        threshold = group.attrs['threshold']
-        mass_ratio = group.attrs['mass_ratio']
-        return cls(threshold, mass_ratio)
+        #backwards compatible read:
+        if 'threshold' in group.attrs:    
+            threshold = group.attrs['threshold']
+            mass_ratio = group.attrs['mass_ratio']
+            mass = 1.0/(1.0/sqrt(mass_ratio)-1.0)
+            q_value = -threshold * sqrt(mass_ratio)
+            return cls(q_value, mass)
+            
+        q_value = group.attrs['q_value']
+        mass = group.attrs['mass']
+        particle = group.attrs['particle']
+        return cls(q_value, mass, particle)
 
     @classmethod
     def from_ace(cls, ace, idx):
@@ -984,9 +1007,12 @@ def from_ace(cls, ace, idx):
             Level inelastic scattering distribution
 
         """
+        particle = {'u':'photon', 'c': 'neutron'}[ace.data_type.value]
         threshold = ace.xss[idx]*EV_PER_MEV
         mass_ratio = ace.xss[idx + 1]
-        return cls(threshold, mass_ratio)
+        mass = 1.0/(1.0/sqrt(mass_ratio)-1.0) if particle == 'neutron' else 1.0-1.0/mass_ratio
+        q_value = -threshold * sqrt(mass_ratio) if particle == 'neutron' else -threshold
+        return cls(threshold, mass_ratio, particle = particle)
 
 
 class ContinuousTabular(EnergyDistribution):
diff --git a/openmc/data/reaction.py b/openmc/data/reaction.py
index 65b59582cf8..4160bf71add 100644
--- a/openmc/data/reaction.py
+++ b/openmc/data/reaction.py
@@ -1200,12 +1200,7 @@ def from_endf(cls, ev, mt):
                 # since it can be calculated analytically. Here we determine the
                 # necessary parameters to create a LevelInelastic object
                 dist = UncorrelatedAngleEnergy()
-
-                A = ev.target['mass']
-                threshold = (A + 1.)/A*abs(rx.q_value)
-                mass_ratio = (A/(A + 1.))**2
-                dist.energy = LevelInelastic(threshold, mass_ratio)
-
+                dist.energy = LevelInelastic(rx.q_value, ev.target['mass'])
                 neutron.distribution.append(dist)
 
             if (4, mt) in ev.section:
diff --git a/src/distribution_energy.cpp b/src/distribution_energy.cpp
index a4a5ce9e1b6..d5ad2b3a89c 100644
--- a/src/distribution_energy.cpp
+++ b/src/distribution_energy.cpp
@@ -1,6 +1,7 @@
 #include "openmc/distribution_energy.h"
 
 #include <algorithm> // for max, min, copy, move
+#include <cmath>     // for sqrt, abs
 #include <cstddef>   // for size_t
 #include <iterator>  // for back_inserter
 
@@ -41,13 +42,35 @@ double DiscretePhoton::sample(double E, uint64_t* seed) const
 
 LevelInelastic::LevelInelastic(hid_t group)
 {
-  read_attribute(group, "threshold", threshold_);
-  read_attribute(group, "mass_ratio", mass_ratio_);
+  // for backwards compatibility:
+  if attribute_exists (group, "mass_ratio") {
+    read_attribute(group, "threshold", b_);
+    read_attribute(group, "mass_ratio", a_);
+    c_ = 0.0;
+  } else {
+    double A, Q;
+    std::string temp;
+    read_attribute(group, "mass", A);
+    read_attribute(group, "q_value", Q);
+    read_attribute(group, "particle", temp);
+    auto particle = str_to_particle_type(temp);
+    if (particle == ParticleType::neutron) {
+      a_ = (A / (A + 1.0)) * (A / (A + 1.0));
+      b_ = (A + 1.0) A * std::abs(Q);
+      c_ = 0.0;
+    } else if (particle_ == ParticleType::photon) {
+      a_ = (A - 1.0) / A;
+      b_ = -Q;
+      c_ = 1.0 / (2.0 * NEUTRON_MASS_EV * A);
+    } else {
+      fatal_error("Unrecognized particle: " + particle);
+    }
+  }
 }
 
 double LevelInelastic::sample(double E, uint64_t* seed) const
 {
-  return mass_ratio_ * (E - threshold_);
+  return a_ * (E - b_ - c_ * (E * E));
 }
 
 //==============================================================================

From ebf9881e21d1277e2db436086f6cf79abc2908a8 Mon Sep 17 00:00:00 2001
From: GuySten <guyste@post.bgu.ac.il>
Date: Tue, 9 Dec 2025 22:19:52 +0200
Subject: [PATCH 02/14] some fixes

---
 src/distribution_energy.cpp | 7 ++++---
 1 file changed, 4 insertions(+), 3 deletions(-)

diff --git a/src/distribution_energy.cpp b/src/distribution_energy.cpp
index d5ad2b3a89c..7a22d40c47e 100644
--- a/src/distribution_energy.cpp
+++ b/src/distribution_energy.cpp
@@ -10,6 +10,7 @@
 #include "openmc/endf.h"
 #include "openmc/hdf5_interface.h"
 #include "openmc/math_functions.h"
+#include "openmc/particle.h"
 #include "openmc/random_dist.h"
 #include "openmc/random_lcg.h"
 #include "openmc/search.h"
@@ -43,7 +44,7 @@ double DiscretePhoton::sample(double E, uint64_t* seed) const
 LevelInelastic::LevelInelastic(hid_t group)
 {
   // for backwards compatibility:
-  if attribute_exists (group, "mass_ratio") {
+  if (attribute_exists(group, "mass_ratio")) {
     read_attribute(group, "threshold", b_);
     read_attribute(group, "mass_ratio", a_);
     c_ = 0.0;
@@ -56,9 +57,9 @@ LevelInelastic::LevelInelastic(hid_t group)
     auto particle = str_to_particle_type(temp);
     if (particle == ParticleType::neutron) {
       a_ = (A / (A + 1.0)) * (A / (A + 1.0));
-      b_ = (A + 1.0) A * std::abs(Q);
+      b_ = (A + 1.0) / A * std::abs(Q);
       c_ = 0.0;
-    } else if (particle_ == ParticleType::photon) {
+    } else if (particle == ParticleType::photon) {
       a_ = (A - 1.0) / A;
       b_ = -Q;
       c_ = 1.0 / (2.0 * NEUTRON_MASS_EV * A);

From a25586d3a05d124520d0fa130f46de644748d4a6 Mon Sep 17 00:00:00 2001
From: GuySten <guyste@post.bgu.ac.il>
Date: Tue, 9 Dec 2025 22:20:30 +0200
Subject: [PATCH 03/14] another fix

---
 src/distribution_energy.cpp | 1 +
 1 file changed, 1 insertion(+)

diff --git a/src/distribution_energy.cpp b/src/distribution_energy.cpp
index 7a22d40c47e..bf3cb13182e 100644
--- a/src/distribution_energy.cpp
+++ b/src/distribution_energy.cpp
@@ -7,6 +7,7 @@
 
 #include "xtensor/xview.hpp"
 
+#include "openmc/constants.h"
 #include "openmc/endf.h"
 #include "openmc/hdf5_interface.h"
 #include "openmc/math_functions.h"

From 7fb6adc167dc9a6e57b7ca0de86fef9f1d01bea2 Mon Sep 17 00:00:00 2001
From: GuySten <guyste@post.bgu.ac.il>
Date: Tue, 9 Dec 2025 22:28:50 +0200
Subject: [PATCH 04/14] fix typo

---
 src/distribution_energy.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/distribution_energy.cpp b/src/distribution_energy.cpp
index bf3cb13182e..5d583b805ee 100644
--- a/src/distribution_energy.cpp
+++ b/src/distribution_energy.cpp
@@ -63,7 +63,7 @@ LevelInelastic::LevelInelastic(hid_t group)
     } else if (particle == ParticleType::photon) {
       a_ = (A - 1.0) / A;
       b_ = -Q;
-      c_ = 1.0 / (2.0 * NEUTRON_MASS_EV * A);
+      c_ = 1.0 / (2.0 * MASS_NEUTRON_EV * A);
     } else {
       fatal_error("Unrecognized particle: " + particle);
     }

From f65f5b473e2ed5631f2d51add67b7b6a409e6abd Mon Sep 17 00:00:00 2001
From: GuySten <guyste@post.bgu.ac.il>
Date: Tue, 9 Dec 2025 22:29:27 +0200
Subject: [PATCH 05/14] another fix

---
 src/distribution_energy.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/distribution_energy.cpp b/src/distribution_energy.cpp
index 5d583b805ee..8385393c20a 100644
--- a/src/distribution_energy.cpp
+++ b/src/distribution_energy.cpp
@@ -65,7 +65,7 @@ LevelInelastic::LevelInelastic(hid_t group)
       b_ = -Q;
       c_ = 1.0 / (2.0 * MASS_NEUTRON_EV * A);
     } else {
-      fatal_error("Unrecognized particle: " + particle);
+      fatal_error("Unrecognized particle: " + temp);
     }
   }
 }

From 123f1b7f77b3354d331fddf0e4cca501434d77b1 Mon Sep 17 00:00:00 2001
From: GuySten <guyste@post.bgu.ac.il>
Date: Tue, 9 Dec 2025 22:46:19 +0200
Subject: [PATCH 06/14] fix typo

---
 openmc/data/energy_distribution.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/openmc/data/energy_distribution.py b/openmc/data/energy_distribution.py
index 463cd85daac..b813393c0b3 100644
--- a/openmc/data/energy_distribution.py
+++ b/openmc/data/energy_distribution.py
@@ -938,7 +938,8 @@ def mass(self):
     def mass(self, mass):
         cv.check_type('level inelastic mass', mass, Real)
         self._mass = mass
-        
+    
+    @property    
     def particle(self):
         return self._particle
     

From 17798bf25c52c71b9abc6266c972e2d8a5715524 Mon Sep 17 00:00:00 2001
From: GuySten <guyste@post.bgu.ac.il>
Date: Tue, 9 Dec 2025 22:50:09 +0200
Subject: [PATCH 07/14] update documentation

---
 docs/source/io_formats/nuclear_data.rst | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/docs/source/io_formats/nuclear_data.rst b/docs/source/io_formats/nuclear_data.rst
index 8174108e1ae..90519d445f8 100644
--- a/docs/source/io_formats/nuclear_data.rst
+++ b/docs/source/io_formats/nuclear_data.rst
@@ -584,9 +584,9 @@ Level Inelastic
 
 :Object type: Group
 :Attributes: - **type** (*char[]*) -- 'level'
-             - **threshold** (*double*) -- Energy threshold in the laboratory
-               system in eV
-             - **mass_ratio** (*double*) -- :math:`(A/(A + 1))^2`
+             - **q_value** (*double*) -- :math:`Q`-value of the reaction.
+             - **mass** (*double*) -- nucleus mass :math:`A` relative to neutron rest mass.  
+             - **particle** (*char[]*) -- incident particle name
 
 Continuous Tabular
 ------------------

From 8ceb9f400802fd6e4a801e184694db78620c516a Mon Sep 17 00:00:00 2001
From: GuySten <guyste@post.bgu.ac.il>
Date: Wed, 10 Dec 2025 00:45:56 +0200
Subject: [PATCH 08/14] fix type conversion

---
 openmc/data/energy_distribution.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/openmc/data/energy_distribution.py b/openmc/data/energy_distribution.py
index b813393c0b3..95be853c8a6 100644
--- a/openmc/data/energy_distribution.py
+++ b/openmc/data/energy_distribution.py
@@ -988,7 +988,7 @@ def from_hdf5(cls, group):
             
         q_value = group.attrs['q_value']
         mass = group.attrs['mass']
-        particle = group.attrs['particle']
+        particle = group.attrs['particle'].decode()
         return cls(q_value, mass, particle)
 
     @classmethod

From 5ce1ffa1cc51d853a5345fb3f22bd8604092bf8a Mon Sep 17 00:00:00 2001
From: GuySten <guyste@post.bgu.ac.il>
Date: Sat, 13 Dec 2025 19:39:05 +0200
Subject: [PATCH 09/14] WIP

---
 docs/source/methods/neutron_physics.rst |  2 +-
 docs/source/methods/photon_physics.rst  | 31 ++++++++++++++++++++++---
 openmc/data/data.py                     |  4 ++++
 openmc/data/energy_distribution.py      | 21 +++++++++++++----
 src/distribution_energy.cpp             |  4 ++--
 5 files changed, 52 insertions(+), 10 deletions(-)

diff --git a/docs/source/methods/neutron_physics.rst b/docs/source/methods/neutron_physics.rst
index 2b797e3dbcf..aa45f51ac39 100644
--- a/docs/source/methods/neutron_physics.rst
+++ b/docs/source/methods/neutron_physics.rst
@@ -568,7 +568,7 @@ and the incoming energy:
 .. math::
     :label: level-scattering
 
-    E' = \left ( \frac{A}{A+1} \right )^2 \left ( E - \frac{A + 1}{A} Q \right )
+    E' = \left ( \frac{A}{A+1} \right )^2 \left ( E - \frac{A + 1}{A} |Q| \right )
 
 where :math:`A` is the mass of the target nucleus measured in neutron masses.
 
diff --git a/docs/source/methods/photon_physics.rst b/docs/source/methods/photon_physics.rst
index d2bd3ac7609..fb5ab594384 100644
--- a/docs/source/methods/photon_physics.rst
+++ b/docs/source/methods/photon_physics.rst
@@ -16,9 +16,9 @@ de-excitation of these atoms can result in the emission of electrons and
 photons. Electrons themselves also can produce photons by means of
 bremsstrahlung radiation.
 
--------------------
-Photon Interactions
--------------------
+--------------------------
+Photon Atomic Interactions
+--------------------------
 
 Coherent (Rayleigh) Scattering
 ------------------------------
@@ -598,6 +598,29 @@ The inverse transform method is used to sample :math:`\mu_{-}` and
 The azimuthal angles for the electron and positron are sampled independently
 and uniformly on :math:`[0, 2\pi)`.
 
+-------------------------
+Photonuclear Interactions
+-------------------------
+
+Photonuclear reactions occur when a nucleus absorbs a high-energy photon (gamma-ray), leading to a change in the nucleus. 
+The absorption of the photon excites the nucleus, and the energy provided by the photon can result in the emission of particles 
+like neutrons, protons, or other light nuclei. Because photonuclear interactions are determined by nuclear rather than atomic properties, 
+their data libraries must be specific to each isotope.
+
+Inelastic Level Scattering
+--------------------------
+
+It can be shown (see Wattenberg_) that in inelastic level scattering, the outgoing
+energy of the neutron :math:`E'` can be related to the Q-value of the reaction
+and the incoming energy of the photon:
+
+.. math::
+    :label: photonuclear-level-scattering
+
+    E' = \left ( \frac{A-1}{A} \right ) \left ( E - |Q| - \frac{E^2}{2 m_n \left (A-1 \right )} \right )
+
+where :math:`A` is the mass of the target nucleus measured in neutron masses, math:`m_n` is the neutron mass in eV.
+
 -------------------
 Secondary Processes
 -------------------
@@ -734,3 +757,5 @@ emitted photon.
 .. _Kaltiaisenaho: https://aaltodoc.aalto.fi/bitstream/handle/123456789/21004/master_Kaltiaisenaho_Toni_2016.pdf
 
 .. _Salvat: https://doi.org/10.1787/32da5043-en
+
+.. _Wattenberg: https://doi.org/10.1103/PhysRev.71.497
diff --git a/openmc/data/data.py b/openmc/data/data.py
index 5ecadd37be0..e075ebbe145 100644
--- a/openmc/data/data.py
+++ b/openmc/data/data.py
@@ -274,6 +274,7 @@
 # Unit conversions
 EV_PER_MEV = 1.0e6
 JOULE_PER_EV = 1.602176634e-19
+EV_PER_AMU = 931.49410372e6 # eV/c^2 per amu
 
 # Avogadro's constant
 AVOGADRO = 6.02214076e23
@@ -281,6 +282,9 @@
 # Neutron mass in units of amu
 NEUTRON_MASS = 1.00866491595
 
+# Neutron mass in units of eV/c^2
+NEUTRON_MASS_EV = NEUTRON_MASS*EV_PER_AMU
+
 # Used in atomic_mass function as a cache
 _ATOMIC_MASS: dict[str, float] = {}
 
diff --git a/openmc/data/energy_distribution.py b/openmc/data/energy_distribution.py
index 95be853c8a6..0ebfaac1d47 100644
--- a/openmc/data/energy_distribution.py
+++ b/openmc/data/energy_distribution.py
@@ -1,6 +1,6 @@
 from abc import ABC, abstractmethod
 from collections.abc import Iterable
-from math import sqrt
+from math import sqrt, isclose
 from numbers import Integral, Real
 from warnings import warn
 
@@ -9,7 +9,7 @@
 import openmc.checkvalue as cv
 from openmc.mixin import EqualityMixin
 from openmc.stats.univariate import Univariate, Tabular, Discrete, Mixture
-from .data import EV_PER_MEV
+from .data import EV_PER_MEV, NEUTRON_MASS_EV
 from .endf import get_tab1_record, get_tab2_record
 from .function import Tabulated1D, INTERPOLATION_SCHEME
 
@@ -902,7 +902,7 @@ class LevelInelastic(EnergyDistribution):
         Q-value of the reaction.
     mass : float
         mass of the nucleus in units of neutron mass.
-    particle : ParticleType
+    particle : {'neutron', 'photon'}
         incident particle type, defaults to neutron
 
     Attributes
@@ -945,8 +945,18 @@ def particle(self):
     
     @particle.setter
     def particle(self, particle):
-        cv.check_type('product particle type', particle, str)
+        cv.check_type('product particle type', particle, ['neutron', 'photon'])
         self._particle = particle
+        
+    @property
+    def threshold(self):
+        A = self.mass
+        Q = self.q_value
+        if particle == 'neutron':
+            return (A+1.0)/A*abs(Q)
+        else:
+            b = 2*NEUTRON_MASS_EV*(A-1)
+            return (b-sqrt(b**2-4.0*b*abs(Q)))/2.0
 
     def to_hdf5(self, group):
         """Write distribution to an HDF5 group
@@ -1012,6 +1022,9 @@ def from_ace(cls, ace, idx):
         threshold = ace.xss[idx]*EV_PER_MEV
         mass_ratio = ace.xss[idx + 1]
         mass = 1.0/(1.0/sqrt(mass_ratio)-1.0) if particle == 'neutron' else 1.0-1.0/mass_ratio
+        if not isclose(ace.atomic_weight_ratio, mass):
+            warn("Level inelastic distribution mass parameter does not match ace table mass.")
+            mass = ace.atomic_weight_ratio
         q_value = -threshold * sqrt(mass_ratio) if particle == 'neutron' else -threshold
         return cls(threshold, mass_ratio, particle = particle)
 
diff --git a/src/distribution_energy.cpp b/src/distribution_energy.cpp
index 8385393c20a..8450826447c 100644
--- a/src/distribution_energy.cpp
+++ b/src/distribution_energy.cpp
@@ -62,8 +62,8 @@ LevelInelastic::LevelInelastic(hid_t group)
       c_ = 0.0;
     } else if (particle == ParticleType::photon) {
       a_ = (A - 1.0) / A;
-      b_ = -Q;
-      c_ = 1.0 / (2.0 * MASS_NEUTRON_EV * A);
+      b_ = std::abs(Q);
+      c_ = 1.0 / (2.0 * MASS_NEUTRON_EV * (A - 1.0));
     } else {
       fatal_error("Unrecognized particle: " + temp);
     }

From bf01db4c93da4ab2e580e931d2046c5418c631d0 Mon Sep 17 00:00:00 2001
From: GuySten <guyste@post.bgu.ac.il>
Date: Sat, 13 Dec 2025 20:29:31 +0200
Subject: [PATCH 10/14] fix typo

---
 openmc/data/energy_distribution.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/openmc/data/energy_distribution.py b/openmc/data/energy_distribution.py
index 0ebfaac1d47..e356b7881d5 100644
--- a/openmc/data/energy_distribution.py
+++ b/openmc/data/energy_distribution.py
@@ -1026,7 +1026,7 @@ def from_ace(cls, ace, idx):
             warn("Level inelastic distribution mass parameter does not match ace table mass.")
             mass = ace.atomic_weight_ratio
         q_value = -threshold * sqrt(mass_ratio) if particle == 'neutron' else -threshold
-        return cls(threshold, mass_ratio, particle = particle)
+        return cls(q_value, mass, particle = particle)
 
 
 class ContinuousTabular(EnergyDistribution):

From ad7f3c60fa446e2d90c73f7ae8df2afa59862101 Mon Sep 17 00:00:00 2001
From: GuySten <guyste@post.bgu.ac.il>
Date: Sat, 13 Dec 2025 21:59:47 +0200
Subject: [PATCH 11/14] fix typo

---
 openmc/data/energy_distribution.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/openmc/data/energy_distribution.py b/openmc/data/energy_distribution.py
index e356b7881d5..249e17f8a36 100644
--- a/openmc/data/energy_distribution.py
+++ b/openmc/data/energy_distribution.py
@@ -945,7 +945,7 @@ def particle(self):
     
     @particle.setter
     def particle(self, particle):
-        cv.check_type('product particle type', particle, ['neutron', 'photon'])
+        cv.check_value('product particle type', particle, ['neutron', 'photon'])
         self._particle = particle
         
     @property

From 4af4576a230c5bfec33db19aefbaefd10df2900a Mon Sep 17 00:00:00 2001
From: GuySten <guyste@post.bgu.ac.il>
Date: Sat, 13 Dec 2025 22:13:11 +0200
Subject: [PATCH 12/14] added documentation and simplified math

---
 docs/source/methods/photon_physics.rst | 6 +++---
 include/openmc/distribution_energy.h   | 6 +++---
 openmc/data/energy_distribution.py     | 8 +++++---
 3 files changed, 11 insertions(+), 9 deletions(-)

diff --git a/docs/source/methods/photon_physics.rst b/docs/source/methods/photon_physics.rst
index fb5ab594384..5eafb439f07 100644
--- a/docs/source/methods/photon_physics.rst
+++ b/docs/source/methods/photon_physics.rst
@@ -16,9 +16,9 @@ de-excitation of these atoms can result in the emission of electrons and
 photons. Electrons themselves also can produce photons by means of
 bremsstrahlung radiation.
 
---------------------------
-Photon Atomic Interactions
---------------------------
+------------------------
+Photoatomic Interactions
+------------------------
 
 Coherent (Rayleigh) Scattering
 ------------------------------
diff --git a/include/openmc/distribution_energy.h b/include/openmc/distribution_energy.h
index b6c8ad2a446..ed9984503c8 100644
--- a/include/openmc/distribution_energy.h
+++ b/include/openmc/distribution_energy.h
@@ -61,9 +61,9 @@ class LevelInelastic : public EnergyDistribution {
   double sample(double E, uint64_t* seed) const override;
 
 private:
-  double a_;
-  double b_;
-  double c_;
+  double a_; //!< a coefficient of the scattering law
+  double b_; //!< a coefficient of the scattering law
+  double c_; //!< a coefficient of the scattering law
 };
 
 //===============================================================================
diff --git a/openmc/data/energy_distribution.py b/openmc/data/energy_distribution.py
index 249e17f8a36..ede97e98ae9 100644
--- a/openmc/data/energy_distribution.py
+++ b/openmc/data/energy_distribution.py
@@ -911,8 +911,10 @@ class LevelInelastic(EnergyDistribution):
         Q-value of the reaction.
     mass : float
         mass of the nucleus in units of neutron mass.
-    particle : ParticleType
+    particle : {'neutron', 'photon'}
         incident particle type.
+    threshold : float
+        Energy threshold in the laboratory system
     """
 
     def __init__(self, q_value, mass, particle = 'neutron'):
@@ -955,8 +957,8 @@ def threshold(self):
         if particle == 'neutron':
             return (A+1.0)/A*abs(Q)
         else:
-            b = 2*NEUTRON_MASS_EV*(A-1)
-            return (b-sqrt(b**2-4.0*b*abs(Q)))/2.0
+            b = NEUTRON_MASS_EV*(A-1)
+            return sqrt(b)*(sqrt(b)-sqrt(b-2*abs(Q)))
 
     def to_hdf5(self, group):
         """Write distribution to an HDF5 group

From cf76c5ffc5d1af354c177f268bd18e55d8b6715d Mon Sep 17 00:00:00 2001
From: GuySten <guyste@post.bgu.ac.il>
Date: Sat, 13 Dec 2025 22:16:29 +0200
Subject: [PATCH 13/14] improve docs

---
 include/openmc/distribution_energy.h | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/include/openmc/distribution_energy.h b/include/openmc/distribution_energy.h
index ed9984503c8..1ca9b905b35 100644
--- a/include/openmc/distribution_energy.h
+++ b/include/openmc/distribution_energy.h
@@ -61,9 +61,10 @@ class LevelInelastic : public EnergyDistribution {
   double sample(double E, uint64_t* seed) const override;
 
 private:
+  //! The scattering law is E_out = a*(E_in-b-c*E_in*E_in)
   double a_; //!< a coefficient of the scattering law
-  double b_; //!< a coefficient of the scattering law
-  double c_; //!< a coefficient of the scattering law
+  double b_; //!< b coefficient of the scattering law
+  double c_; //!< c coefficient of the scattering law
 };
 
 //===============================================================================

From 3fafce82c847f071b632dc23a382d0ffbe5fa26d Mon Sep 17 00:00:00 2001
From: GuySten <guyste@post.bgu.ac.il>
Date: Sat, 13 Dec 2025 23:15:52 +0200
Subject: [PATCH 14/14] improve docs

---
 docs/source/io_formats/nuclear_data.rst | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/docs/source/io_formats/nuclear_data.rst b/docs/source/io_formats/nuclear_data.rst
index 90519d445f8..ade764bb348 100644
--- a/docs/source/io_formats/nuclear_data.rst
+++ b/docs/source/io_formats/nuclear_data.rst
@@ -584,9 +584,9 @@ Level Inelastic
 
 :Object type: Group
 :Attributes: - **type** (*char[]*) -- 'level'
-             - **q_value** (*double*) -- :math:`Q`-value of the reaction.
-             - **mass** (*double*) -- nucleus mass :math:`A` relative to neutron rest mass.  
-             - **particle** (*char[]*) -- incident particle name
+             - **q_value** (*double*) -- Q value in eV
+             - **mass** (*double*) -- Nucleus mass A relative to neutron rest mass  
+             - **particle** (*char[]*) -- Incident particle name
 
 Continuous Tabular
 ------------------