Extend Kalbach Mann systematics to support incident photons

docs/source/io_formats/nuclear_data.rst

@@ -415,6 +415,7 @@ Kalbach-Mann
 
 :Object type: Group
 :Attributes: - **type** (*char[]*) -- 'kalbach-mann'
+             - **particle** (*char[]*) -- Type of particle
 :Datasets: - **energy** (*double[]*) -- Incoming energies at which distributions exist
 
              :Attributes:

include/openmc/secondary_kalbach.h

@@ -44,6 +44,7 @@ class KalbachMann : public AngleEnergy {
     xt::xtensor<double, 1> a;     //!< Parameterized function
   };
 
+  bool is_photon_;                      //!< Whether the projectile is a photon
   int n_region_;                        //!< Number of interpolation regions
   vector<int> breakpoints_;             //!< Breakpoints between regions
   vector<Interpolation> interpolation_; //!< Interpolation laws

openmc/data/data.py

@@ -274,13 +274,17 @@
 # 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
 
 # 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] = {}
 

openmc/data/kalbach_mann.py

@@ -9,7 +9,7 @@
 from openmc.stats import Tabular, Univariate, Discrete, Mixture
 from .function import Tabulated1D, INTERPOLATION_SCHEME
 from .angle_energy import AngleEnergy
-from .data import EV_PER_MEV
+from .data import EV_PER_MEV, NEUTRON_MASS_EV
 from .endf import get_list_record, get_tab2_record
 
 
@@ -204,13 +204,14 @@ def kalbach_slope(energy_projectile, energy_emitted, za_projectile,
         Kalbach-Mann slope given with the same format as ACE file.
 
     """
-    # TODO: develop for photons as projectile
-    # TODO: test for other particles than neutron
-    if za_projectile != 1:
-        raise NotImplementedError(
-            "Developed and tested for neutron projectile only."
-        )
 
+    if za_projectile == 0:
+        # Calculate slope for photons using Eq. 3 in doi:10.1080/18811248.1995.9731830
+        # or ENDF-6 Formats Manual section 6.2.3.2
+        slope_n = kalbach_slope(energy_projectile, energy_emitted, 1,
+                  za_emitted, za_target)
+        return slope_n * np.sqrt(0.5*energy_projectile/NEUTRON_MASS_EV)*np.minimum(4,np.maximum(1,9.3/np.sqrt(energy_emitted/EV_PER_MEV))) 
+
     # Special handling of elemental carbon
     if za_emitted == 6000:
         za_emitted = 6012
@@ -268,6 +269,8 @@ class KalbachMann(AngleEnergy):
     slope : Iterable of openmc.data.Tabulated1D
         Kalbach-Chadwick angular distribution slope value 'a' as a function of
         outgoing energy for each incoming energy
+    particle : {'neutron', 'photon'}
+        Incident particle type, defaults to neutron
 
     Attributes
     ----------
@@ -285,18 +288,21 @@ class KalbachMann(AngleEnergy):
     slope : Iterable of openmc.data.Tabulated1D
         Kalbach-Chadwick angular distribution slope value 'a' as a function of
         outgoing energy for each incoming energy
+    particle : {'neutron', 'photon'} 
+        incident particle type, defaults to neutron        
 
     """
 
     def __init__(self, breakpoints, interpolation, energy, energy_out,
-                 precompound, slope):
+                 precompound, slope, particle = 'neutron'):
         super().__init__()
         self.breakpoints = breakpoints
         self.interpolation = interpolation
         self.energy = energy
         self.energy_out = energy_out
         self.precompound = precompound
         self.slope = slope
+        self.particle = particle
 
     @property
     def breakpoints(self):
@@ -317,6 +323,15 @@ def interpolation(self, interpolation):
         cv.check_type('Kalbach-Mann interpolation', interpolation,
                       Iterable, Integral)
         self._interpolation = interpolation
+
+    @property
+    def particle(self):
+        return self._particle
+
+    @particle.setter
+    def particle(self, particle):
+        cv.check_value('Kalbach-Mann incident particle', particle, ['neutron', 'photon'])
+        self._particle = particle
 
     @property
     def energy(self):
