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: