Add DRESS Monte Carlo integration for DT and DD thermal neutron spectra

pyproject.toml

@@ -40,6 +40,9 @@ docs = [
     "myst-parser>=2.0",
     "nbsphinx>=0.9",
 ]
+dress = [
+    "pydress"
+]
 
 [tool.setuptools.package-data]
 "NeSST.data" = ["*.dat", "*.txt","*.csv","*.json"]

src/NeSST/__init__.py

@@ -1,6 +1,10 @@
 ##################
 from .constants import *
 from .core import *
-from .dress_interface import *
+
+try:
+    from .dress_interface import *
+except ImportError:
+    pass
 from .fitting import *
 from .time_of_flight import *

src/NeSST/core.py

@@ -99,6 +99,41 @@ def QBallabio(Ein: npt.NDArray, mean: float, variance: float) -> npt.NDArray:
     return spec
 
 
+def QDress_DT(Ein: npt.NDArray, T_D: float, T_T: float | None = None, n_samples: int = int(1e6)) -> npt.NDArray:
+    """Calculates the DT primary neutron spectrum using the DRESS Monte Carlo code.
+
+    Args:
+        Ein (numpy.array): array of energy bin edges on which to compute the spectrum (eV)
+        T_D (float): temperature of deuterons in eV
+        T_T (float, optional): temperature of tritons in eV. If None, uses T_D for both.
+        n_samples (int): number of Monte Carlo samples (default 1e6)
+
+    Returns:
+        numpy.array: normalised DT spectrum (1/eV) evaluated at the bin centres of Ein
+    """
+    from NeSST.dress_interface import DRESS_DT_spec
+
+    if T_T is None:
+        T_T = T_D
+    return DRESS_DT_spec(T_D, T_T, n_samples, Ein)
+
+
+def QDress_DD(Ein: npt.NDArray, Tion: float, n_samples: int = int(1e6)) -> npt.NDArray:
+    """Calculates the DD primary neutron spectrum using the DRESS Monte Carlo code.
+
+    Args:
+        Ein (numpy.array): array of energy bin edges on which to compute the spectrum (eV)
+        Tion (float): temperature of the ions in eV
+        n_samples (int): number of Monte Carlo samples (default 1e6)
+
+    Returns:
+        numpy.array: normalised DD spectrum (1/eV) evaluated at the bin centres of Ein
+    """
+    from NeSST.dress_interface import DRESS_DD_spec
+
+    return DRESS_DD_spec(Tion, n_samples, Ein)
+
+
 # TT spectral shape
 def dNdE_TT(E: npt.NDArray, Tion: float) -> npt.NDArray:
     """Calculates the TT primary spectrum with Doppler broadening effect as

src/NeSST/dress_interface.py

@@ -1,22 +1,95 @@
-import dress
+import numpy.typing as npt
 
-def DRESS_DT_spec_single_T(T,n_samples,bins):
-    # Conversion to keV units of DRESS
-    T_keV = T/1e3
-    bins /= 1e3
+try:
+    import dress
+
+    dress_available = True
+except ImportError:
+    dress_available = False
+
+
+def _check_dress():
+    if not dress_available:
+        raise ImportError("pydress is required for DRESS interface functions. Install it with: pip install pydress")
+
+
+def DRESS_DT_spec(T_D: float, T_T: float, n_samples: int, bins: npt.NDArray) -> npt.NDArray:
+    """Computes DT primary neutron spectrum using DRESS with separate D and T temperatures.
+
+    Args:
+        T_D (float): temperature of deuterons in eV
+        T_T (float): temperature of tritons in eV
+        n_samples (int): number of Monte Carlo samples
+        bins (numpy.array): energy bin edges in eV
+
+    Returns:
+        numpy.array: normalised DT spectrum (1/eV) on the bin centres
+    """
+    _check_dress()
+    T_D_keV = T_D / 1e3
+    T_T_keV = T_T / 1e3
+    bins_keV = bins / 1e3
 
     reaction = dress.reactions.DTNHe4Reaction()
     spec_calc = dress.SpectrumCalculator(reaction, n_samples=n_samples)
 
+    dist_a = dress.dists.MaxwellianDistribution(T_D_keV, spec_calc.reactant_a.particle)
+    dist_b = dress.dists.MaxwellianDistribution(T_T_keV, spec_calc.reactant_b.particle)
+
+    spec_calc.reactant_a.v = dist_a.sample(spec_calc.n_samples)
+    spec_calc.reactant_b.v = dist_b.sample(spec_calc.n_samples)
+
+    spec = spec_calc(bins=bins_keV, normalize=True)
+
+    # Convert from unit-sum to unit-area (1/eV)
+    bins_w = bins[1:] - bins[:-1]
+
+    return spec / bins_w
+
+
+def DRESS_DT_spec_single_T(T: float, n_samples: int, bins: npt.NDArray) -> npt.NDArray:
+    """Computes DT primary neutron spectrum using DRESS with a single ion temperature.
+
+    Wrapper around DRESS_DT_spec for the case where D and T share the same temperature.
+
+    Args:
+        T (float): ion temperature in eV
+        n_samples (int): number of Monte Carlo samples
+        bins (numpy.array): energy bin edges in eV
+
+    Returns:
+        numpy.array: normalised DT spectrum (1/eV) on the bin centres
+    """
+    return DRESS_DT_spec(T, T, n_samples, bins)
+
+
+def DRESS_DD_spec(T: float, n_samples: int, bins: npt.NDArray) -> npt.NDArray:
+    """Computes DD primary neutron spectrum using DRESS.
+
+    Args:
+        T (float): temperature of deuterons in eV
+        n_samples (int): number of Monte Carlo samples
+        bins (numpy.array): energy bin edges in eV
+
+    Returns:
+        numpy.array: normalised DD spectrum (1/eV) on the bin centres
+    """
+    _check_dress()
+    T_keV = T / 1e3
+    bins_keV = bins / 1e3
+
+    reaction = dress.reactions.DDNHe3Reaction()
+    spec_calc = dress.SpectrumCalculator(reaction, n_samples=n_samples)
+
     dist_a = dress.dists.MaxwellianDistribution(T_keV, spec_calc.reactant_a.particle)
     dist_b = dress.dists.MaxwellianDistribution(T_keV, spec_calc.reactant_b.particle)
 
     spec_calc.reactant_a.v = dist_a.sample(spec_calc.n_samples)
     spec_calc.reactant_b.v = dist_b.sample(spec_calc.n_samples)
 
-    spec = spec_calc(bins=bins,normalize=True)
+    spec = spec_calc(bins=bins_keV, normalize=True)
 
-    # Ensure unit area for eV
-    bins_w = 1e3*(bins[1:]-bins[:-1])
+    # Convert from unit-sum to unit-area (1/eV)
+    bins_w = bins[1:] - bins[:-1]
 
-    return spec/bins_w
+    return spec / bins_w

tests/test_dress.py

@@ -0,0 +1,94 @@
+import numpy as np
+import pytest
+
+pytest.importorskip("dress", reason="pydress is required for DRESS tests")
+
+import NeSST as nst
+
+# Number of Monte Carlo samples used in tests - large enough for accurate
+# statistics but small enough to keep test runtime reasonable.
+_N_SAMPLES = int(1e5)
+
+# Tolerance for relative comparison of mean energy against Ballabio (%)
+_MEAN_RTOL = 0.01  # 1 %
+# Tolerance for relative comparison of standard deviation against Ballabio (%)
+_STDDEV_RTOL = 0.05  # 5 %
+
+
+def _spectrum_moments(Ein, spec):
+    """Return the mean and standard deviation of a normalised spectrum."""
+    Emid = 0.5 * (Ein[:-1] + Ein[1:])
+    dE = Ein[1:] - Ein[:-1]
+    total = np.sum(spec * dE)
+    mean = np.sum(Emid * spec * dE) / total
+    stddev = np.sqrt(np.sum((Emid - mean) ** 2 * spec * dE) / total)
+    return mean, stddev
+
+
+def test_QDress_DT_mean_matches_Ballabio():
+    """DRESS DT mean neutron energy should match Ballabio within 1 %."""
+    Tion = 5e3  # 5 keV in eV
+    mean_ball, stddev_ball, _ = nst.DTprimspecmoments(Tion)
+
+    Ein = np.linspace(mean_ball - 10 * stddev_ball, mean_ball + 10 * stddev_ball, 500)
+    spec = nst.QDress_DT(Ein, Tion, n_samples=_N_SAMPLES)
+
+    mean_dress, _ = _spectrum_moments(Ein, spec)
+
+    assert mean_dress == pytest.approx(mean_ball, rel=_MEAN_RTOL)
+
+
+def test_QDress_DT_stddev_matches_Ballabio():
+    """DRESS DT spectral width should match Ballabio within 5 %."""
+    Tion = 5e3  # 5 keV in eV
+    mean_ball, stddev_ball, _ = nst.DTprimspecmoments(Tion)
+
+    Ein = np.linspace(mean_ball - 10 * stddev_ball, mean_ball + 10 * stddev_ball, 500)
+    spec = nst.QDress_DT(Ein, Tion, n_samples=_N_SAMPLES)
+
+    _, stddev_dress = _spectrum_moments(Ein, spec)
+
+    assert stddev_dress == pytest.approx(stddev_ball, rel=_STDDEV_RTOL)
+
+
+def test_QDress_DT_separate_temperatures():
+    """Using different D and T temperatures should produce a different spectrum
+    from a single-temperature run."""
+    Tion = 5e3  # 5 keV in eV
+    T_D = 5e3
+    T_T = 10e3  # tritons hotter than deuterons
+
+    mean_ball, stddev_ball, _ = nst.DTprimspecmoments(Tion)
+    Ein = np.linspace(mean_ball - 10 * stddev_ball, mean_ball + 10 * stddev_ball, 500)
+
+    spec_single = nst.QDress_DT(Ein, Tion, n_samples=_N_SAMPLES)
+    spec_diff = nst.QDress_DT(Ein, T_D, T_T, n_samples=_N_SAMPLES)
+
+    # The spectra should not be identical when temperatures differ
+    assert not np.allclose(spec_single, spec_diff)
+
+
+def test_QDress_DD_mean_matches_Ballabio():
+    """DRESS DD mean neutron energy should match Ballabio within 1 %."""
+    Tion = 5e3  # 5 keV in eV
+    mean_ball, stddev_ball, _ = nst.DDprimspecmoments(Tion)
+
+    Ein = np.linspace(mean_ball - 10 * stddev_ball, mean_ball + 10 * stddev_ball, 500)
+    spec = nst.QDress_DD(Ein, Tion, n_samples=_N_SAMPLES)
+
+    mean_dress, _ = _spectrum_moments(Ein, spec)
+
+    assert mean_dress == pytest.approx(mean_ball, rel=_MEAN_RTOL)
+
+
+def test_QDress_DD_stddev_matches_Ballabio():
+    """DRESS DD spectral width should match Ballabio within 5 %."""
+    Tion = 5e3  # 5 keV in eV
+    mean_ball, stddev_ball, _ = nst.DDprimspecmoments(Tion)
+
+    Ein = np.linspace(mean_ball - 10 * stddev_ball, mean_ball + 10 * stddev_ball, 500)
+    spec = nst.QDress_DD(Ein, Tion, n_samples=_N_SAMPLES)
+
+    _, stddev_dress = _spectrum_moments(Ein, spec)
+
+    assert stddev_dress == pytest.approx(stddev_ball, rel=_STDDEV_RTOL)