@@ -367,6 +382,7 @@ def to_hdf5(self, group):
 
         """
         group.attrs['type'] = np.bytes_('kalbach-mann')
+        group.attrs['particle'] = np.bytes_(self.particle)
 
         dset = group.create_dataset('energy', data=self.energy)
         dset.attrs['interpolation'] = np.vstack((self.breakpoints,
@@ -436,6 +452,7 @@ def from_hdf5(cls, group):
             Kalbach-Mann energy distribution
 
         """
+        particle = group.attrs.get("particle", b"neutron").decode()
         interp_data = group['energy'].attrs['interpolation']
         energy_breakpoints = interp_data[0, :]
         energy_interpolation = interp_data[1, :]
@@ -491,7 +508,7 @@ def from_hdf5(cls, group):
             slope.append(km_a)
 
         return cls(energy_breakpoints, energy_interpolation,
-                   energy, energy_out, precompound, slope)
+                   energy, energy_out, precompound, slope, particle = particle)
 
     @classmethod
     def from_ace(cls, ace, idx, ldis):
@@ -514,6 +531,7 @@ def from_ace(cls, ace, idx, ldis):
             Kalbach-Mann energy-angle distribution
 
         """
+        particle = {'u':'photon', 'c':'neutron'}[ace.data_type.value]
         # Read number of interpolation regions and incoming energies
         n_regions = int(ace.xss[idx])
         n_energy_in = int(ace.xss[idx + 1 + 2*n_regions])
@@ -586,10 +604,10 @@ def from_ace(cls, ace, idx, ldis):
             km_r.append(Tabulated1D(data[0], data[3]))
             km_a.append(Tabulated1D(data[0], data[4]))
 
-        return cls(breakpoints, interpolation, energy, energy_out, km_r, km_a)
+        return cls(breakpoints, interpolation, energy, energy_out, km_r, km_a, particle = particle)
 
     @classmethod
-    def from_endf(cls, file_obj, za_emitted, za_target, projectile_mass):
+    def from_endf(cls, file_obj, za_emitted, za_target, za_projectile):
         """Generate Kalbach-Mann distribution from an ENDF evaluation.
 
         If the projectile is a neutron, the slope is calculated when it is
@@ -606,21 +624,16 @@ def from_endf(cls, file_obj, za_emitted, za_target, projectile_mass):
             ZA identifier of the emitted particle
         za_target : int
             ZA identifier of the target
