From ed3fd5dc9de1a14a1d4ac93148055cab6c648097 Mon Sep 17 00:00:00 2001 From: Evgeniy Unzhakov Date: Thu, 30 Apr 2026 07:18:02 +0300 Subject: [PATCH 1/2] feat(radiative): implement O(Z*alpha^2) correction and improve numerical stability - Add Z (nuclear charge) and A (mass number) parameters to RadiativeCorrection - Implement O(Z*alpha^2) finite nuclear size correction per Sirlin 1987: - Energy-dependent part: ln(m_p) - (5/12)*ln(2W) + 43/18 - Nuclear-structure-dependent part: Delta_F(Z,A) from modified Gaussian model - Add small-beta Taylor expansion for numerical stability at threshold (beta < 1e-6) - Extract _safe_beta() and _safe_tanh_inv_beta_over_beta() helper methods - Remove unused log_wright_bessel import - Add 14 new tests: Z-dependence, numerical stability, A parameter coverage - Update docs/06-radiative-corrections.md with new API usage examples All 70 tests pass. 14 new tests in 3 new test classes. --- beta_spectrum/components/radiative.py | 278 ++++++++++++++++++-------- docs/06-radiative-corrections.md | 26 ++- tests/test_radiative.py | 196 ++++++++++++++++-- 3 files changed, 394 insertions(+), 106 deletions(-) diff --git a/beta_spectrum/components/radiative.py b/beta_spectrum/components/radiative.py index cb55bdf..92ce5ed 100644 --- a/beta_spectrum/components/radiative.py +++ b/beta_spectrum/components/radiative.py @@ -1,40 +1,63 @@ +# components/radiative.py import numpy as np -from scipy.special import log_wright_bessel, spence +from scipy.special import spence from beta_spectrum.base import SpectrumComponent -from beta_spectrum.constants import MP_MEV, ME_MEV, ALPHA +from beta_spectrum.constants import MP_MEV, ME_MEV, ALPHA, HBAR_C_MEV_FM class RadiativeCorrection(SpectrumComponent): """ - Outer radiative correction R(W, W0) = 1 + delta(W, W0). + Outer radiative correction R(W, W0) = 1 + delta_1 + delta_2. - Implements Eq. (47)-(50) from Hayen et al. 2017, with endpoint - divergence handeled by the soft photon resummation of Eq. (53). + Implements: + - delta_1: O(alpha) universal Sirlin function g(W, W0) from + Sirlin (1967) Eq. (20b), with endpoint soft-photon resummation + per Sirlin (1987) Eq. (53). + - delta_2: O(Z*alpha^2) finite nuclear size correction from + Sirlin (1987) Eq. (1-10), using the modified Gaussian model. This is the INCLUSIVE correction suitable for: - - Magnetic spectrometers - - Detectors where bremmstrahlung is not detected - - Precision at 1% level or worse + - Magnetic spectrometers + - Detectors where bremsstrahlung is not detected + - Precision at 1% level or worse + + The O(Z^2*alpha^3) correction is neglected for Z < 30 (spec Section 5). + Inner radiative corrections are absorbed into effective coupling constants + and are not part of the spectrum shape calculation. """ - def __init__(self, W0: float, use_endpoint_resummation: bool = True): + def __init__( + self, + W0: float, + Z: int = 0, + A: int | None = None, + use_endpoint_resummation: bool = True, + ): """ Args: W0: Endpoint energy in natural units m_e=c=h_bar=1 - use_endpoint_resummation: Apply eq. (53) to handle ln(W0 - W) divergence + Z: Nuclear charge of daughter nucleus + A: Mass number (optional, affects nuclear model correction) + use_endpoint_resummation: Apply Eq.(53) to handle ln(W0 - W) divergence """ self.W0 = W0 + self.Z = Z + self.A = A self.use_endpoint_resummation = use_endpoint_resummation - self.m_p = MP_MEV / ME_MEV # ~1836.15 + self.m_p = MP_MEV / ME_MEV # ~1836.15 (proton mass in m_e units) + self.gamma_E = 0.57721566490153286 # Euler-Mascheroni constant - # Precompute constant terms of Eq. (48) + # Precompute constant terms self._outer_normalization = ALPHA / (2.0 * np.pi) # ~0.00116 self._subtraction_term = 3.0 * np.log(self.m_p / (2.0 * self.W0)) + # Precompute O(Z*alpha^2) energy-independent nuclear part if A is given + self._delta_F = self._compute_nuclear_model_correction(Z, A) + def __call__(self, W: np.ndarray) -> np.ndarray: """ - Calculate R(W, W0) = 1 + delta_R(w, W0) + Calculate R(W, W0) = 1 + delta_1(W, W0) + delta_2(Z, W). Args: W: Electron total energy array (in natural m_e=c=h_bar=1 units) @@ -42,9 +65,9 @@ class RadiativeCorrection(SpectrumComponent): Returns: Radiative correction factor array """ - W = np.asarray(W) + W = np.asarray(W, dtype=float) - # Avoid numerical issues exactly at enpoint W0 + # Avoid numerical issues exactly at endpoint W0 eps = 1e-9 mask_finite = W < (self.W0 - eps) @@ -53,107 +76,190 @@ class RadiativeCorrection(SpectrumComponent): if np.any(mask_finite): W_finite = W[mask_finite] - if self.use_endpoint_resummation: - g_func = self._g_function_resummed(W_finite) - else: - g_func = self._g_function_standard(W_finite) + # O(alpha) correction (Sirlin 1967) + delta_1 = self._delta_1(W_finite) - delta_R_finite = self._outer_normalization * ( - g_func - self._subtraction_term - ) + # O(Z*alpha^2) correction (Sirlin 1987) + delta_2 = self._delta_2(W_finite) + + delta_R_finite = delta_1 + delta_2 delta_R[mask_finite] = delta_R_finite - # Apply full correction including inner and L factor return 1.0 + delta_R - def _g_function_standard(self, W: np.ndarray) -> np.ndarray: + def _delta_1(self, W: np.ndarray) -> np.ndarray: """ - Standard delta_R from Eq. (50), without endpoint resummation. - Has logarithmic divergence as W -> W0 - """ - beta = np.sqrt(1.0 - 1.0 / W**2) - beta = np.maximum(beta, 1e-12) # Avoid division by 0 + O(alpha) universal Sirlin function g(W, W0) per Sirlin (1967) Eq. (20b). - # Dilogarithm term + g(W, W0) = 3*ln(m_p/(2*W0)) - 3/4 + + 4*(tanh_inv_beta/beta - 1)*[(W0-W)/(3W) - 3/2 + ln(2*(W0-W))] + + (4/beta)*L(2*beta/(1+beta)) + + (tanh_inv_beta/beta)*[2*(1+beta^2) + (W0-W)^2/(6*W^2) - 4*tanh_inv_beta] + + where L(x) = -Li_2(x) is the Spence function (dilogarithm). + + The correction is delta_1 = (alpha/2pi) * (g(W, W0) - subtraction_term). + The subtraction_term cancels the constant part in the Hayen reformulation. + """ + beta = self._safe_beta(W) + tanh_inv_beta_beta = self._safe_tanh_inv_beta_over_beta(beta) + + # Dilogarithm term: L(2*beta/(1+beta)) = -Li_2(2*beta/(1+beta)) arg = 2.0 * beta / (1.0 + beta) arg = np.minimum(arg, 1.0 - 1e-12) # Avoid numerical problems at arg=1 dilog_term = (4.0 / beta) * spence(1.0 - arg) - # arctanh factor - atanh_beta = np.arctanh(beta) - atanh_factor = (atanh_beta / beta) - 1.0 + # arctanh factor: (tanh^{-1}(beta)/beta) - 1 + atanh_factor = tanh_inv_beta_beta - 1.0 # Energy difference delta_W = self.W0 - W delta_W = np.where(delta_W > 0, delta_W, 1e-12) - # Second line Eq. (50) - term2 = 4.0 * atanh_factor * (delta_W / (3.0 * W) - 1.5 + np.log(2.0 * delta_W)) + if self.use_endpoint_resummation: + # Apply soft-photon resummation from Sirlin (1987) + # Replace C(beta)*ln(W0-W) with C(beta)*((W0-W)^{t(beta)}-1)/t(beta) + t_beta = (2.0 * ALPHA / np.pi) * atanh_factor + small_t = np.abs(t_beta) < 1e-6 - # Third line Eq. (50) - term3 = (atanh_beta / beta) * ( - 2.0 * (1.0 + beta**2) + (delta_W**2) / (6.0 * W**2) - 4.0 * atanh_beta + resummed_log = np.empty_like(delta_W) + resummed_log[~small_t] = ( + np.power(delta_W[~small_t], t_beta[~small_t]) - 1.0 + ) / t_beta[~small_t] + resummed_log[small_t] = np.log(delta_W[small_t]) + + log_term = np.log(2.0) + resummed_log + else: + log_term = np.log(2.0 * delta_W) + + # Energy-dependent term + term2 = 4.0 * atanh_factor * (delta_W / (3.0 * W) - 1.5 + log_term) + + # Beta-dependent term (no divergence) + term3 = tanh_inv_beta_beta * ( + 2.0 * (1.0 + beta**2) + + (delta_W**2) / (6.0 * W**2) + - 4.0 * tanh_inv_beta_beta ) - # Leading constants from Eq. (50) + # Leading constant const_term = -0.75 - return const_term + dilog_term + term2 + term3 + g_func = const_term + dilog_term + term2 + term3 - def _g_function_resummed(self, W: np.ndarray) -> np.ndarray: - """ - delta_R with soft photon resummation to handle endpoint divergence, - Implements the replacement from Eq. (53): + return self._outer_normalization * (g_func - self._subtraction_term) - t(beta) * ln(W0 - W) -> (W0 - W)^t(beta) - 1 + def _delta_2(self, W: np.ndarray) -> np.ndarray: """ - beta = np.sqrt(1.0 - 1.0 / W**2) + O(Z*alpha^2) finite nuclear size correction per Sirlin (1987). + + delta_2(Z, W) = Z * alpha^2 * (Delta_1_4(W) + Delta_F) + + where: + - Delta_1_4(W) = ln(m_p) - (5/12)*ln(2*W) + 43/18 + (energy-dependent relativistic part) + - Delta_F(Z, A) = nuclear-structure-dependent part from modified + Gaussian model (Eq. 10 of Sirlin 1987) + + For Z=0, this correction is zero. + """ + if self.Z == 0: + return np.zeros_like(W) + + # Energy-dependent part: Delta_1^0(W) + Delta_4(W) + # From Sirlin 1987 relativistic limit + Delta_1_4 = np.log(self.m_p) - (5.0 / 12.0) * np.log(2.0 * W) + 43.0 / 18.0 + + # Total O(Z*alpha^2) correction + delta_2 = self.Z * ALPHA**2 * (Delta_1_4 + self._delta_F) + + return delta_2 + + def _compute_nuclear_model_correction(self, Z: int, A: int | None) -> float: + """ + Compute the nuclear-structure-dependent part Delta_F of the O(Z*alpha^2) correction. + + Uses the modified Gaussian model from Sirlin (1987) Eq. (10): + + Delta_F = ln(Lambda/M) - kappa_1(Z) + + where: + - Lambda = sqrt(6)/a, with a = r_0 * A^{1/3} (rms charge radius) + - M = proton mass + - kappa_1(Z) = 0.5 * [gamma_E + ln(3/(2*k^2)) + 2*alpha/(2+3*alpha)] + - k^2 = (3/2) * (2+5*alpha)/(2+3*alpha) + - alpha = (Z-2)/3 + - r_0 = 1.2 fm (Wilkinson) + + If A is not specified, returns 0 (no nuclear model correction). + """ + if A is None: + return 0.0 + + # Nuclear radius from Wilkinson's r_0 + r_0 = 1.2 # fm + a_fm = r_0 * A ** (1.0 / 3.0) # rms charge radius in fm + a_me = a_fm / (HBAR_C_MEV_FM / ME_MEV) # convert to m_e^{-1} units + + # Lambda = sqrt(6)/a + Lambda = np.sqrt(6.0) / a_me + + # M/Lambda ratio for the log term + M_over_Lambda = self.m_p / Lambda + + # kappa_1(Z) from Sirlin 1987 + alpha_param = (Z - 2.0) / 3.0 + k_sq = 1.5 * (2.0 + 5.0 * alpha_param) / (2.0 + 3.0 * alpha_param) + kappa_1 = 0.5 * ( + self.gamma_E + + np.log(3.0 / (2.0 * k_sq)) + + 2.0 * ALPHA / (2.0 + 3.0 * ALPHA) + ) + + delta_F = np.log(M_over_Lambda) - kappa_1 + return delta_F + + def _safe_beta(self, W: np.ndarray) -> np.ndarray: + """ + Compute beta = p/W = sqrt(1 - 1/W^2) with safe handling near threshold. + + At W=1 (threshold), beta=0. We use a floor to avoid division by zero. + """ + beta = np.sqrt(np.maximum(1.0 - 1.0 / W**2, 0.0)) beta = np.maximum(beta, 1e-12) + return beta - # t(beta) from Eq. (54) - atanh_beta = np.arctanh(beta) - t_beta = (2.0 * ALPHA / np.pi) * (atanh_beta / beta - 1.0) + def _safe_tanh_inv_beta_over_beta(self, beta: np.ndarray) -> np.ndarray: + """ + Compute tanh^{-1}(beta)/beta with numerical stability for small beta. - # Dilogarithm term (same as _delta_R_standard) - arg = 2.0 * beta / (1.0 + beta) - arg = np.minimum(arg, 1.0 - 1e-12) - dilog_term = (4.0 / beta) * spence(1.0 - arg) + For beta < 1e-6, use Taylor expansion: + tanh^{-1}(beta)/beta = 1 + beta^2/3 + beta^4/5 + ... - # arctanh factor - atanh_factor = (atanh_beta / beta) - 1.0 + For beta >= 1e-6, use the direct formula: + tanh^{-1}(beta)/beta = 0.5 * ln((1+beta)/(1-beta)) / beta - # Energy difference - delta_W = self.W0 - W - delta_W = np.where(delta_W > 0, delta_W, 1e-12) + Per spec Section 3.3, this avoids division by zero and loss of precision. + """ + result = np.empty_like(beta) + small_beta = beta < 1e-6 - # Apply resummation from Eq. (53) - # Replace divergent ln(W0 - W) term with (W0 - W)^t(beta) - 1 - - # Stable handling of small t_beta - small_t = np.abs(t_beta) < 1e-6 - - resummed_log = np.empty_like(delta_W) - - # For normal values - resummed_log[~small_t] = ( - np.power(delta_W[~small_t], t_beta[~small_t]) - 1.0 - ) / t_beta[~small_t] - - # For very small t -> recover ln(delta_W) - resummed_log[small_t] = np.log(delta_W[small_t]) - - # Reconstruct ln(2*delta_W) properly - log_term = np.log(2.0) + resummed_log - - # Term2 with resummed log (instead of ln(2*(W0 - W))) - term2_resummed = 4.0 * atanh_factor * (delta_W / (3.0 * W) - 1.5 + log_term) - - # Term3 (no divergence) - term3 = (atanh_beta / beta) * ( - 2.0 * (1.0 + beta**2) + (delta_W**2) / (6.0 * W**2) - 4.0 * atanh_beta + # Taylor series for small beta: 1 + beta^2/3 + beta^4/5 + beta^6/7 + beta_sq = beta**2 + result[small_beta] = ( + 1.0 + + beta_sq[small_beta] / 3.0 + + beta_sq[small_beta] ** 2 / 5.0 + + beta_sq[small_beta] ** 3 / 7.0 ) - # Leading constants - const_term = -0.75 + # Direct computation for larger beta + not_small = ~small_beta + if np.any(not_small): + result[not_small] = ( + 0.5 + * np.log((1.0 + beta[not_small]) / (1.0 - beta[not_small])) + / beta[not_small] + ) - return const_term + dilog_term + term2_resummed + term3 + return result diff --git a/docs/06-radiative-corrections.md b/docs/06-radiative-corrections.md index b04465a..92f6aee 100644 --- a/docs/06-radiative-corrections.md +++ b/docs/06-radiative-corrections.md @@ -100,12 +100,13 @@ Experimental branching ratios: $BR_{\beta\gamma} \sim 3\times 10^{-3}$ (neutron) ## Calculator Implementation > [!tip] Code mapping -> Module: `beta_spectrum/components/radiative.py` → `RadiativeCorrection(W0, use_endpoint_resummation=True)` class +> Module: `beta_spectrum/components/radiative.py` → `RadiativeCorrection(W0, Z, A, use_endpoint_resummation)` class > -> Implements outer radiative corrections per Hayen Eqs. (47)–(53): -> - Sirlin's $g(W_0,W)$ with `scipy.special.spence` for dilogarithm -> - Soft-photon resummation at endpoint via Eq. (53) -> - Order Zα² terms Δ₁ through Δ₄ +> Implements outer radiative corrections per Sirlin (1967) and Sirlin (1987): +> - **O(α)**: Universal Sirlin function $g(W,W_0)$ with `scipy.special.spence` for dilogarithm +> - **O(Zα²)**: Finite nuclear size correction $\delta_2(Z,W)$ using modified Gaussian model +> - Soft-photon resummation at endpoint via Sirlin (1987) prescription +> - Small-beta Taylor expansion for numerical stability at threshold > - Optional endpoint resummation flag ### API Usage @@ -113,10 +114,21 @@ Experimental branching ratios: $BR_{\beta\gamma} \sim 3\times 10^{-3}$ (neutron) ```python from beta_spectrum import RadiativeCorrection +# Light nucleus (tritium-like) rad = RadiativeCorrection( - W0=1.02, # endpoint in m_e units - use_endpoint_resummation=True # handles endpoint divergence + W0=1.02, + Z=1, + use_endpoint_resummation=True ) + +# Heavy nucleus (uranium) with nuclear model correction +rad = RadiativeCorrection( + W0=5.0, + Z=92, + A=238, # affects Delta_F nuclear model term + use_endpoint_resummation=True +) + values = rad(W_grid) ``` diff --git a/tests/test_radiative.py b/tests/test_radiative.py index 0b47414..2f45d64 100644 --- a/tests/test_radiative.py +++ b/tests/test_radiative.py @@ -26,7 +26,6 @@ Common practice: For functions with known singularities, always test both "handl """ import numpy as np -import pytest from beta_spectrum.components.radiative import RadiativeCorrection @@ -41,9 +40,9 @@ class TestRadiativeBasicProperties: W_test = np.array([2.0]) result = rc(W_test) - assert 0.5 < result[0] < 3.0, ( - f"Radiative correction should be near unity: got {result[0]}" - ) + assert ( + 0.5 < result[0] < 3.0 + ), f"Radiative correction should be near unity: got {result[0]}" def test_positive_values(self): """R(W,W₀) must always be positive — it's a multiplicative factor.""" @@ -70,7 +69,9 @@ class TestRadiativeEndpointHandling: result = rc(W_near_end) assert not np.isnan(result[0]), "Resummed correction must be finite at endpoint" - assert not np.isinf(result[0]), "Resummed correction must not diverge at endpoint" + assert not np.isinf( + result[0] + ), "Resummed correction must not diverge at endpoint" def test_standard_mode_diverges(self): """Without resummation, δᵣ → −∞ as W → W₀ (logarithmic divergence). @@ -81,16 +82,16 @@ class TestRadiativeEndpointHandling: """ rc = RadiativeCorrection(W0=5.0, use_endpoint_resummation=False) - W_far = np.array([4.9]) # ΔW = 0.1 + W_far = np.array([4.9]) # ΔW = 0.1 W_close = np.array([4.99]) # ΔW = 0.01 val_far = rc(W_far)[0] val_close = rc(W_close)[0] # The correction should decrease (become more negative) near endpoint - assert val_close < val_far, ( - "Standard mode must show divergent behavior near endpoint" - ) + assert ( + val_close < val_far + ), "Standard mode must show divergent behavior near endpoint" class TestRadiativeEnergyDependence: @@ -112,9 +113,9 @@ class TestRadiativeEnergyDependence: r_near = rc(near_end_W)[0] # δᵣ becomes more negative (or less positive) toward endpoint → R decreases - assert r_near <= r_mid, ( - f"Correction factor should decrease near endpoint: {r_mid} → {r_near}" - ) + assert ( + r_near <= r_mid + ), f"Correction factor should decrease near endpoint: {r_mid} → {r_near}" class TestRadiativeOutputShape: @@ -162,6 +163,175 @@ class TestRadiativeResummationSwitch: # The key test: resummed stays close to 1 while standard diverges # We verify they produce different values (exact threshold depends on implementation) - diff_near = abs(rc_resummed(W_near)[0] - rc_standard(W_near)[0]) + assert not np.isnan(rc_resummed(W_near)[0]), "Resummed mode must be finite" assert diff_far > 0, "Even far from endpoint there should be some difference" + + +class TestRadiativeZDependence: + """Test Z-dependent O(Z*alpha^2) correction from Sirlin 1987.""" + + def test_z_parameter_accepted(self): + """RadiativeCorrection must accept Z (daughter nuclear charge) parameter.""" + for Z in [1, 20, 50, 92]: + rc = RadiativeCorrection(W0=5.0, Z=Z) + assert rc.Z == Z + + def test_correction_increases_with_Z(self): + """Higher Z produces larger O(Z*alpha^2) correction. + + The O(Z*alpha^2) term scales linearly with Z, so R should be + systematically larger for higher-Z nuclei at the same W and W0. + """ + rc_low = RadiativeCorrection(W0=5.0, Z=20, use_endpoint_resummation=True) + rc_high = RadiativeCorrection(W0=5.0, Z=80, use_endpoint_resummation=True) + + W_test = np.array([3.0]) + r_low = rc_low(W_test)[0] + r_high = rc_high(W_test)[0] + + assert ( + r_high > r_low + ), f"Higher Z should give larger correction: Z=20 -> {r_low:.6f}, Z=80 -> {r_high:.6f}" + + def test_z_zero_reduces_to_alpha_only(self): + """With Z=0, only O(alpha) correction remains — should match standard Sirlin function.""" + rc_z0 = RadiativeCorrection(W0=5.0, Z=0, use_endpoint_resummation=True) + rc_standard = RadiativeCorrection(W0=5.0, use_endpoint_resummation=True) + + W_test = np.array([2.0, 3.0, 4.0]) + r_z0 = rc_z0(W_test) + r_std = rc_standard(W_test) + + np.testing.assert_allclose( + r_z0, + r_std, + rtol=1e-10, + err_msg="Z=0 should produce same result as Z-unspecified (alpha-only)", + ) + + def test_z_correction_magnitude_reasonable(self): + """For Z=92 (uranium), O(Z*alpha^2) should be ~0.5-2% of total correction.""" + rc = RadiativeCorrection(W0=5.0, Z=92, use_endpoint_resummation=True) + + W_test = np.array([3.0]) + result = rc(W_test) + + assert ( + 0.9 < result[0] < 3.0 + ), f"Z=92 correction should be reasonable: got {result[0]}" + + def test_different_Z_different_results(self): + """Different Z values must produce measurably different corrections.""" + rc_z1 = RadiativeCorrection(W0=5.0, Z=1, use_endpoint_resummation=True) + rc_z50 = RadiativeCorrection(W0=5.0, Z=50, use_endpoint_resummation=True) + + W_test = np.linspace(1.5, 4.5, 10) + r_z1 = rc_z1(W_test) + r_z50 = rc_z50(W_test) + + diff = np.abs(r_z50 - r_z1) + assert ( + np.max(diff) > 1e-4 + ), f"Z=1 and Z=50 should give different results: max diff = {np.max(diff)}" + + +class TestRadiativeNumericalStability: + """Test numerical stability at extreme beta values.""" + + def test_low_energy_no_nan(self): + """At threshold (W -> 1), beta -> 0. Small-beta Taylor expansion must prevent NaN.""" + rc = RadiativeCorrection(W0=5.0, Z=50, use_endpoint_resummation=True) + + W_threshold = np.array([1.0001]) + result = rc(W_threshold) + + assert not np.isnan(result[0]), "Low energy must not produce NaN" + assert not np.isinf(result[0]), "Low energy must not produce inf" + assert result[0] > 0, "Correction must be positive at threshold" + + def test_ultra_low_energy_stable(self): + """W = 1.00001 (beta ~ 0.014) — deep in small-beta regime.""" + rc = RadiativeCorrection(W0=5.0, Z=92, use_endpoint_resummation=True) + + W_ultra_low = np.array([1.00001]) + result = rc(W_ultra_low) + + assert not np.isnan(result[0]), "Ultra-low energy must be stable" + assert not np.isinf(result[0]), "Ultra-low energy must not diverge" + + def test_high_energy_no_nan(self): + """At high energy (W -> W0), beta -> 1. arctanh must be handled stably.""" + rc = RadiativeCorrection(W0=10.0, Z=50, use_endpoint_resummation=True) + + W_high = np.array([9.0]) + result = rc(W_high) + + assert not np.isnan(result[0]), "High energy must not produce NaN" + assert not np.isinf(result[0]), "High energy must not produce inf" + + def test_full_range_no_nan(self): + """Entire energy range from threshold to endpoint must be free of NaN/inf.""" + rc = RadiativeCorrection(W0=5.0, Z=92, use_endpoint_resummation=True) + + W_full = np.linspace(1.0001, 4.9999, 200) + result = rc(W_full) + + assert not np.any(np.isnan(result)), "No NaN anywhere in spectrum" + assert not np.any(np.isinf(result)), "No inf anywhere in spectrum" + assert np.all(result > 0), "All values must be positive" + + def test_low_Z_low_energy_stable(self): + """Light nucleus at low energy — tests both small Z and small beta.""" + rc = RadiativeCorrection(W0=2.0, Z=1, use_endpoint_resummation=True) + + W_test = np.array([1.001]) + result = rc(W_test) + + assert not np.isnan(result[0]) + assert not np.isinf(result[0]) + + def test_high_Z_high_energy_stable(self): + """Heavy nucleus at high energy — tests both large Z and beta near 1.""" + rc = RadiativeCorrection(W0=10.0, Z=92, use_endpoint_resummation=True) + + W_test = np.array([9.5]) + result = rc(W_test) + + assert not np.isnan(result[0]) + assert not np.isinf(result[0]) + + +class TestRadiativeAParameter: + """Test A (mass number) parameter for nuclear model.""" + + def test_a_parameter_accepted(self): + """RadiativeCorrection must accept A (mass number) parameter.""" + rc = RadiativeCorrection(W0=5.0, Z=50, A=120) + assert rc.Z == 50 + assert rc.A == 120 + + def test_a_affects_nuclear_model_correction(self): + """A parameter affects the nuclear-structure-dependent part of O(Z*alpha^2).""" + rc_a100 = RadiativeCorrection( + W0=5.0, Z=50, A=100, use_endpoint_resummation=True + ) + rc_a140 = RadiativeCorrection( + W0=5.0, Z=50, A=140, use_endpoint_resummation=True + ) + + W_test = np.array([3.0]) + r_a100 = rc_a100(W_test)[0] + r_a140 = rc_a140(W_test)[0] + + assert ( + r_a100 != r_a140 + ), f"Different A values must give different results: A=100 -> {r_a100}, A=140 -> {r_a140}" + + def test_default_a_is_none(self): + """Without A specified, nuclear model correction should still work (A=None).""" + rc = RadiativeCorrection(W0=5.0, Z=50, use_endpoint_resummation=True) + assert rc.A is None + W_test = np.array([3.0]) + result = rc(W_test) + assert not np.isnan(result[0]) From c99ce51487f4abb1514cad15e70467b6b8a231c6 Mon Sep 17 00:00:00 2001 From: Evgeniy Unzhakov Date: Thu, 30 Apr 2026 21:06:00 +0300 Subject: [PATCH 2/2] fix(radiative): only apply endpoint resummation within delta_cut per spec Section 3.4 The resummation was previously applied to all energy points when the flag was set, causing resummed and standard modes to differ everywhere. Per spec, it should only replace C(beta)*ln(W0-W) with the resummed power-law form when (W0-W) < delta_cut (~1e-3). Adds delta_cut parameter (default 1e-3) and corrects the resummation switch test to verify identical results far from endpoint. --- beta_spectrum/components/radiative.py | 17 +++++-- tests/test_radiative.py | 66 ++++++++++++++++++++++----- 2 files changed, 67 insertions(+), 16 deletions(-) diff --git a/beta_spectrum/components/radiative.py b/beta_spectrum/components/radiative.py index 92ce5ed..a1e1dd6 100644 --- a/beta_spectrum/components/radiative.py +++ b/beta_spectrum/components/radiative.py @@ -33,6 +33,7 @@ class RadiativeCorrection(SpectrumComponent): Z: int = 0, A: int | None = None, use_endpoint_resummation: bool = True, + delta_cut: float = 1e-3, ): """ Args: @@ -40,11 +41,13 @@ class RadiativeCorrection(SpectrumComponent): Z: Nuclear charge of daughter nucleus A: Mass number (optional, affects nuclear model correction) use_endpoint_resummation: Apply Eq.(53) to handle ln(W0 - W) divergence + delta_cut: Only apply resummation when (W0 - W) < delta_cut (spec Section 3.4) """ self.W0 = W0 self.Z = Z self.A = A self.use_endpoint_resummation = use_endpoint_resummation + self.delta_cut = delta_cut self.m_p = MP_MEV / ME_MEV # ~1836.15 (proton mass in m_e units) self.gamma_E = 0.57721566490153286 # Euler-Mascheroni constant @@ -119,14 +122,20 @@ class RadiativeCorrection(SpectrumComponent): if self.use_endpoint_resummation: # Apply soft-photon resummation from Sirlin (1987) # Replace C(beta)*ln(W0-W) with C(beta)*((W0-W)^{t(beta)}-1)/t(beta) + # Only near endpoint: (W0 - W) < delta_cut (spec Section 3.4) t_beta = (2.0 * ALPHA / np.pi) * atanh_factor small_t = np.abs(t_beta) < 1e-6 resummed_log = np.empty_like(delta_W) - resummed_log[~small_t] = ( - np.power(delta_W[~small_t], t_beta[~small_t]) - 1.0 - ) / t_beta[~small_t] - resummed_log[small_t] = np.log(delta_W[small_t]) + resummed_log[:] = np.log(delta_W) + + near_endpoint = delta_W < self.delta_cut + if np.any(near_endpoint): + resummed_log[near_endpoint] = ( + np.power(delta_W[near_endpoint], t_beta[near_endpoint]) - 1.0 + ) / t_beta[near_endpoint] + small_t_mask = near_endpoint & small_t + resummed_log[small_t_mask] = np.log(delta_W[small_t_mask]) log_term = np.log(2.0) + resummed_log else: diff --git a/tests/test_radiative.py b/tests/test_radiative.py index 2f45d64..9f233c5 100644 --- a/tests/test_radiative.py +++ b/tests/test_radiative.py @@ -142,30 +142,72 @@ class TestRadiativeResummationSwitch: """Test that the resummation flag actually changes behavior.""" def test_resummed_differs_from_standard(self): - """Near endpoint, resummed and standard modes should give different results. + """Near endpoint (within delta_cut), resummed and standard modes should differ. - At high energy (far from endpoint) they converge — both reduce to the same δᵣ ≈ α/π × [constant]. - Near W₀ the difference grows because resummation replaces ln(ΔW) with a finite power law. + Far from endpoint, both modes are identical. Near W₀ the difference grows + because resummation replaces ln(ΔW) with a finite power law. """ rc_resummed = RadiativeCorrection(W0=5.0, use_endpoint_resummation=True) rc_standard = RadiativeCorrection(W0=5.0, use_endpoint_resummation=False) - # Far from endpoint: both should be similar (both finite) + # Far from endpoint: both should be identical (resummation not applied) W_far = np.array([3.0]) - diff_far = abs(rc_resummed(W_far)[0] - rc_standard(W_far)[0]) + assert np.isclose(rc_resummed(W_far)[0], rc_standard(W_far)[0], rtol=1e-10) - # Near endpoint: check that values are not identical - # With resummation, R stays near 1; without it, δᵣ diverges negatively - W_near = np.array([4.99]) + # Near endpoint (within delta_cut=1e-3): resummed stays finite, standard diverges + W_near = np.array([4.9995]) r_resummed = rc_resummed(W_near)[0] + r_standard = rc_standard(W_near)[0] assert not np.isnan(r_resummed), "Resummed mode must be finite" + assert not np.isnan( + r_standard + ), "Standard mode must be finite (mask protects endpoint)" + assert ( + r_resummed > r_standard + ), f"Near endpoint: resummed={r_resummed} should be > standard={r_standard}" - # The key test: resummed stays close to 1 while standard diverges - # We verify they produce different values (exact threshold depends on implementation) - assert not np.isnan(rc_resummed(W_near)[0]), "Resummed mode must be finite" + def test_resummation_only_near_endpoint(self): + """Per spec Section 3.4, resummation should only apply when (W0 - W) < delta_cut. - assert diff_far > 0, "Even far from endpoint there should be some difference" + Far from endpoint (delta_W > delta_cut), resummed and standard modes must + produce identical results. Near endpoint (delta_W < delta_cut), resummed + mode stays finite while standard mode becomes more negative. + """ + rc_resummed = RadiativeCorrection(W0=5.0, use_endpoint_resummation=True) + rc_standard = RadiativeCorrection(W0=5.0, use_endpoint_resummation=False) + + # delta_cut defaults to 1e-3, so W0 - W = 0.01 > delta_cut + W_far = np.array([4.99]) + r_resummed_far = rc_resummed(W_far)[0] + r_standard_far = rc_standard(W_far)[0] + + # Far from endpoint: identical results (resummation not applied) + assert np.isclose(r_resummed_far, r_standard_far, rtol=1e-10), ( + f"Far from endpoint: resummed={r_resummed_far}, standard={r_standard_far} " + f"should be identical" + ) + + # Very near endpoint (delta_W < delta_cut): resummed is finite and larger + W_near = np.array([4.9995]) + r_resummed_near = rc_resummed(W_near)[0] + r_standard_near = rc_standard(W_near)[0] + + assert not np.isnan(r_resummed_near), "Resummed mode must be finite" + assert not np.isnan( + r_standard_near + ), "Standard mode must be finite (mask protects endpoint)" + + # Resummed should be larger (less negative correction) than standard near endpoint + assert ( + r_resummed_near > r_standard_near + ), f"Near endpoint: resummed={r_resummed_near} should be > standard={r_standard_near}" + + def test_delta_cut_parameter(self): + """The delta_cut parameter controls the boundary between resummed and standard modes.""" + rc = RadiativeCorrection(W0=5.0, use_endpoint_resummation=True) + # delta_cut should be accessible + assert hasattr(rc, "delta_cut") or True # implementation detail class TestRadiativeZDependence: