unzhakov/beta_shape_pnpi
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

test: complete migration to structured test suite

Browse Source 
Phase 1: Removed 6 duplicate old test files (fermi, finite_size, screening,
exchange, radiative, phase_space) — their tests are now in physics/ and quality/
using shared helpers from common/.

Phase 2: Migrated remaining 7 old test files:
- test_spectrum.py → integration/test_spectrum_pipeline.py (integration tests)
- test_fitter.py → physics/test_fitter.py (fit quality tests)
- test_cw_extractor.py → physics/test_cw_extractor.py (extraction accuracy)
- test_multi_branch.py → physics/test_multi_branch.py (branching physics)
- test_detector_response.py → physics/ + quality/ (split into physics + API)
- test_nuclear_data.py → quality/test_nuclear_data.py (parsing/API)
- test_logging_utils.py → quality/test_logging_utils.py (API tests)

Removed 12 old flat files, created 4 new files.
Total: 251 tests passing (6 skipped), down from 230 → 251 (net +21
from better coverage, not duplicates).

Test structure:
- tests/common/ — shared property/API helpers (eliminates ~18 duplicates)
- tests/physics/ — 10 files, unique physics behavior tests
- tests/quality/ — 8 files, API contract + parsing tests
- tests/integration/ — 1 file, end-to-end pipeline tests
This commit is contained in:
unzhakov
2026-05-15 15:52:50 +03:00
parent fea151d8d7
commit c78d0655a0
15 changed files with 145 additions and 1325 deletions
Download Patch File Download Diff File Expand all files Collapse all files
tests/physics/test_cw_extractor.py
+675
View File
@@ -0,0 +1,675 @@
# tests/test_cw_extractor.py — Tests for C(W) shape factor extraction
"""
Why test this?
--------------
The C(W) extractor is the core scientific tool for extracting the shape factor
from measured beta spectra. It must correctly handle:
1. Energy unit conversions (keV ↔ m_e units)
2. Theoretical factor computation (Fermi, screening, radiative, etc.)
3. Normalization of extracted C(W)
4. Kurie plot endpoint determination
5. g_V, g_A extraction from C(W) data
We test with synthetic data where the true answer is known, so we can verify
the extraction recovers the correct values.
"""
import numpy as np
import pytest
from beta_spectrum import (
BetaSpectrum,
CWExtractor,
CWExtractionResult,
GVAExtractionResult,
SpectrumConfig,
T_to_W,
W_to_T,
)
class TestCWExtractionBasic:
"""Basic C(W) extraction tests."""
def _make_synthetic_data(
self,
endpoint_keV: float = 294.0,
n_points: int = 100,
seed: int = 42,
):
"""Create synthetic beta spectrum with C(W) = 1 (allowed decay)."""
np.random.seed(seed)
# Energy grid
energies_keV = np.linspace(5, endpoint_keV * 0.98, n_points)
# Build theoretical spectrum with C(W) = 1
config = SpectrumConfig(
Z_parent=43,
Z_daughter=44,
A_number=99,
endpoint_MeV=endpoint_keV / 1000.0,
transition_type="F1",
use_screening=True,
use_exchange=True,
use_finite_size=True,
use_radiative=True,
)
spectrum = BetaSpectrum.from_config(config)
W, E_MeV = spectrum.get_energy_grid(config)
# Interpolate to our energy grid
from scipy.interpolate import interp1d
flux_interp = interp1d(
E_MeV * 1000, spectrum(W), kind="linear", fill_value=0.0, bounds_error=False
)
flux_theory = flux_interp(energies_keV)
# Add noise (Poisson-like)
noise_std = np.sqrt(np.maximum(flux_theory, 1e-10))
flux = flux_theory + np.random.normal(0, noise_std, len(flux_theory))
flux = np.maximum(flux, 0.0) # Flux cannot be negative
flux_errors = np.sqrt(np.maximum(flux, 1e-10))
return energies_keV, flux, flux_errors, config
def test_CW_extraction_allowed_decay(self):
"""For allowed decay (C(W)=1), extracted C(W) should be ~1."""
energies_keV, flux, flux_errors, config = self._make_synthetic_data()
extractor = CWExtractor(config, flux, energies_keV, flux_errors)
result = extractor.extract_CW(endpoint_keV=294.0)
# C(W) should be close to 1 in mid-energy region
mid_mask = (energies_keV > 50) & (energies_keV < 200)
if np.any(mid_mask):
avg_cw = np.mean(result.cw_values[mid_mask])
assert (
abs(avg_cw - 1.0) < 0.3
), f"Allowed decay C(W) should be ~1, got {avg_cw}"
def test_CW_extraction_preserves_endpoint(self):
"""Extracted C(W) should report correct endpoint."""
energies_keV, flux, flux_errors, config = self._make_synthetic_data()
extractor = CWExtractor(config, flux, energies_keV, flux_errors)
result = extractor.extract_CW(endpoint_keV=294.0)
assert (
abs(result.endpoint_keV - 294.0) < 1.0
), f"Endpoint should be ~294 keV, got {result.endpoint_keV}"
def test_CW_extraction_errors_positive(self):
"""C(W) uncertainties should be positive."""
energies_keV, flux, flux_errors, config = self._make_synthetic_data()
extractor = CWExtractor(config, flux, energies_keV, flux_errors)
result = extractor.extract_CW(endpoint_keV=294.0)
assert np.all(result.cw_errors > 0), "All errors should be positive"
def test_CW_extraction_clipped(self):
"""C(W) values should be clipped to [0.01, 10.0]."""
energies_keV, flux, flux_errors, config = self._make_synthetic_data()
extractor = CWExtractor(config, flux, energies_keV, flux_errors)
result = extractor.extract_CW(endpoint_keV=294.0)
assert np.all(result.cw_values >= 0.01), "C(W) should be >= 0.01"
assert np.all(result.cw_values <= 10.0), "C(W) should be <= 10.0"
class TestCWParametrization:
"""Test C(W) parametrization fitting."""
def _make_CW_data(
self,
endpoint_keV: float = 294.0,
n_points: int = 60,
seed: int = 42,
):
"""Create synthetic C(W) data with known parametrization."""
np.random.seed(seed)
config = SpectrumConfig(
Z_parent=43,
Z_daughter=44,
A_number=99,
endpoint_MeV=endpoint_keV / 1000.0,
transition_type="F1",
use_screening=True,
use_exchange=True,
use_finite_size=True,
use_radiative=True,
)
# Generate C(W) values with known form
W = np.linspace(T_to_W(0.05), T_to_W(endpoint_keV * 0.95), n_points)
E_keV = np.array([W_to_T(w) for w in W])
# Linear C(W): C(W) = 1.0 + 0.5*(W-1)
true_a0, true_a1 = 1.0, 0.5
cw_true = true_a0 + true_a1 * (W - 1.0)
cw_noise = cw_true + np.random.normal(0, 0.05, len(W))
# C(W) = flux / (C0 * phase_space * Fermi * corrections * (W0-W)^2)
# So: flux = C(W) * C0 * phase_space * Fermi * corrections * (W0-W)^2
p = np.sqrt(np.maximum(W**2 - 1.0, 0.0))
from beta_spectrum.components.fermi import FermiFunction
F = FermiFunction(Z=config.Z_parent, A=config.A_number)(W)
C0 = 1.0 # Normalization
endpoint_factor = (T_to_W(endpoint_keV / 1000.0) - W) ** 2
flux_synthetic = cw_noise * C0 * p * W * F * endpoint_factor
flux_synthetic = np.maximum(flux_synthetic, 0)
flux_errors = np.sqrt(np.maximum(flux_synthetic, 1.0))
return E_keV, flux_synthetic, flux_errors, config, true_a0, true_a1
def test_constant_parametrization(self):
"""Fit C(W) = a₀ with known constant value."""
np.random.seed(42)
energies_keV = np.linspace(10, 250, 50)
W = np.array([T_to_W(e / 1000.0) for e in energies_keV])
cw_true = np.full_like(W, 1.0)
cw = cw_true + np.random.normal(0, 0.05, len(W))
config = SpectrumConfig(
Z_parent=43,
Z_daughter=44,
A_number=99,
endpoint_MeV=0.294,
transition_type="F1",
use_screening=True,
use_exchange=True,
use_finite_size=True,
use_radiative=True,
)
extractor = CWExtractor(
config, cw, energies_keV, flux_errors=np.ones_like(cw) * 0.05
)
cw_result = CWExtractionResult(
energies_W=W.copy(),
energies_keV=energies_keV.copy(),
cw_values=cw.copy(),
cw_errors=np.ones_like(cw) * 0.05,
flux=cw.copy(),
flux_errors=np.ones_like(cw) * 0.05,
endpoint_W=T_to_W(0.294),
endpoint_keV=294.0,
)
result = extractor.fit_CW_parametrization(cw_result, parametrization="constant")
assert (
abs(result.parameters[0] - 1.0) < 0.1
), f"Constant C(W) should be ~1, got {result.parameters[0]}"
def test_linear_parametrization(self):
"""Fit C(W) = a₀ + a₁·(W-1) with known coefficients."""
np.random.seed(42)
energies_keV = np.linspace(10, 250, 50)
W = np.array([T_to_W(e / 1000.0) for e in energies_keV])
true_a0, true_a1 = 1.0, 0.5
cw = true_a0 + true_a1 * (W - 1.0) + np.random.normal(0, 0.05, len(W))
config = SpectrumConfig(
Z_parent=43,
Z_daughter=44,
A_number=99,
endpoint_MeV=0.294,
transition_type="F1",
use_screening=True,
use_exchange=True,
use_finite_size=True,
use_radiative=True,
)
extractor = CWExtractor(
config, cw, energies_keV, flux_errors=np.ones_like(cw) * 0.05
)
cw_result = CWExtractionResult(
energies_W=W.copy(),
energies_keV=energies_keV.copy(),
cw_values=cw.copy(),
cw_errors=np.ones_like(cw) * 0.05,
flux=cw.copy(),
flux_errors=np.ones_like(cw) * 0.05,
endpoint_W=T_to_W(0.294),
endpoint_keV=294.0,
)
result = extractor.fit_CW_parametrization(cw_result, parametrization="linear")
assert (
abs(result.parameters[0] - true_a0) < 0.15
), f"a₀: expected {true_a0}, got {result.parameters[0]}"
assert (
abs(result.parameters[1] - true_a1) < 0.5
), f"a₁: expected {true_a1}, got {result.parameters[1]}"
def test_quadratic_parametrization(self):
"""Fit C(W) = a₀ + a₁·(W-1) + a₂·(W-1)² with known coefficients."""
np.random.seed(42)
energies_keV = np.linspace(10, 250, 60)
W = np.array([T_to_W(e / 1000.0) for e in energies_keV])
true_a0, true_a1, true_a2 = 1.0, 0.3, -0.1
cw = true_a0 + true_a1 * (W - 1.0) + true_a2 * (W - 1.0) ** 2
cw += np.random.normal(0, 0.03, len(W))
config = SpectrumConfig(
Z_parent=43,
Z_daughter=44,
A_number=99,
endpoint_MeV=0.294,
transition_type="F1",
use_screening=True,
use_exchange=True,
use_finite_size=True,
use_radiative=True,
)
extractor = CWExtractor(
config, cw, energies_keV, flux_errors=np.ones_like(cw) * 0.03
)
cw_result = CWExtractionResult(
energies_W=W.copy(),
energies_keV=energies_keV.copy(),
cw_values=cw.copy(),
cw_errors=np.ones_like(cw) * 0.03,
flux=cw.copy(),
flux_errors=np.ones_like(cw) * 0.03,
endpoint_W=T_to_W(0.294),
endpoint_keV=294.0,
)
result = extractor.fit_CW_parametrization(
cw_result, parametrization="quadratic"
)
assert abs(result.parameters[0] - true_a0) < 0.15
assert abs(result.parameters[1] - true_a1) < 0.5
assert abs(result.parameters[2] - true_a2) < 0.5
def test_invalid_parametrization(self):
"""Invalid parametrization should raise ValueError."""
np.random.seed(42)
energies_keV = np.linspace(10, 250, 30)
W = np.array([T_to_W(e / 1000.0) for e in energies_keV])
cw = np.ones_like(W)
config = SpectrumConfig(
Z_parent=43,
Z_daughter=44,
A_number=99,
endpoint_MeV=0.294,
)
extractor = CWExtractor(config, cw, energies_keV)
cw_result = CWExtractionResult(
energies_W=W.copy(),
energies_keV=energies_keV.copy(),
cw_values=cw.copy(),
cw_errors=np.ones_like(cw) * 0.05,
flux=cw.copy(),
flux_errors=np.ones_like(cw) * 0.05,
endpoint_W=T_to_W(0.294),
endpoint_keV=294.0,
)
with pytest.raises(ValueError, match="Unknown parametrization"):
extractor.fit_CW_parametrization(cw_result, parametrization="invalid")
class TestGVAExtraction:
"""Test g_V, g_A extraction from C(W) data."""
def test_gV_gA_pure_fermi(self):
"""For pure Fermi transition, g_A ≈ 0."""
np.random.seed(42)
energies_keV = np.linspace(10, 250, 50)
W = np.array([T_to_W(e / 1000.0) for e in energies_keV])
# Pure Fermi: C(W) = g_V² (constant)
true_gV = 0.5
cw = true_gV**2 + np.random.normal(0, 0.02, len(W))
config = SpectrumConfig(
Z_parent=43,
Z_daughter=44,
A_number=99,
endpoint_MeV=0.294,
transition_type="F1",
use_screening=True,
use_exchange=True,
use_finite_size=True,
use_radiative=True,
)
extractor = CWExtractor(
config, cw, energies_keV, flux_errors=np.ones_like(cw) * 0.02
)
cw_result = CWExtractionResult(
energies_W=W.copy(),
energies_keV=energies_keV.copy(),
cw_values=cw.copy(),
cw_errors=np.ones_like(cw) * 0.02,
flux=cw.copy(),
flux_errors=np.ones_like(cw) * 0.02,
endpoint_W=T_to_W(0.294),
endpoint_keV=294.0,
)
result = extractor.fit_gV_gA(cw_result, M_F=1.0, M_GT=0.0)
assert isinstance(result, GVAExtractionResult)
assert (
abs(result.g_V - true_gV) < 0.15
), f"g_V: expected {true_gV}, got {result.g_V}"
def test_gV_gA_pure_gt(self):
"""For pure Gamow-Teller, g_V ≈ 0."""
np.random.seed(42)
energies_keV = np.linspace(10, 250, 50)
W = np.array([T_to_W(e / 1000.0) for e in energies_keV])
# Pure GT: C(W) = g_A² (constant)
true_gA = 0.8
cw = true_gA**2 + np.random.normal(0, 0.02, len(W))
config = SpectrumConfig(
Z_parent=43,
Z_daughter=44,
A_number=99,
endpoint_MeV=0.294,
transition_type="F1",
use_screening=True,
use_exchange=True,
use_finite_size=True,
use_radiative=True,
)
extractor = CWExtractor(
config, cw, energies_keV, flux_errors=np.ones_like(cw) * 0.02
)
cw_result = CWExtractionResult(
energies_W=W.copy(),
energies_keV=energies_keV.copy(),
cw_values=cw.copy(),
cw_errors=np.ones_like(cw) * 0.02,
flux=cw.copy(),
flux_errors=np.ones_like(cw) * 0.02,
endpoint_W=T_to_W(0.294),
endpoint_keV=294.0,
)
result = extractor.fit_gV_gA(cw_result, M_F=0.0, M_GT=1.0)
assert isinstance(result, GVAExtractionResult)
assert (
abs(result.g_A - true_gA) < 0.15
), f"g_A: expected {true_gA}, got {result.g_A}"
def test_gV_gA_mixed_transition(self):
"""For mixed transition, extract both g_V and g_A."""
np.random.seed(42)
energies_keV = np.linspace(10, 250, 50)
W = np.array([T_to_W(e / 1000.0) for e in energies_keV])
# Mixed: C(W) = g_V² + g_A² (constant for simplicity)
true_gV, true_gA = 0.4, 0.6
cw = true_gV**2 + true_gA**2 + np.random.normal(0, 0.02, len(W))
config = SpectrumConfig(
Z_parent=43,
Z_daughter=44,
A_number=99,
endpoint_MeV=0.294,
transition_type="F1",
use_screening=True,
use_exchange=True,
use_finite_size=True,
use_radiative=True,
)
extractor = CWExtractor(
config, cw, energies_keV, flux_errors=np.ones_like(cw) * 0.02
)
cw_result = CWExtractionResult(
energies_W=W.copy(),
energies_keV=energies_keV.copy(),
cw_values=cw.copy(),
cw_errors=np.ones_like(cw) * 0.02,
flux=cw.copy(),
flux_errors=np.ones_like(cw) * 0.02,
endpoint_W=T_to_W(0.294),
endpoint_keV=294.0,
)
result = extractor.fit_gV_gA(cw_result, M_F=1.0, M_GT=1.0)
assert isinstance(result, GVAExtractionResult)
# Check that sum of squares is correct
reconstructed = result.g_V**2 + result.g_A**2
expected = true_gV**2 + true_gA**2
assert (
abs(reconstructed - expected) < 0.15
), f"g_V²+g_A²: expected {expected}, got {reconstructed}"
def test_gV_gA_summary(self):
"""summary() should return a non-empty string with key info."""
np.random.seed(42)
energies_keV = np.linspace(10, 250, 50)
W = np.array([T_to_W(e / 1000.0) for e in energies_keV])
cw = np.full_like(W, 0.5)
config = SpectrumConfig(
Z_parent=43,
Z_daughter=44,
A_number=99,
endpoint_MeV=0.294,
)
extractor = CWExtractor(config, cw, energies_keV)
cw_result = CWExtractionResult(
energies_W=W.copy(),
energies_keV=energies_keV.copy(),
cw_values=cw.copy(),
cw_errors=np.ones_like(cw) * 0.02,
flux=cw.copy(),
flux_errors=np.ones_like(cw) * 0.02,
endpoint_W=T_to_W(0.294),
endpoint_keV=294.0,
)
result = extractor.fit_gV_gA(cw_result, M_F=1.0, M_GT=0.0)
summary = result.summary()
assert "g_V" in summary
assert "g_A" in summary
assert "χ²" in summary
assert "p-value" in summary
class TestKurieAnalysis:
"""Test Kurie plot analysis."""
def test_kurie_analysis_allowed_decay(self):
"""Kurie analysis on allowed decay should produce reasonable C(W)."""
np.random.seed(42)
energies_keV = np.linspace(5, 280, 80)
config = SpectrumConfig(
Z_parent=43,
Z_daughter=44,
A_number=99,
endpoint_MeV=0.294,
transition_type="F1",
use_screening=True,
use_exchange=True,
use_finite_size=True,
use_radiative=True,
)
spectrum = BetaSpectrum.from_config(config)
W, E_MeV = spectrum.get_energy_grid(config)
flux_theory = spectrum(W)
# Interpolate to our energy grid
from scipy.interpolate import interp1d
flux_interp = interp1d(
E_MeV * 1000, flux_theory, kind="linear", fill_value=0.0, bounds_error=False
)
flux = flux_interp(energies_keV)
# Add noise
flux += np.random.normal(0, 0.02 * np.max(flux), len(flux))
extractor = CWExtractor(config, flux, energies_keV)
cw_result, fit_result = extractor.kurie_analysis(endpoint_keV=294.0)
assert isinstance(cw_result, CWExtractionResult)
assert cw_result.fit_result is not None
assert np.all(cw_result.cw_values > 0)
def test_kurie_plot_values(self):
"""Kurie plot should produce valid values."""
np.random.seed(42)
energies_keV = np.linspace(10, 250, 50)
W = np.array([T_to_W(e / 1000.0) for e in energies_keV])
cw = np.ones_like(W) * 1.0
cw_result = CWExtractionResult(
energies_W=W.copy(),
energies_keV=energies_keV.copy(),
cw_values=cw.copy(),
cw_errors=np.ones_like(cw) * 0.05,
flux=cw.copy(),
flux_errors=np.ones_like(cw) * 0.05,
endpoint_W=T_to_W(0.294),
endpoint_keV=294.0,
)
E_plot, kurie = cw_result.kurie_plot()
assert len(E_plot) == len(energies_keV)
assert len(kurie) == len(energies_keV)
assert np.all(kurie >= 0), "Kurie values should be non-negative"
class TestEdgeCases:
"""Test edge cases and error handling."""
def test_not_enough_points_for_fit(self):
"""Fitting with too few points should raise ValueError."""
np.random.seed(42)
energies_keV = np.array([10.0, 20.0, 30.0])
W = np.array([T_to_W(e / 1000.0) for e in energies_keV])
cw = np.ones_like(W)
config = SpectrumConfig(
Z_parent=43,
Z_daughter=44,
A_number=99,
endpoint_MeV=0.294,
)
extractor = CWExtractor(config, cw, energies_keV)
cw_result = CWExtractionResult(
energies_W=W.copy(),
energies_keV=energies_keV.copy(),
cw_values=cw.copy(),
cw_errors=np.ones_like(cw) * 0.05,
flux=cw.copy(),
flux_errors=np.ones_like(cw) * 0.05,
endpoint_W=T_to_W(0.294),
endpoint_keV=294.0,
)
with pytest.raises(ValueError, match="Not enough data points"):
extractor.fit_CW_parametrization(cw_result, parametrization="quadratic")
def test_not_enough_points_for_gV_gA(self):
"""g_V/g_A fit with too few points should raise ValueError."""
np.random.seed(42)
energies_keV = np.array([10.0, 20.0, 30.0])
W = np.array([T_to_W(e / 1000.0) for e in energies_keV])
cw = np.ones_like(W)
config = SpectrumConfig(
Z_parent=43,
Z_daughter=44,
A_number=99,
endpoint_MeV=0.294,
)
extractor = CWExtractor(config, cw, energies_keV)
cw_result = CWExtractionResult(
energies_W=W.copy(),
energies_keV=energies_keV.copy(),
cw_values=cw.copy(),
cw_errors=np.ones_like(cw) * 0.05,
flux=cw.copy(),
flux_errors=np.ones_like(cw) * 0.05,
endpoint_W=T_to_W(0.294),
endpoint_keV=294.0,
)
with pytest.raises(ValueError, match="Not enough data points"):
extractor.fit_gV_gA(cw_result, M_F=1.0, M_GT=1.0)
def test_default_uncertainties(self):
"""Default uncertainties should be sqrt(flux)."""
np.random.seed(42)
energies_keV = np.linspace(10, 250, 30)
flux = np.array(
[
100.0,
200.0,
150.0,
180.0,
120.0,
90.0,
80.0,
70.0,
60.0,
50.0,
45.0,
40.0,
35.0,
30.0,
25.0,
20.0,
18.0,
15.0,
12.0,
10.0,
8.0,
7.0,
6.0,
5.0,
4.0,
3.0,
2.0,
1.5,
1.0,
0.5,
]
)
config = SpectrumConfig(
Z_parent=43,
Z_daughter=44,
A_number=99,
endpoint_MeV=0.294,
)
extractor = CWExtractor(config, flux, energies_keV)
expected_unc = np.sqrt(np.maximum(flux, 1.0))
np.testing.assert_allclose(extractor.flux_errors, expected_unc)
tests/physics/test_detector_response.py
+328
View File
@@ -0,0 +1,328 @@
"""
Physics tests for DetectorResponse — unique physical behavior.
Tests detector response shapes, resolution, and convolution physics.
"""
import numpy as np
import pytest
from beta_spectrum.components.detector_response import DetectorResponse
class TestAnalyticalGaussianResponse:
"""Test pure Gaussian detector response."""
def test_gaussian_peak_position(self):
"""Gaussian response should peak at the input energy."""
channel_energies = np.linspace(0.0, 5.0, 4096)
resp = DetectorResponse(channel_energies, model="gaussian", sigma_a=0.005)
energy = 2.0
response = resp.get_response(energy)
peak_idx = np.argmax(response)
assert np.isclose(
channel_energies[peak_idx], energy, atol=0.01
), "Gaussian peak should be at input energy"
def test_gaussian_normalization(self):
"""Gaussian response should integrate to 1 (unit area)."""
channel_energies = np.linspace(0.0, 5.0, 8192)
resp = DetectorResponse(channel_energies, model="gaussian", sigma_a=0.005)
energy = 2.0
response = resp.get_response(energy)
integral = np.trapezoid(response, channel_energies)
assert np.isclose(
integral, 1.0, rtol=1e-2
), f"Gaussian should normalize to 1, got {integral}"
def test_gaussian_sigma_dependence(self):
"""Larger sigma_a should produce wider Gaussian."""
channel_energies = np.linspace(0.0, 5.0, 4096)
resp_narrow = DetectorResponse(
channel_energies, model="gaussian", sigma_a=0.002
)
resp_wide = DetectorResponse(channel_energies, model="gaussian", sigma_a=0.02)
energy = 2.0
resp_n = resp_narrow.get_response(energy)
resp_w = resp_wide.get_response(energy)
half_max_n = np.max(resp_n) / 2
half_max_w = np.max(resp_w) / 2
fwhm_n = np.sum(resp_n > half_max_n)
fwhm_w = np.sum(resp_w > half_max_w)
assert fwhm_w > fwhm_n, "Larger sigma should produce wider peak"
def test_delta_limit(self):
"""Zero sigma should produce delta function at peak channel."""
channel_energies = np.linspace(0.0, 5.0, 4096)
resp = DetectorResponse(
channel_energies, model="gaussian", sigma_a=0.0, fano_factor=0.0
)
energy = 2.0
response = resp.get_response(energy)
assert (
np.sum(response > 0) == 1
), "Delta function should have single non-zero element"
assert np.isclose(response[np.argmax(response)], 1.0)
class TestAnalyticalGaussianTailResponse:
"""Test Gaussian + exponential tail detector response."""
def test_tail_increases_low_energy_counts(self):
"""Adding tail should produce non-zero tail contribution below peak."""
channel_energies = np.linspace(0.0, 5.0, 4096)
resp_with_tail = DetectorResponse(
channel_energies,
model="gaussian_tail",
sigma_a=0.005,
tail_fraction=0.3,
tau=0.01,
)
energy = 2.0
response = resp_with_tail.get_response(energy)
below_peak = channel_energies < energy
assert (
np.sum(response[below_peak]) > 0
), "Response below peak should be non-zero with tail enabled"
peak_idx = np.argmax(response)
assert peak_idx > 0, "Peak should not be at first channel"
assert (
response[peak_idx] >= response[peak_idx - 1]
), "Response at peak should be >= response just below peak"
def test_tail_fraction_effect(self):
"""Higher tail_fraction should produce more tail contribution."""
channel_energies = np.linspace(0.0, 5.0, 4096)
resp_low = DetectorResponse(
channel_energies,
model="gaussian_tail",
sigma_a=0.005,
tail_fraction=0.1,
tau=0.01,
)
resp_high = DetectorResponse(
channel_energies,
model="gaussian_tail",
sigma_a=0.005,
tail_fraction=0.5,
tau=0.01,
)
resp_gaussian_only = DetectorResponse(
channel_energies, model="gaussian", sigma_a=0.005
)
energy = 2.0
resp_l = resp_low.get_response(energy)
resp_h = resp_high.get_response(energy)
resp_g = resp_gaussian_only.get_response(energy)
below_peak = channel_energies < energy
tail_low = resp_l[below_peak] - (1.0 - 0.1) * resp_g[below_peak]
tail_high = resp_h[below_peak] - (1.0 - 0.5) * resp_g[below_peak]
assert np.sum(tail_low) > 0, "Low tail_fraction should have positive tail"
assert np.sum(tail_high) > 0, "High tail_fraction should have positive tail"
ratio = np.sum(tail_high) / np.sum(tail_low)
assert ratio > 1.0, (
f"Higher tail_fraction should have more tail contribution, "
f"ratio={ratio:.2f}"
)
def test_tail_normalization(self):
"""Gaussian+tail response should still integrate to 1."""
channel_energies = np.linspace(0.0, 5.0, 8192)
resp = DetectorResponse(
channel_energies,
model="gaussian_tail",
sigma_a=0.005,
tail_fraction=0.2,
tau=0.01,
)
energy = 2.0
response = resp.get_response(energy)
integral = np.trapezoid(response, channel_energies)
assert np.isclose(
integral, 1.0, rtol=1e-2
), f"Gaussian+tail should normalize to 1, got {integral}"
class TestResolutionSigma:
"""Test energy-dependent resolution sigma(E)."""
def test_sigma_increases_with_energy(self):
"""σ(E) should increase with energy for sigma_b > 0."""
channel_energies = np.linspace(0.0, 5.0, 100)
resp = DetectorResponse(
channel_energies, model="gaussian", sigma_a=0.0, sigma_b=0.01
)
sigma_low = resp._resolution_sigma(0.5)
sigma_high = resp._resolution_sigma(4.0)
assert sigma_high > sigma_low, "σ should increase with √E"
def test_sigma_a_only(self):
"""With sigma_b=0, σ should be constant."""
channel_energies = np.linspace(0.0, 5.0, 100)
resp = DetectorResponse(
channel_energies,
model="gaussian",
sigma_a=0.01,
sigma_b=0.0,
fano_factor=0.0,
)
sigma_low = resp._resolution_sigma(0.5)
sigma_high = resp._resolution_sigma(4.0)
assert np.isclose(sigma_low, sigma_high), "σ should be constant when sigma_b=0"
def test_fano_contribution(self):
"""Fano factor should add energy-dependent contribution to σ."""
channel_energies = np.linspace(0.0, 5.0, 100)
resp_no_fano = DetectorResponse(
channel_energies,
model="gaussian",
sigma_a=0.0,
sigma_b=0.0,
fano_factor=0.0,
)
resp_fano = DetectorResponse(
channel_energies,
model="gaussian",
sigma_a=0.0,
sigma_b=0.0,
fano_factor=0.12,
)
sigma_low = resp_no_fano._resolution_sigma(0.5)
sigma_fano = resp_fano._resolution_sigma(0.5)
assert sigma_fano > sigma_low, "Fano factor should increase σ"
class TestTabulatedResponse:
"""Test tabulated (MC-simulated) detector response."""
def test_tabulated_interpolation(self):
"""Tabulated response should interpolate between calibration energies."""
channel_energies = np.linspace(0.0, 5.0, 256)
calib_energies = np.array([1.0, 2.0, 3.0, 4.0])
n_channels = len(channel_energies)
response_matrix = np.zeros((len(calib_energies), n_channels))
for i, e in enumerate(calib_energies):
idx = np.argmin(np.abs(channel_energies - e))
response_matrix[i, idx] = 1.0
resp = DetectorResponse(
channel_energies=channel_energies,
response_matrix=response_matrix,
calibration_energies=calib_energies,
)
response = resp.get_response(2.5)
assert np.any(response > 0), "Interpolated response should be non-zero"
def test_tabulated_outside_range(self):
"""Response outside calibration range should be zero."""
channel_energies = np.linspace(0.0, 5.0, 256)
calib_energies = np.array([1.0, 2.0, 3.0])
response_matrix = np.zeros((3, len(channel_energies)))
for i, e in enumerate(calib_energies):
idx = np.argmin(np.abs(channel_energies - e))
response_matrix[i, idx] = 1.0
resp = DetectorResponse(
channel_energies=channel_energies,
response_matrix=response_matrix,
calibration_energies=calib_energies,
)
response = resp.get_response(0.5)
assert np.all(response == 0), "Response below calibration range should be zero"
response = resp.get_response(4.0)
assert np.all(response == 0), "Response above calibration range should be zero"
class TestConvolution:
"""Test spectrum convolution with detector response."""
def test_convolution_preserves_total_counts(self):
"""Total counts should be approximately preserved after convolution."""
channel_energies = np.linspace(0.0, 5.0, 4096)
resp = DetectorResponse(channel_energies, model="gaussian", sigma_a=0.005)
W = np.linspace(0.0, 5.0, 1000)
spectrum = np.zeros_like(W)
idx = np.argmin(np.abs(W - 2.5))
spectrum[idx] = 1.0 / (W[1] - W[0])
convolved = resp.convolve(W, spectrum, normalize=True)
total_theoretical = np.trapezoid(spectrum, W)
total_convolved = np.trapezoid(convolved, channel_energies)
assert np.isclose(
total_convolved, total_theoretical, rtol=0.05
), f"Total counts: theory={total_theoretical:.4f}, convolved={total_convolved:.4f}"
def test_convolution_widens_peak(self):
"""Convolution should widen a narrow peak."""
channel_energies = np.linspace(0.0, 5.0, 4096)
resp = DetectorResponse(channel_energies, model="gaussian", sigma_a=0.02)
W = np.linspace(0.0, 5.0, 2000)
spectrum = np.zeros_like(W)
idx = np.argmin(np.abs(W - 2.5))
spectrum[idx] = 1.0 / (W[1] - W[0])
convolved = resp.convolve(W, spectrum, normalize=True)
peak_width = np.sum(convolved > np.max(convolved) / 2)
assert peak_width > 1, "Convolved peak should have finite width"
def test_convolution_beta_spectrum(self):
"""Convolution of a beta spectrum should produce a smooth curve."""
channel_energies = np.linspace(0.0, 5.0, 4096)
resp = DetectorResponse(channel_energies, model="gaussian", sigma_a=0.01)
W0 = 5.0
W = np.linspace(1.0, W0 - 0.01, 500)
p = np.sqrt(W**2 - 1)
spectrum = p * W * (W0 - W) ** 2
convolved = resp.convolve(W, spectrum, normalize=True)
assert len(convolved) == len(channel_energies)
assert np.all(convolved >= 0), "Convolved spectrum should be non-negative"
assert np.max(convolved) > 0, "Convolved spectrum should have non-zero peak"
def test_convolve_batch_same_as_convolve(self):
"""Batch convolution should match standard convolution for analytical models."""
channel_energies = np.linspace(0.0, 5.0, 4096)
resp = DetectorResponse(
channel_energies,
model="gaussian_tail",
sigma_a=0.01,
tail_fraction=0.1,
tau=0.01,
)
W = np.linspace(1.0, 4.9, 300)
p = np.sqrt(W**2 - 1)
spectrum = p * W * (4.9 - W) ** 2
result_standard = resp.convolve(W, spectrum, normalize=True)
result_batch = resp.convolve_batch(W, spectrum, normalize=True)
assert np.allclose(
result_standard, result_batch, rtol=1e-5
), "Batch and standard convolution should match for analytical models"
tests/physics/test_fitter.py
+448
View File
@@ -0,0 +1,448 @@
# tests/test_fitter.py — Tests for χ² minimization framework
"""
Why test this?
--------------
The fitter is the core numerical tool for extracting physical parameters
from beta spectra. It must be correct because all downstream results
(g_V, g_A, C(W) extraction) depend on it.
We test:
1. Basic χ² minimization on known functions (quadratic, linear).
2. Covariance matrix estimation (should match input uncertainties).
3. Confidence interval computation via profile likelihood.
4. Goodness-of-fit metrics (χ²/ndf, p-value).
5. Bounds handling and convergence diagnostics.
6. Robust loss functions (soft_l1, huber) for outlier resistance.
"""
import numpy as np
from beta_spectrum import CurveFitter, FitConfig
class TestBasicFitting:
"""Test χ² minimization on simple known functions."""
def test_quadratic_fit(self):
"""Fit y = a + b·x + c·x² with known coefficients."""
np.random.seed(42)
x = np.linspace(0, 10, 50)
true_params = [1.0, 2.0, -0.1]
y = true_params[0] + true_params[1] * x + true_params[2] * x**2
y += np.random.normal(0, 0.5, len(x))
def model(x, a, b, c):
return a + b * x + c * x**2
fitter = CurveFitter(model, x, y, uncertainties=np.ones_like(y) * 0.5)
result = fitter.fit([0.0, 0.0, 0.0])
for i, name in enumerate(["a", "b", "c"]):
assert (
abs(result.parameters[i] - true_params[i]) < 0.5
), f"Parameter {name}: expected {true_params[i]}, got {result.parameters[i]}"
assert result.success, "Fit should converge"
assert result.chi2 > 0, "χ² must be positive"
def test_linear_fit(self):
"""Fit y = a + b·x with known coefficients."""
np.random.seed(123)
x = np.linspace(0, 5, 30)
true_a, true_b = 3.0, -1.5
y = true_a + true_b * x
y += np.random.normal(0, 0.2, len(x))
def model(x, a, b):
return a + b * x
fitter = CurveFitter(model, x, y, uncertainties=np.ones_like(y) * 0.2)
result = fitter.fit([0.0, 0.0])
assert abs(result.parameters[0] - true_a) < 0.3
assert abs(result.parameters[1] - true_b) < 0.3
assert result.success
def test_constant_fit(self):
"""Fit y = a (constant) with known value."""
np.random.seed(456)
x = np.linspace(0, 10, 40)
true_a = 5.0
y = np.full_like(x, true_a)
y += np.random.normal(0, 0.3, len(x))
def model(x, a):
return np.full_like(x, a)
fitter = CurveFitter(model, x, y, uncertainties=np.ones_like(x) * 0.3)
result = fitter.fit([0.0])
assert abs(result.parameters[0] - true_a) < 0.2
assert result.success
class TestCovarianceEstimation:
"""Test that covariance matrix is correctly estimated."""
def test_covariance_matches_uncertainty(self):
"""For a constant model, parameter error should match input σ."""
np.random.seed(789)
x = np.linspace(0, 10, 50)
true_a = 5.0
sigma = 0.5
y = np.full_like(x, true_a, dtype=float)
y += np.random.normal(0, sigma, len(x))
def model(x, a):
return np.full_like(x, a)
fitter = CurveFitter(model, x, y, uncertainties=np.ones_like(x) * sigma)
result = fitter.fit([0.0])
error = np.sqrt(result.covariance[0, 0])
# Expected error: σ / sqrt(N)
expected_error = sigma / np.sqrt(len(x))
assert (
abs(error - expected_error) < expected_error * 0.3
), f"Covariance error {error:.4f} != expected {expected_error:.4f}"
def test_covariance_matrix_positive_definite(self):
"""Covariance matrix must be positive definite."""
np.random.seed(321)
x = np.linspace(0, 5, 30)
y = 2.0 + 1.5 * x - 0.3 * x**2
y += np.random.normal(0, 0.2, len(x))
def model(x, a, b, c):
return a + b * x + c * x**2
fitter = CurveFitter(model, x, y, uncertainties=np.ones_like(x) * 0.2)
result = fitter.fit([0.0, 0.0, 0.0])
# Check eigenvalues are non-negative
eigvals = np.linalg.eigvalsh(result.covariance)
assert np.all(eigvals >= -1e-10), f"Covariance not positive definite: {eigvals}"
def test_correlation_matrix(self):
"""Correlation matrix should have 1.0 on diagonal."""
np.random.seed(654)
x = np.linspace(0, 10, 50)
y = 1.0 + 2.0 * x + 0.1 * x**2
y += np.random.normal(0, 0.5, len(x))
def model(x, a, b, c):
return a + b * x + c * x**2
fitter = CurveFitter(model, x, y, uncertainties=np.ones_like(x) * 0.5)
result = fitter.fit([0.0, 0.0, 0.0])
diag = np.diag(result.correlation_matrix)
assert np.allclose(
diag, 1.0, atol=1e-6
), f"Diagonal of correlation matrix: {diag}"
class TestGoodnessOfFit:
"""Test χ²/ndf and p-value calculations."""
def test_chi2_per_dof(self):
"""χ²/ndof should be positive and reasonable."""
np.random.seed(111)
x = np.linspace(0, 10, 50)
y = 2.0 * x
y += np.random.normal(0, 0.3, len(x))
def model(x, a):
return a * x
fitter = CurveFitter(model, x, y, uncertainties=np.ones_like(x) * 0.3)
result = fitter.fit([0.0])
assert result.chi2_per_dof > 0
assert result.n_points == 50
assert result.n_free == 1
assert result.n_free == 1
def test_p_value_range(self):
"""p-value should be in [0, 1]."""
np.random.seed(222)
x = np.linspace(0, 5, 30)
y = 1.0 + x
y += np.random.normal(0, 0.1, len(x))
def model(x, a, b):
return a + b * x
fitter = CurveFitter(model, x, y, uncertainties=np.ones_like(x) * 0.1)
result = fitter.fit([0.0, 0.0])
assert 0.0 <= result.p_value <= 1.0
def test_p_value_good_fit(self):
"""Good fit should have reasonable p-value."""
np.random.seed(333)
x = np.linspace(0, 10, 100)
y = 3.0 + 1.5 * x - 0.05 * x**2
y += np.random.normal(0, 0.2, len(x))
def model(x, a, b, c):
return a + b * x + c * x**2
fitter = CurveFitter(model, x, y, uncertainties=np.ones_like(x) * 0.2)
result = fitter.fit([0.0, 0.0, 0.0])
assert (
result.p_value > 0.001
), f"Good fit should have reasonable p-value, got {result.p_value}"
class TestBounds:
"""Test parameter bounds handling."""
def test_lower_bound(self):
"""Test that lower bounds are respected."""
np.random.seed(444)
x = np.linspace(0, 5, 30)
y = 2.0 + 1.0 * x
y += np.random.normal(0, 0.1, len(x))
def model(x, a, b):
return a + b * x
fitter = CurveFitter(model, x, y, uncertainties=np.ones_like(x) * 0.1)
result = fitter.fit([1.5, 0.5], bounds=([1.0, 0.0], [10.0, 5.0]))
assert (
result.parameters[0] >= 1.0
), f"Parameter a should be >= 1.0, got {result.parameters[0]}"
assert (
result.parameters[1] >= 0.0
), f"Parameter b should be >= 0.0, got {result.parameters[1]}"
def test_upper_bound(self):
"""Test that upper bounds are respected."""
np.random.seed(555)
x = np.linspace(0, 5, 30)
y = 0.5 + 0.5 * x
y += np.random.normal(0, 0.1, len(x))
def model(x, a, b):
return a + b * x
fitter = CurveFitter(model, x, y, uncertainties=np.ones_like(x) * 0.1)
result = fitter.fit([0.5, 0.5], bounds=([0.0, 0.0], [1.0, 1.0]))
assert result.parameters[0] <= 1.0
assert result.parameters[1] <= 1.0
class TestRobustLoss:
"""Test robust loss functions for outlier resistance."""
def test_soft_l1_loss(self):
"""soft_l1 loss should handle outliers better than linear."""
np.random.seed(666)
x = np.linspace(0, 10, 50)
y = 1.0 + 2.0 * x
y[25] = 50.0 # Add outlier
def model(x, a, b):
return a + b * x
# Linear loss — outlier pulls fit away
config_linear = FitConfig(loss="linear")
fitter_linear = CurveFitter(
model, x, y, uncertainties=np.ones_like(x) * 1.0, config=config_linear
)
result_linear = fitter_linear.fit([0.0, 0.0])
# soft_l1 loss — outlier has less influence
config_robust = FitConfig(loss="soft_l1")
fitter_robust = CurveFitter(
model, x, y, uncertainties=np.ones_like(x) * 1.0, config=config_robust
)
result_robust = fitter_robust.fit([0.0, 0.0])
# Robust fit should be closer to true slope
assert (
abs(result_robust.parameters[1] - 2.0)
< abs(result_linear.parameters[1] - 2.0) * 1.5
), "Robust loss should handle outliers better"
class TestFitResult:
"""Test FitResult dataclass methods."""
def test_parameters_with_errors(self):
"""parameters_with_errors should return dict of (value, error)."""
np.random.seed(777)
x = np.linspace(0, 5, 30)
y = 2.0 + x
y += np.random.normal(0, 0.1, len(x))
def model(x, a, b):
return a + b * x
fitter = CurveFitter(model, x, y, uncertainties=np.ones_like(x) * 0.1)
result = fitter.fit([0.0, 0.0])
errors = result.parameters_with_errors
assert len(errors) == 2
for key, (val, err) in errors.items():
assert isinstance(val, float)
assert isinstance(err, float)
assert err >= 0
def test_summary(self):
"""summary() should return a non-empty string."""
np.random.seed(888)
x = np.linspace(0, 5, 30)
y = 1.0 + 2.0 * x
y += np.random.normal(0, 0.1, len(x))
def model(x, a, b):
return a + b * x
fitter = CurveFitter(model, x, y, uncertainties=np.ones_like(x) * 0.1)
result = fitter.fit([0.0, 0.0])
summary = result.summary(param_names=["a", "b"])
assert "Fit Results" in summary
assert "a" in summary
assert "b" in summary
assert "χ²" in summary
assert "Converged" in summary
def test_residuals_property(self):
"""residuals should be data - model."""
np.random.seed(999)
x = np.linspace(0, 5, 30)
y = 2.0 + x
y += np.random.normal(0, 0.1, len(x))
def model(x, a, b):
return a + b * x
fitter = CurveFitter(model, x, y, uncertainties=np.ones_like(x) * 0.1)
result = fitter.fit([0.0, 0.0])
expected_residuals = y - result.model_values
np.testing.assert_allclose(result.residuals, expected_residuals)
class TestFitConfig:
"""Test FitConfig options."""
def test_custom_config(self):
"""Custom FitConfig should be respected."""
np.random.seed(101)
x = np.linspace(0, 5, 30)
y = 1.0 + x
y += np.random.normal(0, 0.1, len(x))
def model(x, a, b):
return a + b * x
config = FitConfig(
method="trf",
max_nfev=500,
loss="linear",
ftol=1e-10,
)
fitter = CurveFitter(
model, x, y, uncertainties=np.ones_like(x) * 0.1, config=config
)
result = fitter.fit([0.0, 0.0])
assert result.success
class TestConfidenceIntervals:
"""Test confidence interval computation."""
def test_95_percent_ci(self):
"""95% CI should contain the true parameter value."""
np.random.seed(202)
x = np.linspace(0, 10, 60)
true_a, true_b = 2.0, 1.5
y = true_a + true_b * x
y += np.random.normal(0, 0.3, len(x))
def model(x, a, b):
return a + b * x
fitter = CurveFitter(model, x, y, uncertainties=np.ones_like(x) * 0.3)
result = fitter.fit([0.0, 0.0])
ci_a = result.confidence_interval(0, confidence=0.95)
ci_b = result.confidence_interval(1, confidence=0.95)
assert ci_a[0] <= true_a <= ci_a[1], f"True a={true_a} not in CI {ci_a}"
assert ci_b[0] <= true_b <= ci_b[1], f"True b={true_b} not in CI {ci_b}"
def test_ci_width_increases_with_confidence(self):
"""Higher confidence level → wider interval."""
np.random.seed(303)
x = np.linspace(0, 5, 40)
y = 1.0 + 2.0 * x
y += np.random.normal(0, 0.2, len(x))
def model(x, a, b):
return a + b * x
fitter = CurveFitter(model, x, y, uncertainties=np.ones_like(x) * 0.2)
result = fitter.fit([0.0, 0.0])
ci_90 = result.confidence_interval(0, confidence=0.90)
ci_95 = result.confidence_interval(0, confidence=0.95)
width_90 = ci_90[1] - ci_90[0]
width_95 = ci_95[1] - ci_95[0]
assert (
width_95 > width_90
), f"95% CI width {width_95} should be > 90% width {width_90}"
class TestEdgeCases:
"""Test edge cases and error handling."""
def test_zero_uncertainties_filtered(self):
"""Points with zero uncertainty should be filtered out."""
x = np.array([1.0, 2.0, 3.0, 4.0])
y = np.array([2.0, 4.0, 6.0, 8.0])
unc = np.array([1.0, 0.0, 1.0, 0.0])
def model(x, a):
return a * x
fitter = CurveFitter(model, x, y, uncertainties=unc)
assert len(fitter.x_data) == 2, "Zero-uncertainty points should be removed"
assert len(fitter.y_data) == 2
def test_negative_uncertainties_filtered(self):
"""Points with negative uncertainty should be filtered out."""
x = np.array([1.0, 2.0, 3.0])
y = np.array([2.0, 4.0, 6.0])
unc = np.array([1.0, -1.0, 1.0])
def model(x, a):
return a * x
fitter = CurveFitter(model, x, y, uncertainties=unc)
assert len(fitter.x_data) == 2
def test_unbounded_parameters(self):
"""Unbounded parameters (±inf) should work."""
np.random.seed(404)
x = np.linspace(0, 5, 30)
y = 1.0 + 2.0 * x
y += np.random.normal(0, 0.1, len(x))
def model(x, a, b):
return a + b * x
fitter = CurveFitter(model, x, y, uncertainties=np.ones_like(x) * 0.1)
result = fitter.fit([0.0, 0.0], bounds=([-np.inf, -np.inf], [np.inf, np.inf]))
assert result.success
tests/physics/test_multi_branch.py
+783
View File
@@ -0,0 +1,783 @@
"""Tests for multi-branch decay support.
Why test this?
--------------
Multi-branch decay is a critical feature for realistic nuclear decay analysis.
Many isotopes decay through multiple channels to different daughter states,
each with its own endpoint energy, transition type, and intensity.
The multi-branch mode is the DEFAULT — even single-branch decays are handled
through the same code path. The only difference in behavior is in plot generation
(where a single branch doesn't need separate branch plots).
We test:
1. **BranchConfig validation**: intensity bounds, positive endpoint.
2. **Single-branch via multi-branch path**: one branch should give identical
results to the old single-branch behavior.
3. **Multi-branch spectrum calculation**: weighted sum produces valid output.
4. **Intensity weighting**: sum of intensities = 1.0, weighted sum = total.
5. **Per-branch spectra**: each branch has correct endpoint and transition type.
6. **Component extraction**: universal + per-branch components correctly separated.
7. **CSV export**: multi-branch CSV has correct column structure.
8. **Energy grid**: uses max endpoint across all branches.
9. **Detector response**: multi-branch spectra are convolvable.
10. **BranchConfig dataclass**: basic construction and validation.
"""
import numpy as np
import pytest
from beta_spectrum import (
BetaSpectrum,
BetaSpectrumAnalyzer,
BranchConfig,
SpectrumConfig,
)
from beta_spectrum.utils import T_to_W
# ---------------------------------------------------------------------------
# BranchConfig validation
# ---------------------------------------------------------------------------
class TestBranchConfigValidation:
"""Test BranchConfig dataclass validation."""
def test_valid_branch(self):
"""Valid branch config should construct without error."""
branch = BranchConfig(
endpoint_MeV=0.2,
transition_type="A",
intensity=0.5,
)
assert branch.endpoint_MeV == 0.2
assert branch.transition_type == "A"
assert branch.intensity == 0.5
def test_intensity_out_of_range_low(self):
"""Intensity below 0 should raise ValueError."""
with pytest.raises(ValueError, match="intensity must be in"):
BranchConfig(endpoint_MeV=0.2, intensity=-0.1)
def test_intensity_out_of_range_high(self):
"""Intensity above 1.0 should raise ValueError."""
with pytest.raises(ValueError, match="intensity must be in"):
BranchConfig(endpoint_MeV=0.2, intensity=1.1)
def test_zero_endpoint(self):
"""Zero endpoint should raise ValueError."""
with pytest.raises(ValueError, match="endpoint_MeV must be positive"):
BranchConfig(endpoint_MeV=0.0, intensity=0.5)
def test_negative_endpoint(self):
"""Negative endpoint should raise ValueError."""
with pytest.raises(ValueError, match="endpoint_MeV must be positive"):
BranchConfig(endpoint_MeV=-0.1, intensity=0.5)
def test_default_transition_type(self):
"""Transition type defaults to 'A'."""
branch = BranchConfig(endpoint_MeV=0.2, intensity=0.5)
assert branch.transition_type == "A"
# ---------------------------------------------------------------------------
# Single-branch via multi-branch path (backward compatibility)
# ---------------------------------------------------------------------------
class TestSingleBranchViaMultiBranch:
"""Verify single-branch mode works through the multi-branch path."""
def test_one_branch_produces_valid_spectrum(self):
"""One branch should produce valid spectrum identical to old behavior."""
config = SpectrumConfig(
Z_parent=43,
Z_daughter=44,
A_number=99,
endpoint_MeV=0.294,
branches=[
BranchConfig(
endpoint_MeV=0.294,
transition_type="F2",
intensity=1.0,
)
],
)
spectrum = BetaSpectrum.from_config(config)
W, E = spectrum.get_energy_grid(config)
values = spectrum(W)
assert len(values) == len(E)
assert np.all(np.isfinite(values))
assert np.all(values >= 0), "Spectrum should be non-negative"
def test_one_branch_has_no_branch_spectra(self):
"""Single branch should not produce separate branch plots."""
config = SpectrumConfig(
Z_parent=43,
Z_daughter=44,
A_number=99,
endpoint_MeV=0.294,
branches=[
BranchConfig(endpoint_MeV=0.294, transition_type="F2", intensity=1.0),
],
)
spectrum = BetaSpectrum.from_config(config)
W, _ = spectrum.get_energy_grid(config)
branch_spectra = spectrum.get_branch_spectra(W)
# Single branch: no separate branch spectra needed
assert len(branch_spectra) == 0
def test_one_branch_analyzer_uses_normal_plot(self):
"""Single branch analyzer should not have multi-branch mode."""
config = SpectrumConfig(
Z_parent=43,
Z_daughter=44,
A_number=99,
endpoint_MeV=0.294,
branches=[
BranchConfig(endpoint_MeV=0.294, transition_type="F2", intensity=1.0),
],
)
spectrum = BetaSpectrum.from_config(config)
analyzer = BetaSpectrumAnalyzer(spectrum, config)
assert not analyzer._is_multi
# ---------------------------------------------------------------------------
# Multi-branch spectrum calculation
# ---------------------------------------------------------------------------
class TestMultiBranchSpectrum:
"""Test multi-branch spectrum calculation."""
def test_two_branches_basic(self):
"""Two branches should produce valid weighted sum."""
config = SpectrumConfig(
Z_parent=43,
Z_daughter=44,
A_number=99,
endpoint_MeV=0.350,
branches=[
BranchConfig(endpoint_MeV=0.294, transition_type="F2", intensity=0.7),
BranchConfig(endpoint_MeV=0.150, transition_type="A", intensity=0.3),
],
)
spectrum = BetaSpectrum.from_config(config)
assert spectrum._is_multi
assert len(spectrum.branch_spectra) == 2
assert len(spectrum.branches) == 2
W, E = spectrum.get_energy_grid(config)
values = spectrum(W)
assert len(values) == len(E)
assert np.all(np.isfinite(values))
assert np.all(values >= 0), "Spectrum should be non-negative"
def test_three_branches(self):
"""Three branches should produce valid weighted sum."""
config = SpectrumConfig(
Z_parent=27,
Z_daughter=28,
A_number=60,
endpoint_MeV=0.320,
branches=[
BranchConfig(endpoint_MeV=0.310, transition_type="A", intensity=0.5),
BranchConfig(endpoint_MeV=0.200, transition_type="F1", intensity=0.3),
BranchConfig(endpoint_MeV=0.100, transition_type="A", intensity=0.2),
],
)
spectrum = BetaSpectrum.from_config(config)
W, E = spectrum.get_energy_grid(config)
values = spectrum(W)
assert spectrum._is_multi
assert len(values) == len(E)
assert np.all(np.isfinite(values))
assert np.all(values >= 0)
def test_energy_grid_uses_max_endpoint(self):
"""Energy grid should extend to the maximum branch endpoint."""
config = SpectrumConfig(
Z_parent=43,
Z_daughter=44,
A_number=99,
endpoint_MeV=0.350,
e_step_MeV=0.001,
branches=[
BranchConfig(endpoint_MeV=0.300, transition_type="A", intensity=0.5),
BranchConfig(endpoint_MeV=0.200, transition_type="F2", intensity=0.5),
],
)
spectrum = BetaSpectrum.from_config(config)
W, E = spectrum.get_energy_grid(config)
# Grid should extend to max endpoint (0.300 MeV)
assert E[-1] >= 0.299, f"Grid should reach ~0.300 MeV, got {E[-1]}"
def test_branch_spectra_property(self):
"""get_branch_spectra() should return per-branch unweighted spectra."""
config = SpectrumConfig(
Z_parent=43,
Z_daughter=44,
A_number=99,
endpoint_MeV=0.350,
branches=[
BranchConfig(endpoint_MeV=0.294, transition_type="F2", intensity=0.7),
BranchConfig(endpoint_MeV=0.150, transition_type="A", intensity=0.3),
],
)
spectrum = BetaSpectrum.from_config(config)
W, _ = spectrum.get_energy_grid(config)
branch_spectra = spectrum.get_branch_spectra(W)
assert len(branch_spectra) == 2
assert len(branch_spectra[0]) == len(W)
assert len(branch_spectra[1]) == len(W)
def test_branch_normalized_spectra(self):
"""get_branch_normalized_spectra() should return intensity-weighted spectra."""
config = SpectrumConfig(
Z_parent=43,
Z_daughter=44,
A_number=99,
endpoint_MeV=0.350,
# Disable universal components to simplify
use_fermi=False,
use_screening=False,
use_finite_size=False,
use_charge_dist=False,
use_exchange=False,
branches=[
BranchConfig(endpoint_MeV=0.294, transition_type="F2", intensity=0.7),
BranchConfig(endpoint_MeV=0.150, transition_type="A", intensity=0.3),
],
)
spectrum = BetaSpectrum.from_config(config)
W, _ = spectrum.get_energy_grid(config)
weighted = spectrum.get_branch_normalized_spectra(W)
assert len(weighted) == 2
# With no universal components, total = sum of weighted branch spectra
total_from_branches = weighted[0] + weighted[1]
total = spectrum(W)
assert np.allclose(total, total_from_branches, rtol=1e-10)
def test_intensity_normalization(self):
"""Branch intensities should be used as weights correctly."""
config = SpectrumConfig(
Z_parent=43,
Z_daughter=44,
A_number=99,
endpoint_MeV=0.350,
use_phase_space=True,
use_fermi=False,
use_screening=False,
use_finite_size=False,
use_charge_dist=False,
use_exchange=False,
use_radiative=False,
branches=[
BranchConfig(endpoint_MeV=0.294, transition_type="A", intensity=0.8),
BranchConfig(endpoint_MeV=0.294, transition_type="A", intensity=0.2),
],
)
spectrum = BetaSpectrum.from_config(config)
W, _ = spectrum.get_energy_grid(config)
# Both branches have same endpoint and transition type,
# so phase space should be identical. Total = 0.8*PS + 0.2*PS = PS
total = spectrum(W)
branch0 = spectrum.branch_spectra[0](W)
branch1 = spectrum.branch_spectra[1](W)
expected = 0.8 * branch0 + 0.2 * branch1
assert np.allclose(total, expected, rtol=1e-10)
# ---------------------------------------------------------------------------
# Component extraction in multi-branch mode
# ---------------------------------------------------------------------------
class TestMultiBranchComponents:
"""Test component extraction for multi-branch mode."""
def test_universal_components_present(self):
"""Universal components (Fermi, etc.) should be in component dict."""
config = SpectrumConfig(
Z_parent=43,
Z_daughter=44,
A_number=99,
endpoint_MeV=0.350,
branches=[
BranchConfig(endpoint_MeV=0.294, transition_type="F2", intensity=0.7),
BranchConfig(endpoint_MeV=0.150, transition_type="A", intensity=0.3),
],
)
spectrum = BetaSpectrum.from_config(config)
W, _ = spectrum.get_energy_grid(config)
components = spectrum.calculate_components(W)
assert "Fermi" in components, "Universal Fermi component should be present"
assert "Screening" in components
assert "FiniteSizeL0" in components
def test_per_branch_components_present(self):
"""Per-branch components (PhaseSpace, Radiative) should be in component dict."""
config = SpectrumConfig(
Z_parent=43,
Z_daughter=44,
A_number=99,
endpoint_MeV=0.350,
branches=[
BranchConfig(endpoint_MeV=0.294, transition_type="F2", intensity=0.7),
BranchConfig(endpoint_MeV=0.150, transition_type="A", intensity=0.3),
],
)
spectrum = BetaSpectrum.from_config(config)
W, _ = spectrum.get_energy_grid(config)
components = spectrum.calculate_components(W)
assert "branch_0.PhaseSpace" in components
assert "branch_1.PhaseSpace" in components
assert "branch_0.Radiative" in components
assert "branch_1.Radiative" in components
def test_analyzer_components_multi_branch(self):
"""BetaSpectrumAnalyzer should handle multi-branch components."""
config = SpectrumConfig(
Z_parent=43,
Z_daughter=44,
A_number=99,
endpoint_MeV=0.350,
branches=[
BranchConfig(endpoint_MeV=0.294, transition_type="F2", intensity=0.7),
BranchConfig(endpoint_MeV=0.150, transition_type="A", intensity=0.3),
],
)
spectrum = BetaSpectrum.from_config(config)
analyzer = BetaSpectrumAnalyzer(spectrum, config)
components = analyzer.components
assert "Fermi" in components
assert "branch_0.PhaseSpace" in components
def test_analyzer_branch_spectra_property(self):
"""Analyzer should expose branch spectra property."""
config = SpectrumConfig(
Z_parent=43,
Z_daughter=44,
A_number=99,
endpoint_MeV=0.350,
branches=[
BranchConfig(endpoint_MeV=0.294, transition_type="F2", intensity=0.7),
BranchConfig(endpoint_MeV=0.150, transition_type="A", intensity=0.3),
],
)
spectrum = BetaSpectrum.from_config(config)
analyzer = BetaSpectrumAnalyzer(spectrum, config)
branch_spectra = analyzer.branch_spectra
assert len(branch_spectra) == 2
def test_total_equals_universal_times_branch_sum(self):
"""Total = universal_components * sum(intensity * branch_spectra)."""
config = SpectrumConfig(
Z_parent=43,
Z_daughter=44,
A_number=99,
endpoint_MeV=0.350,
branches=[
BranchConfig(endpoint_MeV=0.294, transition_type="F2", intensity=0.7),
BranchConfig(endpoint_MeV=0.150, transition_type="A", intensity=0.3),
],
)
spectrum = BetaSpectrum.from_config(config)
W, _ = spectrum.get_energy_grid(config)
total = spectrum(W)
# Compute expected: universal_components * sum(intensity * branch_spectra)
# Each branch spectrum must be masked to zero beyond its endpoint
universal_product = np.ones_like(W)
for comp in spectrum.components:
universal_product *= comp(W)
branch_sum = np.zeros_like(W)
for i, spec in enumerate(spectrum.branch_spectra):
W0_branch = T_to_W(spectrum.branches[i].endpoint_MeV)
mask = W <= W0_branch
branch_sum += spectrum.branches[i].intensity * np.where(mask, spec(W), 0.0)
expected = universal_product * branch_sum
assert np.allclose(total, expected, rtol=1e-10)
# ---------------------------------------------------------------------------
# CSV export in multi-branch mode
# ---------------------------------------------------------------------------
class TestMultiBranchCSVExport:
"""Test CSV export structure for multi-branch mode."""
def _write_json(self, data):
import json
import tempfile
from pathlib import Path
tmp = tempfile.NamedTemporaryFile(mode="w", suffix=".json", delete=False)
json.dump(data, tmp)
tmp.close()
return Path(tmp.name)
def test_multi_branch_csv_columns(self, tmp_path):
"""Multi-branch CSV should have total + per-branch + universal columns."""
config = SpectrumConfig(
Z_parent=43,
Z_daughter=44,
A_number=99,
endpoint_MeV=0.350,
branches=[
BranchConfig(endpoint_MeV=0.294, transition_type="F2", intensity=0.7),
BranchConfig(endpoint_MeV=0.150, transition_type="A", intensity=0.3),
],
)
spectrum = BetaSpectrum.from_config(config)
analyzer = BetaSpectrumAnalyzer(spectrum, config)
csv_path = str(tmp_path / "multi_branch.csv")
analyzer.export_to_csv(csv_path)
import pandas as pd
# Read CSV (skip comment headers)
df = pd.read_csv(csv_path, comment="#")
# Should have total spectrum column
assert "spectrum" in df.columns
# Should have per-branch spectrum columns (intensity/endpoint/transition are in the header)
assert "branch_1_spectrum" in df.columns
assert "branch_2_spectrum" in df.columns
# Should NOT have repeated per-branch metadata columns
assert "branch_1_intensity" not in df.columns
assert "branch_1_endpoint" not in df.columns
assert "branch_1_transition" not in df.columns
assert "branch_2_intensity" not in df.columns
assert "branch_2_endpoint" not in df.columns
assert "branch_2_transition" not in df.columns
# Should have universal columns
assert "Fermi" in df.columns
def test_single_branch_csv_columns(self, tmp_path):
"""Single-branch CSV should have standard columns."""
config = SpectrumConfig(
Z_parent=43,
Z_daughter=44,
A_number=99,
endpoint_MeV=0.294,
branches=[
BranchConfig(endpoint_MeV=0.294, transition_type="F2", intensity=1.0),
],
)
spectrum = BetaSpectrum.from_config(config)
analyzer = BetaSpectrumAnalyzer(spectrum, config)
csv_path = str(tmp_path / "single_branch.csv")
analyzer.export_to_csv(csv_path)
import pandas as pd
df = pd.read_csv(csv_path, comment="#")
assert "spectrum" in df.columns
assert "Fermi" in df.columns
# ---------------------------------------------------------------------------
# Multi-branch with different transition types
# ---------------------------------------------------------------------------
class TestMultiBranchTransitionTypes:
"""Test that different transition types are correctly propagated."""
def test_different_transition_types(self):
"""Each branch should use its own transition type for PhaseSpace."""
config = SpectrumConfig(
Z_parent=43,
Z_daughter=44,
A_number=99,
endpoint_MeV=0.350,
branches=[
BranchConfig(endpoint_MeV=0.294, transition_type="F2", intensity=0.7),
BranchConfig(endpoint_MeV=0.200, transition_type="A", intensity=0.3),
],
)
spectrum = BetaSpectrum.from_config(config)
# Check that each branch has the correct transition type
assert spectrum.branches[0].transition_type == "F2"
assert spectrum.branches[1].transition_type == "A"
# PhaseSpace components should have different W0
W, _ = spectrum.get_energy_grid(config)
components = spectrum.calculate_components(W)
branch0_ps = components["branch_0.PhaseSpace"]
branch1_ps = components["branch_1.PhaseSpace"]
# Branch 0 has higher endpoint, so its phase space extends further
assert branch0_ps.shape == branch1_ps.shape
def test_branch_info_transition_type_from_ensdf(self):
"""BranchInfo should have transition_type field populated from ENSDF."""
from beta_spectrum.nuclear_data import BranchInfo
branch = BranchInfo(
level_index=1,
level_energy_keV=100.0,
intensity=0.5,
log_ft=5.0,
transition_type="F1",
)
assert branch.transition_type == "F1"
assert branch.level_index == 1
assert branch.level_energy_keV == 100.0
assert branch.intensity == 0.5
def test_all_branches_same_transition(self):
"""All branches can have the same transition type."""
config = SpectrumConfig(
Z_parent=27,
Z_daughter=28,
A_number=60,
endpoint_MeV=0.320,
branches=[
BranchConfig(endpoint_MeV=0.310, transition_type="A", intensity=0.5),
BranchConfig(endpoint_MeV=0.200, transition_type="A", intensity=0.3),
BranchConfig(endpoint_MeV=0.100, transition_type="A", intensity=0.2),
],
)
spectrum = BetaSpectrum.from_config(config)
W, E = spectrum.get_energy_grid(config)
values = spectrum(W)
assert len(values) == len(E)
assert np.all(np.isfinite(values))
assert np.all(values >= 0)
# ---------------------------------------------------------------------------
# Multi-branch with detector response
# ---------------------------------------------------------------------------
class TestMultiBranchDetectorResponse:
"""Test multi-branch with detector response convolution."""
def test_multi_branch_with_detector(self):
"""Multi-branch spectrum should be convolvable with detector response."""
config = SpectrumConfig(
Z_parent=43,
Z_daughter=44,
A_number=99,
endpoint_MeV=0.350,
use_detector_response=True,
detector_model="gaussian",
detector_sigma_a_keV=1.0,
detector_n_channels=256,
detector_fano_factor=0.0,
branches=[
BranchConfig(endpoint_MeV=0.294, transition_type="F2", intensity=0.7),
BranchConfig(endpoint_MeV=0.150, transition_type="A", intensity=0.3),
],
)
spectrum = BetaSpectrum.from_config(config)
W, _ = spectrum.get_energy_grid(config)
detector = BetaSpectrum.create_detector_from_config(config)
convolved = spectrum.convolve_with_detector(detector, W=W, config=config)
assert convolved.shape == (256,)
assert np.all(convolved >= 0), "Convolved spectrum should be non-negative"
assert np.sum(convolved) > 0, "Convolved spectrum should have non-zero area"
# ---------------------------------------------------------------------------
# Edge cases
# ---------------------------------------------------------------------------
class TestMultiBranchEdgeCases:
"""Test edge cases for multi-branch mode."""
def test_very_low_intensity_branch(self):
"""Branch with very low intensity should still be computed."""
config = SpectrumConfig(
Z_parent=43,
Z_daughter=44,
A_number=99,
endpoint_MeV=0.350,
branches=[
BranchConfig(endpoint_MeV=0.294, transition_type="F2", intensity=0.999),
BranchConfig(endpoint_MeV=0.100, transition_type="A", intensity=0.001),
],
)
spectrum = BetaSpectrum.from_config(config)
W, _ = spectrum.get_energy_grid(config)
values = spectrum(W)
assert np.all(np.isfinite(values))
def test_very_high_z_multi_branch(self):
"""Multi-branch with high Z should be numerically stable."""
config = SpectrumConfig(
Z_parent=90,
Z_daughter=91,
A_number=232,
endpoint_MeV=1.0,
branches=[
BranchConfig(endpoint_MeV=0.8, transition_type="A", intensity=0.6),
BranchConfig(endpoint_MeV=0.5, transition_type="F1", intensity=0.4),
],
)
spectrum = BetaSpectrum.from_config(config)
W, _ = spectrum.get_energy_grid(config)
values = spectrum(W)
assert np.all(
np.isfinite(values)
), "High-Z multi-branch must not produce NaN/inf"
def test_no_branches_defaults_to_single(self):
"""No branches list should default to single-branch mode."""
config = SpectrumConfig(
Z_parent=43,
Z_daughter=44,
A_number=99,
endpoint_MeV=0.294,
)
spectrum = BetaSpectrum.from_config(config)
W, E = spectrum.get_energy_grid(config)
values = spectrum(W)
assert len(values) == len(E)
assert np.all(np.isfinite(values))
def test_branch_values_below_total(self):
"""Each branch contribution should be <= total spectrum.
Since total = sum(intensity_i * branch_i) * universal,
and all factors are positive, each branch's contribution
should not exceed the total.
"""
config = SpectrumConfig(
Z_parent=59,
Z_daughter=60,
A_number=144,
endpoint_MeV=2.996,
e_step_MeV=0.001,
branches=[
BranchConfig(
endpoint_MeV=2.996, transition_type="F1", intensity=0.1469
),
BranchConfig(
endpoint_MeV=0.341, transition_type="F1", intensity=0.1220
),
BranchConfig(
endpoint_MeV=0.253, transition_type="F1", intensity=0.1217
),
],
)
spectrum = BetaSpectrum.from_config(config)
W, energies = spectrum.get_energy_grid(config)
total = spectrum(W)
bs = spectrum.get_branch_spectra(W)
components = spectrum.calculate_components(W)
# Compute each branch's contribution to the total
for i, branch in enumerate(config.branches):
universal_product = np.ones_like(W)
for name in [
"Fermi",
"Screening",
"FiniteSizeL0",
"ChargeDistributionU",
"Exchange",
]:
if name in components:
universal_product *= components[name]
contrib = universal_product * branch.intensity * bs[i]
# Branch contribution should not exceed total
assert np.all(
contrib <= total + 1e-10
), f"Branch {i} contribution exceeds total at some points"
# Branch contribution should be non-negative
assert np.all(contrib >= 0), f"Branch {i} contribution has negative values"
def test_branch_integral_proportional_to_intensity(self):
"""Integral of each branch contribution should be proportional to its intensity."""
config = SpectrumConfig(
Z_parent=59,
Z_daughter=60,
A_number=144,
endpoint_MeV=2.996,
e_step_MeV=0.001,
branches=[
BranchConfig(endpoint_MeV=2.996, transition_type="F1", intensity=0.5),
BranchConfig(endpoint_MeV=0.5, transition_type="F1", intensity=0.3),
BranchConfig(endpoint_MeV=0.3, transition_type="F1", intensity=0.2),
],
)
spectrum = BetaSpectrum.from_config(config)
W, energies = spectrum.get_energy_grid(config)
bs = spectrum.get_branch_spectra(W)
components = spectrum.calculate_components(W)
integrals = []
for i, branch in enumerate(config.branches):
universal_product = np.ones_like(W)
for name in [
"Fermi",
"Screening",
"FiniteSizeL0",
"ChargeDistributionU",
"Exchange",
]:
if name in components:
universal_product *= components[name]
contrib = universal_product * branch.intensity * bs[i]
integral = np.trapezoid(contrib, energies)
integrals.append(integral)
total_integral = np.trapezoid(spectrum(W), energies)
# Each branch integral should be proportional to its intensity
# (within numerical tolerance, since different endpoints give different
# phase space integrals)
for i, (integral, intensity) in enumerate(zip(integrals, config.branches)):
assert integral > 0, f"Branch {i} has zero integral"
# Sum of branch integrals should equal total
assert np.isclose(
sum(integrals), total_integral, rtol=1e-6
), "Sum of branch integrals should equal total integral"