-        projectile_mass : float
-            Mass of the projectile
-
-        Warns
-        -----
-        UserWarning
-            If the mass of the projectile is not equal to 1 (other than
-            a neutron), the slope is not calculated and set to 0 if missing.
+        za_projectile : int
+            ZA identifier of the projectile
 
         Returns
         -------
         openmc.data.KalbachMann
             Kalbach-Mann energy-angle distribution
 
         """
+        particle = {0: 'photon', 1: 'neutron'}[za_projectile]        
         params, tab2 = get_tab2_record(file_obj)
         lep = params[3]
         ne = params[5]
@@ -654,31 +667,19 @@ def from_endf(cls, file_obj, za_emitted, za_target, projectile_mass):
                 a_i = values[:, 3]
                 calculated_slope.append(False)
             else:
-                # Check if the projectile is not a neutron
-                if not np.isclose(projectile_mass, 1.0, atol=1.0e-12, rtol=0.):
-                    warn(
-                        "Kalbach-Mann slope calculation is only available with "
-                        "neutrons as projectile. Slope coefficients are set to 0."
-                    )
-                    a_i = np.zeros_like(r_i)
-                    calculated_slope.append(False)
-
-                else:
-                    # TODO: retrieve ZA of the projectile
-                    za_projectile = 1
-                    a_i = [kalbach_slope(energy_projectile=energy[i],
-                                         energy_emitted=e,
-                                         za_projectile=za_projectile,
-                                         za_emitted=za_emitted,
-                                         za_target=za_target)
-                           for e in eout_i]
-                    calculated_slope.append(True)
+                a_i = [kalbach_slope(energy_projectile=energy[i],
+                                     energy_emitted=e,
+                                     za_projectile=za_projectile,
+                                     za_emitted=za_emitted,
+                                     za_target=za_target)
+                       for e in eout_i]
+                calculated_slope.append(True)
 
             precompound.append(Tabulated1D(eout_i, r_i))
             slope.append(Tabulated1D(eout_i, a_i))
 
         km_distribution = cls(tab2.breakpoints, tab2.interpolation, energy,
-                              energy_out, precompound, slope)
+                              energy_out, precompound, slope, particle = particle)
 
         # List of bool to indicate slope calculation by OpenMC
         km_distribution._calculated_slope = calculated_slope

openmc/data/reaction.py

@@ -87,6 +87,8 @@ def _get_products(ev, mt):
         is not defined in the 'center-of-mass' system. The breakup logic
         is not implemented which can lead to this error being raised while
         the definition of the product is correct.
+    NotImplementedError
+        When the projectile is not a neutron and not a photon.
 
     Returns
     -------
@@ -163,11 +165,16 @@ def _get_products(ev, mt):
                     )
 
                 zat = ev.target["atomic_number"] * 1000 + ev.target["mass_number"]
-                projectile_mass = ev.projectile["mass"]
+                if ev.projectile['mass'] == 0.0:
+                    za_projectile = 0
+                elif np.isclose(ev.projectile['mass'], 1.0, atol=1.0e-12, rtol=0.):
+                    za_projectile = 1
+                else:
+                    raise NotImplementedError('Unknown projectile')
                 p.distribution = [KalbachMann.from_endf(file_obj,
                                                         za,
                                                         zat,
-                                                        projectile_mass)]
+                                                        za_projectile)]
 
         elif law == 2:
             # Discrete two-body scattering

src/secondary_kalbach.cpp

@@ -22,6 +22,14 @@ namespace openmc {
 
 KalbachMann::KalbachMann(hid_t group)
 {
+  is_photon_ = false;
+  // Check if projectile is a photon
+  if (attribute_exists(group, "particle")) {
+    std::string temp;
+    read_attribute(group, "particle", temp);
+    if (temp == "photon")
+      is_photon_ = true;
+  }
   // Open incoming energy dataset
   hid_t dset = open_dataset(group, "energy");
 
@@ -229,12 +237,19 @@ void KalbachMann::sample(
   }
 
   // Sampled correlated angle from Kalbach-Mann parameters
-  if (prn(seed) > km_r) {
-    double T = uniform_distribution(-1., 1., seed) * std::sinh(km_a);
-    mu = std::log(T + std::sqrt(T * T + 1.0)) / km_a;
+  if (is_photon_) {
+    double T = uniform_distribution(0., 1., seed);
+    double sinha = std::sinh(km_a);
+    mu = sinha * ((1 + T) - 2 * km_r) /
+         (sinha * T + std::cosh(km_a) - std::exp(km_a * T));
   } else {
-    double r1 = prn(seed);
-    mu = std::log(r1 * std::exp(km_a) + (1.0 - r1) * std::exp(-km_a)) / km_a;
+    if (prn(seed) > km_r) {
+      double T = uniform_distribution(-1., 1., seed) * std::sinh(km_a);
+      mu = std::log(T + std::sqrt(T * T + 1.0)) / km_a;
+    } else {
+      double r1 = prn(seed);
+      mu = std::log(r1 * std::exp(km_a) + (1.0 - r1) * std::exp(-km_a)) / km_a;
+    }
   }
 }
 

tests/unit_tests/test_data_kalbach_mann.py

@@ -103,16 +103,13 @@ def test_kalbach_slope():
     energy_projectile = 10.2  # [eV]
     energy_emitted = 5.4  # [eV]
 
-    # Check that NotImplementedError is raised if the projectile is not
-    # a neutron
-    with pytest.raises(NotImplementedError):
-        kalbach_slope(
-            energy_projectile=energy_projectile,
-            energy_emitted=energy_emitted,
-            za_projectile=1000,
-            za_emitted=1,
-            za_target=6012
-        )
+    kalbach_slope(
+        energy_projectile=energy_projectile,
+        energy_emitted=energy_emitted,
+        za_projectile=0,
+        za_emitted=1,
+        za_target=6012
+    )
 
     assert kalbach_slope(
         energy_projectile=energy_projectile,