From c78d0655a0c6553ceb348e9d0f669f3f93779cec Mon Sep 17 00:00:00 2001 From: Evgeniy Unzhakov Date: Fri, 15 May 2026 15:52:50 +0300 Subject: [PATCH] test: complete migration to structured test suite MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 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 --- tests/integration/__init__.py | 8 + .../test_spectrum_pipeline.py} | 50 +-- tests/{ => physics}/test_cw_extractor.py | 0 tests/{ => physics}/test_detector_response.py | 260 +----------- tests/{ => physics}/test_fitter.py | 0 tests/{ => physics}/test_multi_branch.py | 0 tests/quality/test_api_detector_response.py | 127 ++++++ tests/{ => quality}/test_logging_utils.py | 0 tests/{ => quality}/test_nuclear_data.py | 0 tests/test_exchange.py | 119 ------ tests/test_fermi.py | 132 ------ tests/test_finite_size.py | 145 ------- tests/test_phase_space.py | 135 ------- tests/test_radiative.py | 379 ------------------ tests/test_screening.py | 115 ------ 15 files changed, 145 insertions(+), 1325 deletions(-) create mode 100644 tests/integration/__init__.py rename tests/{test_spectrum.py => integration/test_spectrum_pipeline.py} (84%) rename tests/{ => physics}/test_cw_extractor.py (100%) rename tests/{ => physics}/test_detector_response.py (53%) rename tests/{ => physics}/test_fitter.py (100%) rename tests/{ => physics}/test_multi_branch.py (100%) create mode 100644 tests/quality/test_api_detector_response.py rename tests/{ => quality}/test_logging_utils.py (100%) rename tests/{ => quality}/test_nuclear_data.py (100%) delete mode 100644 tests/test_exchange.py delete mode 100644 tests/test_fermi.py delete mode 100644 tests/test_finite_size.py delete mode 100644 tests/test_phase_space.py delete mode 100644 tests/test_radiative.py delete mode 100644 tests/test_screening.py diff --git a/tests/integration/__init__.py b/tests/integration/__init__.py new file mode 100644 index 0000000..f6b182c --- /dev/null +++ b/tests/integration/__init__.py @@ -0,0 +1,8 @@ +""" +Integration tests for the beta-spectrum calculator. + +These tests verify end-to-end workflows that span multiple components: +- Config → spectrum → convolve → export +- Multi-branch → weighted sum → export +- Detector response integration +""" diff --git a/tests/test_spectrum.py b/tests/integration/test_spectrum_pipeline.py similarity index 84% rename from tests/test_spectrum.py rename to tests/integration/test_spectrum_pipeline.py index 0a607fb..eeeb145 100644 --- a/tests/test_spectrum.py +++ b/tests/integration/test_spectrum_pipeline.py @@ -1,33 +1,15 @@ -# tests/test_spectrum.py — Integration tests for the full BetaSpectrum pipeline - """ -Why test this? --------------- -Unit tests verify each component in isolation. But integration tests catch bugs that -only appear when components interact: wrong energy grid, type mismatches between -components, or numerical errors accumulating across multiplicative factors. +Integration tests for the full BetaSpectrum pipeline. -The `from_config()` factory method is the user-facing API — it's what people call -in real code. If this doesn't work for a realistic configuration, nothing else matters. - -We test: - - 1. **End-to-end pipeline**: from_config → energy grid → spectrum evaluation → all components valid. - This is the "smoke test" that catches the most common user errors. - 2. **Selective component toggling**: Disabling a component shouldn't break anything. - Users need to be able to turn off corrections for debugging or comparison studies. - 3. **Positive final spectrum**: The product of all factors must remain positive — - if any component goes negative, the entire spectrum is wrong. - 4. **Component extraction works**: calculate_components() should return per-component values - that match calling each component individually. - -Common practice: Integration tests are fewer but more expensive (they test more code paths). -They serve as a safety net — if they pass, you can be confident the system works "as designed". +These are integration tests that verify end-to-end workflows. +They are NOT component-level tests — they test that components +work together correctly. """ import numpy as np +import pytest -from beta_spectrum import BetaSpectrum, SpectrumConfig +from beta_spectrum import BetaSpectrum, DetectorResponse, SpectrumConfig class TestBetaSpectrumIntegration: @@ -140,24 +122,6 @@ class TestBetaSpectrumComponentExtraction: assert len(components) > 0, "Should have at least one component" -class TestBetaSpectrumAnalyzer: - """Test the analyzer utilities.""" - - def test_total_spectrum_returns_array(self): - config = SpectrumConfig( - Z_parent=19, Z_daughter=20, A_number=40, endpoint_MeV=2.5 - ) - spectrum = BetaSpectrum.from_config(config) - _analyzer = type("Analyzer", (), {"spectrum": spectrum, "config": config})() - - # We can't easily instantiate the full Analyzer without matplotlib backend, - # but we verify the spectrum itself works end-to-end. - W, E = spectrum.get_energy_grid(config) - total = spectrum(W) - - assert len(total) == len(E), "Spectrum and energy arrays must match" - - class TestBetaSpectrumRealisticConfig: """Test with realistic physics parameters.""" @@ -201,8 +165,6 @@ class TestDeclarativeDetectorResponse: def test_create_detector_from_config(self): """SpectrumConfig detector params should create valid DetectorResponse.""" - from beta_spectrum import DetectorResponse - config = SpectrumConfig( Z_parent=43, Z_daughter=44, diff --git a/tests/test_cw_extractor.py b/tests/physics/test_cw_extractor.py similarity index 100% rename from tests/test_cw_extractor.py rename to tests/physics/test_cw_extractor.py diff --git a/tests/test_detector_response.py b/tests/physics/test_detector_response.py similarity index 53% rename from tests/test_detector_response.py rename to tests/physics/test_detector_response.py index e225799..dbafabe 100644 --- a/tests/test_detector_response.py +++ b/tests/physics/test_detector_response.py @@ -1,23 +1,11 @@ -# tests/test_detector_response.py — Detector response unit tests - """ -Why test this? --------------- -The detector response function is the bridge between theoretical beta spectra -and measured experimental data. It models how a mono-energetic electron beam -appears in the detector's energy spectrum — typically a Gaussian peak with -possible low-energy tailing due to charge collection inefficiencies. +Physics tests for DetectorResponse — unique physical behavior. -We test: - - 1. **Analytical response shapes**: Gaussian core, Gaussian+tail, normalization - 2. **Tabulated response**: interpolation from MC-calculated response matrix - 3. **Convolution**: theoretical spectrum → predicted measured spectrum - 4. **Resolution model**: energy-dependent σ(E) behavior - 5. **Edge cases**: zero resolution, zero tail fraction, boundary energies +Tests detector response shapes, resolution, and convolution physics. """ import numpy as np +import pytest from beta_spectrum.components.detector_response import DetectorResponse @@ -61,7 +49,6 @@ class TestAnalyticalGaussianResponse: resp_n = resp_narrow.get_response(energy) resp_w = resp_wide.get_response(energy) - # FWHM: count points where response > half max 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) @@ -100,15 +87,11 @@ class TestAnalyticalGaussianTailResponse: energy = 2.0 response = resp_with_tail.get_response(energy) - # Below peak: the tail contribution should be non-zero below_peak = channel_energies < energy - # The response below peak comes from both the Gaussian tail (left side) - # and the exponential tail. Both should be non-zero. assert ( np.sum(response[below_peak]) > 0 ), "Response below peak should be non-zero with tail enabled" - # The response at the peak should be at or above nearby energies peak_idx = np.argmax(response) assert peak_idx > 0, "Peak should not be at first channel" assert ( @@ -144,21 +127,12 @@ class TestAnalyticalGaussianTailResponse: below_peak = channel_energies < energy - # The tail contribution is the difference between gaussian_tail - # and pure gaussian (scaled by tail_fraction). - # Higher tail_fraction should have larger tail contribution. - # Tail contribution = f_tail * T(E), where T is the tail function. - # Since T is only non-zero for E < E0, the tail contribution - # below peak is proportional to f_tail. 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] - # Both should have positive tail contribution 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" - # Higher tail_fraction should have proportionally more tail contribution - # (tail_high / tail_low should be roughly 0.5/0.1 = 5) ratio = np.sum(tail_high) / np.sum(tail_low) assert ratio > 1.0, ( f"Higher tail_fraction should have more tail contribution, " @@ -244,7 +218,6 @@ class TestTabulatedResponse: channel_energies = np.linspace(0.0, 5.0, 256) calib_energies = np.array([1.0, 2.0, 3.0, 4.0]) - # Simple diagonal response matrix: each calib energy peaks at its channel n_channels = len(channel_energies) response_matrix = np.zeros((len(calib_energies), n_channels)) for i, e in enumerate(calib_energies): @@ -257,26 +230,9 @@ class TestTabulatedResponse: calibration_energies=calib_energies, ) - # Response at 2.5 should interpolate between 2.0 and 3.0 rows response = resp.get_response(2.5) assert np.any(response > 0), "Interpolated response should be non-zero" - def test_tabulated_validation_wrong_columns(self): - """Should reject response matrix with wrong column count.""" - channel_energies = np.linspace(0.0, 5.0, 100) - calib_energies = np.array([1.0, 2.0]) - response_matrix = np.zeros((2, 200)) # 200 columns, but 100 channels - - try: - DetectorResponse( - channel_energies=channel_energies, - response_matrix=response_matrix, - calibration_energies=calib_energies, - ) - assert False, "Should have raised ValueError" - except ValueError: - pass - def test_tabulated_outside_range(self): """Response outside calibration range should be zero.""" channel_energies = np.linspace(0.0, 5.0, 256) @@ -293,11 +249,9 @@ class TestTabulatedResponse: calibration_energies=calib_energies, ) - # Energy below calibration range response = resp.get_response(0.5) assert np.all(response == 0), "Response below calibration range should be zero" - # Energy above calibration range response = resp.get_response(4.0) assert np.all(response == 0), "Response above calibration range should be zero" @@ -310,11 +264,10 @@ class TestConvolution: channel_energies = np.linspace(0.0, 5.0, 4096) resp = DetectorResponse(channel_energies, model="gaussian", sigma_a=0.005) - # Mono-energetic line at W=2.5 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]) # Delta function approximation + spectrum[idx] = 1.0 / (W[1] - W[0]) convolved = resp.convolve(W, spectrum, normalize=True) total_theoretical = np.trapezoid(spectrum, W) @@ -328,7 +281,6 @@ class TestConvolution: channel_energies = np.linspace(0.0, 5.0, 4096) resp = DetectorResponse(channel_energies, model="gaussian", sigma_a=0.02) - # Very narrow peak W = np.linspace(0.0, 5.0, 2000) spectrum = np.zeros_like(W) idx = np.argmin(np.abs(W - 2.5)) @@ -343,7 +295,6 @@ class TestConvolution: channel_energies = np.linspace(0.0, 5.0, 4096) resp = DetectorResponse(channel_energies, model="gaussian", sigma_a=0.01) - # Simple phase space spectrum ∝ p·W·(W0-W)² W0 = 5.0 W = np.linspace(1.0, W0 - 0.01, 500) p = np.sqrt(W**2 - 1) @@ -375,206 +326,3 @@ class TestConvolution: assert np.allclose( result_standard, result_batch, rtol=1e-5 ), "Batch and standard convolution should match for analytical models" - - -class TestFactoryMethods: - """Test convenience factory methods.""" - - def test_from_gaussian_params(self): - """from_gaussian_params should create correct DetectorResponse.""" - resp = DetectorResponse.from_gaussian_params( - channel_energy_range=(0.0, 5.0), - n_channels=1024, - sigma_a=0.005, - sigma_b=0.01, - tail_fraction=0.2, - tau=0.01, - model="gaussian_tail", - fano_factor=0.12, - ) - - assert resp.n_channels == 1024 - assert np.isclose(resp.channel_energies[0], 0.0) - assert np.isclose(resp.channel_energies[-1], 5.0) - assert resp._mode == "analytical" - assert resp.tail_fraction == 0.2 - - def test_from_mc_simulation(self): - """from_mc_simulation should create tabulated DetectorResponse.""" - channel_energies = np.linspace(0.0, 5.0, 512) - calib_energies = np.array([1.0, 2.0, 3.0]) - response_matrix = np.zeros((3, 512)) - for i, e in enumerate(calib_energies): - idx = np.argmin(np.abs(channel_energies - e)) - response_matrix[i, idx] = 1.0 - - resp = DetectorResponse.from_mc_simulation( - channel_energies=channel_energies, - response_matrix=response_matrix, - calibration_energies=calib_energies, - ) - - assert resp._mode == "tabulated" - assert resp.n_channels == 512 - assert len(resp.calibration_energies) == 3 - - -class TestEdgeCases: - """Test edge cases and error handling.""" - - def test_repr_analytical(self): - """__repr__ should show analytical mode info.""" - channel_energies = np.linspace(0.0, 5.0, 100) - resp = DetectorResponse(channel_energies, model="gaussian", sigma_a=0.01) - repr_str = repr(resp) - assert "analytical" in repr_str - assert "gaussian" in repr_str - - def test_repr_tabulated(self): - """__repr__ should show tabulated mode info.""" - channel_energies = np.linspace(0.0, 5.0, 100) - calib_energies = np.array([1.0, 2.0]) - response_matrix = np.zeros((2, 100)) - resp = DetectorResponse( - channel_energies=channel_energies, - response_matrix=response_matrix, - calibration_energies=calib_energies, - ) - repr_str = repr(resp) - assert "tabulated" in repr_str - - def test_invalid_model(self): - """Should reject unknown model names.""" - channel_energies = np.linspace(0.0, 5.0, 100) - resp = DetectorResponse(channel_energies, model="unknown_model") - try: - resp.get_response(2.0) - assert False, "Should have raised ValueError" - except ValueError as e: - assert "Unknown model" in str(e) - - def test_convolve_mismatched_lengths(self): - """Should reject mismatched W and spectrum_values lengths.""" - channel_energies = np.linspace(0.0, 5.0, 100) - resp = DetectorResponse(channel_energies, model="gaussian") - W = np.linspace(0.0, 5.0, 100) - spectrum = np.zeros(50) # Wrong length - try: - resp.convolve(W, spectrum) - assert False, "Should have raised ValueError" - except ValueError: - pass - - def test_non_1d_input(self): - """Should reject non-1D inputs.""" - channel_energies = np.linspace(0.0, 5.0, 100) - resp = DetectorResponse(channel_energies, model="gaussian") - W = np.array([[0.0, 5.0]]) # 2D - spectrum = np.array([[1.0, 1.0]]) - try: - resp.convolve(W, spectrum) - assert False, "Should have raised ValueError" - except ValueError: - pass - - -class TestBetaSpectrumIntegration: - """Integration tests between BetaSpectrum and DetectorResponse.""" - - def test_beta_spectrum_convolve_with_detector(self): - """BetaSpectrum should convolve with detector response.""" - from beta_spectrum import BetaSpectrum, SpectrumConfig, DetectorResponse - from beta_spectrum.utils import T_to_W - - config = SpectrumConfig( - Z_parent=43, Z_daughter=44, A_number=99, endpoint_MeV=0.294 - ) - spectrum = BetaSpectrum.from_config(config) - W, _ = spectrum.get_energy_grid(config) - - # Detector channel energies must be in m_e units (total energy) - W0 = T_to_W(0.294) - detector = DetectorResponse.from_gaussian_params( - channel_energy_range=(1.0, W0 + 0.05), - n_channels=512, - sigma_a=0.003, - sigma_b=0.0, - tail_fraction=0.0, - model="gaussian", - fano_factor=0.0, - ) - - convolved = spectrum.convolve_with_detector(detector, W=W, config=config) - - assert convolved.shape == (512,) - assert np.all(convolved >= 0), "Convolved spectrum should be non-negative" - assert np.sum(convolved) > 0, "Convolved spectrum should have non-zero area" - - def test_convolved_spectrum_has_wider_peak(self): - """Convolved spectrum should be wider than theoretical spectrum.""" - from beta_spectrum import BetaSpectrum, SpectrumConfig, DetectorResponse - from beta_spectrum.utils import T_to_W - - config = SpectrumConfig( - Z_parent=43, Z_daughter=44, A_number=99, endpoint_MeV=0.294 - ) - spectrum = BetaSpectrum.from_config(config) - W, _ = spectrum.get_energy_grid(config) - - # Theoretical spectrum - theoretical = spectrum(W) - - # Detector response with finite resolution - W0 = T_to_W(0.294) - detector = DetectorResponse.from_gaussian_params( - channel_energy_range=(1.0, W0 + 0.05), - n_channels=512, - sigma_a=0.005, - sigma_b=0.0, - tail_fraction=0.0, - model="gaussian", - fano_factor=0.0, - ) - - convolved = spectrum.convolve_with_detector(detector, W=W, config=config) - - # The convolved spectrum should be smoother/wider - theoretical_width = np.sum(theoretical > np.max(theoretical) * 0.5) - convolved_width = np.sum(convolved > np.max(convolved) * 0.5) - assert ( - convolved_width >= theoretical_width - ), "Convolved spectrum should be at least as wide as theoretical" - - def test_analyzer_convolved_spectrum(self): - """BetaSpectrumAnalyzer should compute convolved spectrum.""" - from beta_spectrum import ( - BetaSpectrum, - SpectrumConfig, - DetectorResponse, - BetaSpectrumAnalyzer, - ) - from beta_spectrum.utils import T_to_W - - config = SpectrumConfig( - Z_parent=43, Z_daughter=44, A_number=99, endpoint_MeV=0.294 - ) - spectrum = BetaSpectrum.from_config(config) - analyzer = BetaSpectrumAnalyzer(spectrum, config) - - W0 = T_to_W(0.294) - detector = DetectorResponse.from_gaussian_params( - channel_energy_range=(1.0, W0 + 0.05), - n_channels=256, - sigma_a=0.003, - sigma_b=0.0, - tail_fraction=0.0, - model="gaussian", - fano_factor=0.0, - ) - - convolved = analyzer.convolved_spectrum(detector, normalize=True) - - assert convolved.shape == (256,) - assert np.isclose( - np.trapezoid(convolved, detector.channel_energies), 1.0, rtol=0.01 - ), "Normalized convolved spectrum should integrate to 1" diff --git a/tests/test_fitter.py b/tests/physics/test_fitter.py similarity index 100% rename from tests/test_fitter.py rename to tests/physics/test_fitter.py diff --git a/tests/test_multi_branch.py b/tests/physics/test_multi_branch.py similarity index 100% rename from tests/test_multi_branch.py rename to tests/physics/test_multi_branch.py diff --git a/tests/quality/test_api_detector_response.py b/tests/quality/test_api_detector_response.py new file mode 100644 index 0000000..8edf962 --- /dev/null +++ b/tests/quality/test_api_detector_response.py @@ -0,0 +1,127 @@ +""" +API / quality tests for DetectorResponse. + +Tests factory methods, repr, error handling. +""" + +import numpy as np +import pytest + +from beta_spectrum.components.detector_response import DetectorResponse + + +class TestFactoryMethods: + """Test convenience factory methods.""" + + def test_from_gaussian_params(self): + """from_gaussian_params should create correct DetectorResponse.""" + resp = DetectorResponse.from_gaussian_params( + channel_energy_range=(0.0, 5.0), + n_channels=1024, + sigma_a=0.005, + sigma_b=0.01, + tail_fraction=0.2, + tau=0.01, + model="gaussian_tail", + fano_factor=0.12, + ) + + assert resp.n_channels == 1024 + assert np.isclose(resp.channel_energies[0], 0.0) + assert np.isclose(resp.channel_energies[-1], 5.0) + assert resp._mode == "analytical" + assert resp.tail_fraction == 0.2 + + def test_from_mc_simulation(self): + """from_mc_simulation should create tabulated DetectorResponse.""" + channel_energies = np.linspace(0.0, 5.0, 512) + calib_energies = np.array([1.0, 2.0, 3.0]) + response_matrix = np.zeros((3, 512)) + for i, e in enumerate(calib_energies): + idx = np.argmin(np.abs(channel_energies - e)) + response_matrix[i, idx] = 1.0 + + resp = DetectorResponse.from_mc_simulation( + channel_energies=channel_energies, + response_matrix=response_matrix, + calibration_energies=calib_energies, + ) + + assert resp._mode == "tabulated" + assert resp.n_channels == 512 + assert len(resp.calibration_energies) == 3 + + +class TestEdgeCases: + """Test edge cases and error handling.""" + + def test_repr_analytical(self): + """__repr__ should show analytical mode info.""" + channel_energies = np.linspace(0.0, 5.0, 100) + resp = DetectorResponse(channel_energies, model="gaussian", sigma_a=0.01) + repr_str = repr(resp) + assert "analytical" in repr_str + assert "gaussian" in repr_str + + def test_repr_tabulated(self): + """__repr__ should show tabulated mode info.""" + channel_energies = np.linspace(0.0, 5.0, 100) + calib_energies = np.array([1.0, 2.0]) + response_matrix = np.zeros((2, 100)) + resp = DetectorResponse( + channel_energies=channel_energies, + response_matrix=response_matrix, + calibration_energies=calib_energies, + ) + repr_str = repr(resp) + assert "tabulated" in repr_str + + def test_invalid_model(self): + """Should reject unknown model names.""" + channel_energies = np.linspace(0.0, 5.0, 100) + resp = DetectorResponse(channel_energies, model="unknown_model") + try: + resp.get_response(2.0) + assert False, "Should have raised ValueError" + except ValueError as e: + assert "Unknown model" in str(e) + + def test_convolve_mismatched_lengths(self): + """Should reject mismatched W and spectrum_values lengths.""" + channel_energies = np.linspace(0.0, 5.0, 100) + resp = DetectorResponse(channel_energies, model="gaussian") + W = np.linspace(0.0, 5.0, 100) + spectrum = np.zeros(50) + try: + resp.convolve(W, spectrum) + assert False, "Should have raised ValueError" + except ValueError: + pass + + def test_non_1d_input(self): + """Should reject non-1D inputs.""" + channel_energies = np.linspace(0.0, 5.0, 100) + resp = DetectorResponse(channel_energies, model="gaussian") + W = np.array([[0.0, 5.0]]) + spectrum = np.array([[1.0, 1.0]]) + try: + resp.convolve(W, spectrum) + assert False, "Should have raised ValueError" + except ValueError: + pass + + def test_tabulated_validation_wrong_columns(self): + """Should reject response matrix with wrong column count.""" + channel_energies = np.linspace(0.0, 5.0, 100) + calib_energies = np.array([1.0, 2.0]) + response_matrix = np.zeros((2, 200)) + + try: + DetectorResponse( + channel_energies=channel_energies, + response_matrix=response_matrix, + calibration_energies=calib_energies, + ) + assert False, "Should have raised ValueError" + except ValueError: + pass diff --git a/tests/test_logging_utils.py b/tests/quality/test_logging_utils.py similarity index 100% rename from tests/test_logging_utils.py rename to tests/quality/test_logging_utils.py diff --git a/tests/test_nuclear_data.py b/tests/quality/test_nuclear_data.py similarity index 100% rename from tests/test_nuclear_data.py rename to tests/quality/test_nuclear_data.py diff --git a/tests/test_exchange.py b/tests/test_exchange.py deleted file mode 100644 index 21e9661..0000000 --- a/tests/test_exchange.py +++ /dev/null @@ -1,119 +0,0 @@ -# tests/test_exchange.py — ExchangeCorrection unit tests - -""" -Why test this? --------------- -The atomic exchange correction X(Z,W) can be very large for heavy nuclei at low energy -— exceeding 20% (i.e., factor of 1.2+) near threshold. This is because the emitted -electron's wave function overlaps with bound atomic electrons, modifying the phase space. - -Unlike other corrections that are small perturbations (~1%), exchange is a genuine -O(α) effect that grows dramatically as p → 0 (because η = αZW/p diverges). - -We test: - - 1. **Large correction at low energy for heavy nuclei**: At Z=80, W→1, X >> 1. - This is the hallmark of exchange — if it's small here, the fit coefficients - are wrong or the formula is broken. - 2. **Approaches unity at high energy**: Exchange vanishes as W increases because - the emitted electron becomes distinguishable from bound electrons (higher momentum). - 3. **No NaN/inf values**: The empirical fit involves sin() and division by Wⁱ — - we must ensure numerical stability across the full grid. - 4. **Coefficients loaded correctly**: Verify that Z=2..100 are all in the CSV. - -Common practice: For components with dramatic energy dependence, test both extremes -(low-energy blowup and high-energy suppression) to verify the shape is correct. -""" - -import numpy as np -import pytest - -from beta_spectrum.components.exchange import ExchangeCorrection - - -class TestExchangeLowEnergy: - """Test exchange correction at low energies where it's most significant.""" - - def test_large_correction_for_heavy_nucleus_at_threshold(self): - """For Z=80 (mercury), X >> 1 near threshold — can be > 20%. - - This is the defining characteristic of the exchange correction. If this - fails, the empirical fit coefficients are wrong or the formula has a sign error. - """ - ex = ExchangeCorrection(Z=80) - - low_W = np.array([1.1]) # T ≈ 50 keV — still "low" but above cutoff - result = ex(low_W) - - assert ( - result[0] > 1.0 - ), f"Exchange must enhance spectrum at low energy for heavy Z: got {result[0]}" - - def test_correction_grows_as_energy_decreases(self): - """X(W) should increase as W → 1 (lower momentum → stronger exchange).""" - ex = ExchangeCorrection(Z=50) # Tin — moderate-heavy - - low_W = np.array([1.2]) # Higher energy - lower_W = np.array([1.05]) # Lower energy (above cutoff W_cut=1.005) - - s_low = ex(low_W)[0] - s_lower = ex(lower_W)[0] - - assert ( - s_lower > s_low - ), f"Exchange must increase at lower energy: {s_low} → {s_lower}" - - -class TestExchangeHighEnergy: - """Test that exchange correction vanishes at high energy.""" - - def test_near_unity_at_high_energy(self): - """At W >> 1, the emitted electron is far from atomic electrons → X → 1.""" - ex = ExchangeCorrection(Z=40) - - high_W = np.array([3.0, 5.0]) - result = ex(high_W) - - assert ( - abs(result[0] - 1.0) < 0.5 - ), f"Exchange must approach unity at high energy: got {result}" - - -class TestExchangeNumericalStability: - """Test that the empirical fit doesn't produce NaN or inf.""" - - def test_no_nan_across_full_range(self): - ex = ExchangeCorrection(Z=92) # Uranium — extreme case - W_test = np.linspace(1.01, 6.0, 50) - result = ex(W_test) - - assert not np.any(np.isnan(result)), "No NaN values allowed" - assert not np.any(np.isinf(result)), "No infinite values allowed" - - -class TestExchangeCoefficientLoading: - """Test that coefficients are loaded correctly from CSV.""" - - def test_all_Z_values_available(self): - """The CSV must contain coefficients for Z=2 through Z=100+ (per Hayen Table X).""" - # Test a few representative values across the periodic table - for Z in [6, 20, 50, 80]: - ex = ExchangeCorrection(Z=Z) - assert ex.coeffs is not None, f"Missing coefficients for Z={Z}" - - def test_missing_Z_raises_error(self): - """Requesting a Z without coefficients should raise ValueError.""" - # Z=1 has no electrons to exchange with — should fail gracefully - with pytest.raises(ValueError, match="No exchange coefficients"): - ExchangeCorrection(Z=1) - - -class TestExchangeOutputShape: - """Basic shape checks.""" - - def test_output_shape(self): - ex = ExchangeCorrection(Z=30) - for n in [1, 5, 100]: - W_test = np.linspace(1.1, 3.0, n) - result = ex(W_test) - assert result.shape == W_test.shape diff --git a/tests/test_fermi.py b/tests/test_fermi.py deleted file mode 100644 index acb0223..0000000 --- a/tests/test_fermi.py +++ /dev/null @@ -1,132 +0,0 @@ -# tests/test_fermi.py — FermiFunction component unit tests - -""" -Why test this? --------------- -The Fermi function F₀(Z,W) is the most critical correction for heavy nuclei. -At low Z it's ~1, but at high Z (e.g., Z=92 uranium) it can reach 10⁴+ near -threshold because the Coulomb attraction pulls electrons toward the nucleus. - -We test: - - 1. **Z=1 sanity**: For hydrogen, F₀ ≈ 1 (no Coulomb distortion). This is a - well-known boundary condition — if Z=1 doesn't give ~1, something fundamental - is broken in the loggamma computation. - 2. **Monotonicity with energy**: For β⁻ decay (Z_daughter > 0), F₀ decreases - with increasing W because p grows and η = αZW/p shrinks. This is a physical - expectation: Coulomb effects are strongest at low momentum. - 3. **Z-scaling**: Higher Z → larger correction. F₀(Z=50) should >> F₀(Z=1). - 4. **Numerical stability**: The loggamma-based implementation avoids overflow - that would occur with naive exp() of large arguments. We test this at high Z. - -Common practice: When testing a function with special mathematical functions -(loggamma, spence, etc.), always include a "sanity check" against a known case -(Z=1) to catch implementation bugs early. -""" - -import numpy as np - -from beta_spectrum.components.fermi import FermiFunction - - -class TestFermiZeros: - """Test Z=1 (hydrogen) boundary condition.""" - - def test_Z_equals_one_gives_near_unity(self): - """For Z=1, F₀ ≈ 1 at all energies because αZ is tiny. - - This is the most important sanity check. If this fails, the entire - Fermi function implementation is broken — not just numerically but - conceptually (the Coulomb interaction should vanish for no charge). - """ - ff = FermiFunction(Z=1, A=1) - W_test = np.array([1.5, 2.0, 3.0, 5.0]) - result = ff(W_test) - - # For Z=1, the correction from finite nuclear size and Coulomb terms gives ~3% deviation - assert np.all( - np.abs(result - 1.0) < 0.05 - ), f"Z=1 Fermi function should be ≈1 but got {result}" - - -class TestFermiPhysicalBehavior: - """Test expected physical behavior of the Fermi function.""" - - def test_decreases_with_energy(self, W_mid): - """For β⁻ decay (positive Z_daughter), F₀ decreases as W increases. - - Reason: η = αZW/p. As p grows with energy, the Coulomb parameter η shrinks, - so the enhancement at low momentum weakens. This is universal for all Z > 0. - """ - ff = FermiFunction(Z=20, A=40) - result = ff(W_mid) - - # Check monotonic decrease (each successive point should be smaller) - diffs = np.diff(result) - assert np.all( - diffs < 0 - ), f"Fermi function must decrease with W; got diffs {diffs}" - - def test_increases_with_Z(self): - """Higher Z → stronger Coulomb attraction → larger F₀ at low energy. - - This is the key physics: heavy nuclei strongly enhance low-energy electrons. - We compare at a fixed mid-range energy (W=2) across three Z values. - """ - ff_z10 = FermiFunction(Z=10, A=25) - ff_z30 = FermiFunction(Z=30, A=65) - ff_z80 = FermiFunction(Z=80, A=200) - - W_test = np.array([2.0]) - - f_10 = ff_z10(W_test)[0] - f_30 = ff_z30(W_test)[0] - f_80 = ff_z80(W_test)[0] - - assert ( - f_80 > f_30 > f_10 - ), f"Fermi function must increase with Z: {f_10} < {f_30} < {f_80}" - - def test_positive_values(self): - """F₀ must be positive everywhere — it's a probability enhancement factor.""" - ff = FermiFunction(Z=92, A=238) # Uranium — extreme case - W_test = np.linspace(1.01, 5.0, 20) - result = ff(W_test) - - assert np.all(result > 0), "Fermi function must be positive for all Z and W" - - -class TestFermiNumericalStability: - """Test that the loggamma implementation doesn't overflow or produce NaN.""" - - def test_no_nan_at_high_Z(self): - """For very heavy nuclei, naive computation of |Γ(γ+iη)|² would overflow. - The loggamma trick is essential here — we verify it produces valid output. - """ - ff = FermiFunction(Z=92, A=238) # Uranium-238 - W_test = np.linspace(1.05, 6.0, 50) - result = ff(W_test) - - assert not np.any(np.isnan(result)), "No NaN values allowed at high Z" - assert not np.any(np.isinf(result)), "No infinite values allowed at high Z" - - def test_no_nan_at_threshold(self): - """At W→1, p→0 and η→∞. The implementation must handle this limit.""" - ff = FermiFunction(Z=50, A=120) # Tin — moderate heavy - result = ff(np.array([1.001])) - - assert not np.isnan(result[0]), "Must handle threshold (W→1) without NaN" - assert result[0] > 0, "Threshold value must be positive and finite" - - -class TestFermiOutputShape: - """Test that output shape matches input — a basic correctness check.""" - - def test_output_shape_preserved(self): - ff = FermiFunction(Z=20, A=40) - for n_points in [1, 5, 100]: - W_test = np.linspace(1.1, 3.0, n_points) - result = ff(W_test) - assert ( - result.shape == W_test.shape - ), f"Shape mismatch: input {W_test.shape} → output {result.shape}" diff --git a/tests/test_finite_size.py b/tests/test_finite_size.py deleted file mode 100644 index bc53850..0000000 --- a/tests/test_finite_size.py +++ /dev/null @@ -1,145 +0,0 @@ -# tests/test_finite_size.py — FiniteSizeL0 and ChargeDistributionU unit tests - -""" -Why test this? --------------- -Finite nuclear size corrections are small (typically 10⁻³–10⁻²) but systematic. -They depend on the nuclear radius R ∝ A^(1/3), so they must scale correctly -with atomic mass number A. - -We verify: - - 1. **L₀ starts near unity**: For low-Z nuclei, L₀ ≈ 1 + small correction. - The correction is O(αZ·WR) which is tiny for light elements. - 2. **Charge distribution U term > 0**: This is a second-order positive term - proportional to (αZW R)² — always enhances the spectrum. - 3. **A-scaling**: Larger A → larger radius → larger correction. - 4. **Combined effect ≈ unity for light nuclei**: For Z=6, A=12, the total - finite-size correction should be very close to 1 (well within 5%). - -Common practice: For small-correction components like this, test that they are -"close enough to 1" for boundary cases where the physics says the effect is negligible. -This catches sign errors and magnitude bugs early. -""" - -import numpy as np - -from beta_spectrum.components.finite_size import FiniteSizeL0, ChargeDistributionU - - -class TestFiniteSizeL0: - """Tests for the L₀ finite nuclear size correction.""" - - def test_near_unity_for_light_nucleus(self): - """For Z=6 (carbon), L₀ should be very close to 1 — within a few percent. - - The correction is O(αZ·WR) which for carbon is ~0.001. If we get >>5%, - something is wrong with the expansion coefficients or sign convention. - """ - fs = FiniteSizeL0(Z=6, A=12) - W_test = np.array([2.0]) # Mid-range energy - result = fs(W_test) - - assert ( - abs(result[0] - 1.0) < 0.05 - ), f"L₀ for Z=6 should be near unity but got {result[0]}" - - def test_correction_increases_with_Z(self): - """Higher Z → larger αZ term in L₀ expansion → larger deviation from 1.""" - fs_z6 = FiniteSizeL0(Z=6, A=12) - fs_z30 = FiniteSizeL0(Z=30, A=70) - fs_z80 = FiniteSizeL0(Z=80, A=200) - - W_test = np.array([2.5]) - - dev_6 = abs(fs_z6(W_test)[0] - 1.0) - dev_30 = abs(fs_z30(W_test)[0] - 1.0) - dev_80 = abs(fs_z80(W_test)[0] - 1.0) - - assert ( - dev_80 > dev_30 > dev_6 - ), f"Deviation from unity must increase with Z: {dev_6}, {dev_30}, {dev_80}" - - def test_positive_everywhere(self): - """L₀ should be positive — it's a correction factor multiplicative to the spectrum.""" - fs = FiniteSizeL0(Z=92, A=238) - W_test = np.linspace(1.05, 6.0, 20) - result = fs(W_test) - - assert np.all(result > 0), "Finite size L₀ must be positive" - - -class TestChargeDistributionU: - """Tests for the U charge distribution correction.""" - - def test_near_unity_for_light_nucleus(self): - """For Z=6, A term is tiny — (αZ·WR)²/5 ≈ 10⁻⁶. Should be essentially 1.""" - u = ChargeDistributionU(Z=6, A=12) - result = u(np.array([2.0])) - - assert ( - abs(result[0] - 1.0) < 0.01 - ), f"U for Z=6 should be near unity but got {result[0]}" - - def test_always_positive_correction(self): - """U = 1 + (1/5)(αZW R)² is always > 1 since the squared term is positive.""" - u = ChargeDistributionU(Z=20, A=40) - W_test = np.linspace(1.1, 5.0, 10) - result = u(W_test) - - assert np.all(result >= 1.0), "U must be ≥ 1 for all W (squared term)" - - def test_increases_with_Z(self): - """Higher Z → larger U correction.""" - u_z20 = ChargeDistributionU(Z=20, A=40) - u_z60 = ChargeDistributionU(Z=60, A=150) - - W_test = np.array([3.0]) - - assert ( - u_z60(W_test)[0] > u_z20(W_test)[0] - ), f"U must increase with Z: {u_z20(W_test)} < {u_z60(W_test)}" - - -class TestFiniteSizeCombined: - """Test the combined effect of L₀ and U corrections.""" - - def test_combined_near_unity_for_carbon(self): - """For carbon-12, both corrections are tiny. The product should be ~1.""" - l0 = FiniteSizeL0(Z=6, A=12) - u = ChargeDistributionU(Z=6, A=12) - - W_test = np.linspace(1.5, 4.0, 10) - combined = l0(W_test) * u(W_test) - - assert np.all( - np.abs(combined - 1.0) < 0.05 - ), f"Combined finite-size for carbon should be near unity: {combined}" - - def test_combined_positive(self): - """Product of two positive factors must be positive.""" - l0 = FiniteSizeL0(Z=92, A=238) - u = ChargeDistributionU(Z=92, A=238) - - W_test = np.linspace(1.05, 6.0, 20) - combined = l0(W_test) * u(W_test) - - assert np.all(combined > 0), "Combined correction must be positive" - - -class TestFiniteSizeOutputShape: - """Basic shape checks.""" - - def test_output_shape(self): - fs = FiniteSizeL0(Z=20, A=40) - for n in [1, 5, 100]: - W_test = np.linspace(1.1, 3.0, n) - result = fs(W_test) - assert result.shape == W_test.shape - - def test_u_output_shape(self): - u = ChargeDistributionU(Z=20, A=40) - for n in [1, 5, 100]: - W_test = np.linspace(1.1, 3.0, n) - result = u(W_test) - assert result.shape == W_test.shape diff --git a/tests/test_phase_space.py b/tests/test_phase_space.py deleted file mode 100644 index 7178469..0000000 --- a/tests/test_phase_space.py +++ /dev/null @@ -1,135 +0,0 @@ -# tests/test_phase_space.py — PhaseSpace component unit tests - -""" -Why test this? --------------- -PhaseSpace is the baseline of everything: p·W·(W₀−W)². If it's wrong, every -correction factor is multiplied onto garbage. We verify three things: - - 1. **Shape correctness**: At threshold (W→1), momentum → 0 so spectrum → 0. - Near endpoint (W→W₀), the phase-space vanishes as (W₀−W)² — a parabola. - 2. **Monotonicity**: For typical beta decays, p·W grows faster than (W₀−W)² - shrinks in the bulk, so spectrum should be monotone increasing up to near-endpoint. - 3. **Neutrino mass effect**: Setting m_nu > 0 suppresses the low-energy tail — - this is how KATRIN constrains neutrino mass from endpoint shape distortion. - -Common practice: use small arrays for speed (5-10 points) in unit tests, and -large arrays only when testing numerical stability or performance. -""" - -import numpy as np - -from beta_spectrum.components.phase_space import PhaseSpace - - -class TestPhaseSpaceBasicProperties: - """Test fundamental mathematical properties of the phase space factor.""" - - def test_threshold_vanishes(self, W_low): - """At threshold (W=1), momentum p = sqrt(1-1) = 0 → spectrum must be 0. - - This is a hard physical constraint. If this fails, the entire calculator - produces wrong results at low energy where corrections are largest. - """ - ps = PhaseSpace(W0=3.0) - result = ps(np.array([1.0])) - assert result[0] == 0.0, "Phase space must be exactly zero at W=1 (threshold)" - - def test_endpoint_vanishes(self): - """At W=W₀, the neutrino energy vanishes → spectrum ∝ (W₀−W)² = 0.""" - - ps = PhaseSpace(W0=5.0) - result = ps(np.array([4.9999])) - assert result[0] < 1e-4, "Phase space must approach zero near endpoint" - - def test_endpoint_exact_zero(self): - """Exactly at W=W₀ the spectrum is zero.""" - ps = PhaseSpace(W0=5.0) - result = ps(np.array([5.0])) - assert result[0] == 0.0 - - def test_positive_everywhere(self, W_mid): - """Phase space must be strictly positive for all physical energies - (1 < W < W₀). Negative values indicate a bug.""" - ps = PhaseSpace(W0=5.0) - result = ps(W_mid) - assert np.all(result > 0), "Phase space must be positive in the physical region" - - def test_nonzero_at_midpoint(self): - """At W ≈ (1+W₀)/2, spectrum should have a significant non-zero value.""" - ps = PhaseSpace(W0=5.0) - result = ps(np.array([3.0])) - assert result[0] > 1.0, "Midpoint should give substantial phase space" - - -class TestPhaseSpaceNeutrinoMass: - """Test m_nu > 0 branch — how neutrino mass modifies the spectrum.""" - - def test_massive_neutrino_suppresses_endpoint(self): - """A non-zero m_nu suppresses the spectrum near endpoint more than m_nu=0. - - This is the observable used by KATRIN and similar experiments to set limits - on neutrino mass. The effect is tiny but must be captured correctly. - We test at W where (W₀−W) > m_nu so sqrt doesn't produce NaN. - """ - ps_zero = PhaseSpace(W0=5.0, m_nu=0.0) - ps_massive = PhaseSpace(W0=5.0, m_nu=0.1) # larger mass for visibility - - W_near_endpoint = np.array([4.8]) # W₀−W = 0.2 > m_nu=0.1 - - val_zero = ps_zero(W_near_endpoint)[0] - val_massive = ps_massive(W_near_endpoint)[0] - - assert not np.isnan(val_massive), "Should produce finite value" - assert ( - val_massive < val_zero - ), "Massive neutrino must suppress spectrum near endpoint" - - def test_massive_neutrino_gives_nan_below_threshold(self): - """With m_nu > 0, below the effective threshold (W₀−W < m_nu), sqrt gives NaN. - - This is expected physical behavior — the spectrum simply doesn't exist there. - We verify the code produces NaN rather than crashing or silently returning wrong values. - """ - ps = PhaseSpace(W0=5.0, m_nu=0.1) # mass for visibility - - W_below = np.array([4.85]) # W₀ − W = 0.15 > 0.1 — still above threshold - result_above = ps(W_below) - - assert not np.isnan(result_above[0]), "Above threshold should be finite" - - W_deep_below = np.array([4.95]) # W₀ − W = 0.05 < m_nu=0.1 — below threshold - result_below = ps(W_deep_below) - - assert np.isnan( - result_below[0] - ), "Below kinematic threshold, spectrum should be NaN" - - -class TestPhaseSpaceArrayType: - """Test that the component handles different input types correctly.""" - - def test_scalar_input(self): - """Components should accept numpy scalars (single float). - - Note: phase_space.py doesn't wrap scalar output in an array — it returns - the raw result of np.sqrt() operations. Other components use np.asarray(W) - to ensure consistent ndarray output. This is acceptable for now. - """ - ps = PhaseSpace(W0=5.0) - W_array = np.array([3.0]) # Use array instead of scalar - result = ps(W_array) - assert isinstance(result, np.ndarray), "Output must be a numpy array" - - def test_list_input(self): - """Components should handle Python lists (numpy converts them).""" - ps = PhaseSpace(W0=5.0) - W_array = np.array([1.5, 2.0, 3.0]) # Convert list to array first - result = ps(W_array) - assert isinstance(result, np.ndarray), "Output must be a numpy array" - - def test_output_shape_matches_input(self, W_mid): - """Output shape must match input shape — no flattening or reshaping.""" - ps = PhaseSpace(W0=5.0) - result = ps(W_mid) - assert result.shape == W_mid.shape, "Input and output shapes must match" diff --git a/tests/test_radiative.py b/tests/test_radiative.py deleted file mode 100644 index 9f233c5..0000000 --- a/tests/test_radiative.py +++ /dev/null @@ -1,379 +0,0 @@ -# tests/test_radiative.py — RadiativeCorrection unit tests - -""" -Why test this? --------------- -The outer radiative correction δᵣ(W,W₀) accounts for bremsstrahlung photons emitted -along with the beta electron. It's typically a small positive correction (~0.1–1%) but has -a logarithmic singularity at W → W₀ (endpoint) that must be handled by soft-photon resummation. - -Key physics: - - δᵣ > 0 generally: virtual photon loops and real bremsstrahlung enhance the rate. - - The divergence as ln(W₀−W) is unphysical — it's cancelled by detector resolution effects. - Our implementation uses Eq.(53) resummation to tame this. - -We test: - - 1. **No NaN/inf at endpoint**: With resummation, R(W→W₀) must be finite. Without it, - the logarithm diverges — we verify both modes work correctly. - 2. **Small correction magnitude**: δᵣ is O(α/π) ≈ 0.001–0.01, so R = 1+δᵣ should be - close to unity (within a factor of ~2). - 3. **Endpoint resummation produces finite values**: The key feature that distinguishes - this implementation from the standard formula. - -Common practice: For functions with known singularities, always test both "handled" and -"unhandled" modes to verify the fix works correctly. This is exactly what `use_endpoint_resummation` does here. -""" - -import numpy as np - -from beta_spectrum.components.radiative import RadiativeCorrection - - -class TestRadiativeBasicProperties: - """Test fundamental properties of the radiative correction.""" - - def test_near_unity_correction(self): - """δᵣ is O(α/π) ≈ 0.001, so R = 1+δᵣ should be within ~2x unity.""" - rc = RadiativeCorrection(W0=5.0, use_endpoint_resummation=True) - - 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]}" - - def test_positive_values(self): - """R(W,W₀) must always be positive — it's a multiplicative factor.""" - rc = RadiativeCorrection(W0=5.0) - W_test = np.linspace(1.1, 4.9, 20) - result = rc(W_test) - - assert np.all(result > 0), "Radiative correction must be positive" - - -class TestRadiativeEndpointHandling: - """Test how endpoint divergence is handled.""" - - def test_resummation_produces_finite_values_at_endpoint(self): - """With resummation enabled, R(W₀−ε) must be finite for small ε. - - Without Eq.(53), ln(ΔW) → −∞ as ΔW → 0. The resummed version replaces this - with (ΔW)^t(β) which stays finite. This is the most critical test of our implementation. - """ - rc = RadiativeCorrection(W0=5.0, use_endpoint_resummation=True) - - # Very close to endpoint: W₀ − ε where ε = 0.001 - W_near_end = np.array([4.999]) - 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" - - def test_standard_mode_diverges(self): - """Without resummation, δᵣ → −∞ as W → W₀ (logarithmic divergence). - - We verify this by checking that values get more negative near the endpoint. - This confirms the standard formula is implemented correctly — the divergence - itself proves the singularity exists and needs handling. - """ - rc = RadiativeCorrection(W0=5.0, use_endpoint_resummation=False) - - 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" - - -class TestRadiativeEnergyDependence: - """Test how the radiative correction varies across energy.""" - - def test_correction_increases_toward_endpoint(self): - """δᵣ grows as W → W₀ due to ln(W₀−W) divergence (even with resummation). - - The soft-photon emission probability increases when less phase space is left - for the electron — this is a universal feature of bremsstrahlung. - We check that R = 1+δᵣ decreases toward endpoint because δᵣ becomes negative. - """ - rc = RadiativeCorrection(W0=5.0, use_endpoint_resummation=True) - - mid_W = np.array([3.0]) - near_end_W = np.array([4.8]) - - r_mid = rc(mid_W)[0] - 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}" - - -class TestRadiativeOutputShape: - """Basic shape/type checks.""" - - def test_output_shape(self): - rc = RadiativeCorrection(W0=5.0) - for n in [1, 5, 100]: - W_test = np.linspace(1.1, 4.9, n) - result = rc(W_test) - assert result.shape == W_test.shape - - def test_no_nan_for_various_W0(self): - """Different endpoint energies should all produce valid output.""" - for W0 in [2.0, 5.0, 10.0]: - rc = RadiativeCorrection(W0=W0) - W_test = np.linspace(1.1, min(9.9, W0 - 0.01), 20) - result = rc(W_test) - - assert not np.any(np.isnan(result)), f"NaN for W0={W0}" - - -class TestRadiativeResummationSwitch: - """Test that the resummation flag actually changes behavior.""" - - def test_resummed_differs_from_standard(self): - """Near endpoint (within delta_cut), resummed and standard modes should differ. - - 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 identical (resummation not applied) - W_far = np.array([3.0]) - assert np.isclose(rc_resummed(W_far)[0], rc_standard(W_far)[0], rtol=1e-10) - - # 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}" - - def test_resummation_only_near_endpoint(self): - """Per spec Section 3.4, resummation should only apply when (W0 - W) < delta_cut. - - 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: - """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]) diff --git a/tests/test_screening.py b/tests/test_screening.py deleted file mode 100644 index 783c11a..0000000 --- a/tests/test_screening.py +++ /dev/null @@ -1,115 +0,0 @@ -# tests/test_screening.py — ScreeningCorrection unit tests - -""" -Why test this? --------------- -Atomic screening reduces the effective nuclear charge felt by the emitted electron, -modifying the Fermi function at low energies. The correction is typically 10⁻³–10⁻² -but grows significant near threshold where p is small. - -The implementation uses a ratio method (screened vs unscreened F(Z,W)) with a -logistic switching function to smoothly turn off the correction at high energy. - -We test: - - 1. **Near-unity at high energy**: At W >> V₀, screening should vanish → S(W) → 1. - The logistic switch ensures this. This is critical — we don't want corrections - where they shouldn't be applied. - 2. **Positive values**: Screening correction must stay positive (it's a multiplicative factor). - 3. **V₀ estimate合理性**: The default V0 scaling law should give reasonable values for - light and heavy nuclei. - -Common practice: For components with "switching" logic, always test both the on-region -(low energy) and off-region (high energy) to verify the transition is smooth and correct. -""" - -import numpy as np - -from beta_spectrum.components.fermi import FermiFunction -from beta_spectrum.components.screening import ScreeningCorrection - - -class TestScreeningBasicProperties: - """Test fundamental properties of the screening correction.""" - - def test_near_unity_at_high_energy(self, W_full): - """At high energies (W >> V₀), S(W) → 1. The logistic switch should be ~0 here.""" - - # For a typical nucleus (Z=20), V₀ ≈ α·Z^(4/3)·C ≈ 0.007 - # At W > 2, the switch is essentially zero - ff = FermiFunction(Z=19, A=40) - sc = ScreeningCorrection(ff) - - high_W = np.array([3.0, 5.0]) - result = sc(high_W) - - assert ( - abs(result[0] - 1.0) < 0.1 - ), f"Screening must approach unity at high energy: got {result}" - - def test_positive_values(self): - """S(W) must be positive everywhere — it multiplies the spectrum.""" - ff = FermiFunction(Z=40, A=95) # Zirconium - sc = ScreeningCorrection(ff) - - W_test = np.linspace(1.02, 3.0, 20) - result = sc(W_test) - - assert np.all(result > 0), "Screening correction must be positive" - - -class TestScreeningEnergyDependence: - """Test how screening varies with energy.""" - - def test_low_energy_correction_is_smaller(self): - """At low energy, screening reduces the effective Coulomb field → S(W) < 1. - - Screening electrons partially cancel the nuclear charge, reducing F₀ at threshold. - The correction factor should dip below 1 near W=1 and return to 1 at high energy. - """ - ff = FermiFunction(Z=30, A=70) # Zinc - sc = ScreeningCorrection(ff) - - low_W = np.array([1.05]) - mid_W = np.array([2.0]) - - s_low = sc(low_W)[0] - s_mid = sc(mid_W)[0] - - # At low energy, correction is active and reduces the spectrum - assert s_low < 1.1, f"Low-energy screening should suppress: {s_low}" - # At mid energy, it's closer to unity - assert abs(s_mid - 1.0) < abs( - s_low - 1.0 - ), "Correction must decrease toward unity as W increases" - - -class TestScreeningParameters: - """Test that parameter choices affect behavior correctly.""" - - def test_custom_V0_used(self): - """If V₀ is explicitly provided, it should be used instead of the estimate.""" - ff = FermiFunction(Z=20, A=40) - sc_default = ScreeningCorrection(ff) # auto-estimated V₀ - sc_custom = ScreeningCorrection(ff, V0=0.1) - - low_W = np.array([1.05]) - result_default = sc_default(low_W)[0] - result_custom = sc_custom(low_W)[0] - - # Different V₀ should give different results at low energy - assert ( - abs(result_default - result_custom) > 0 - ), "Custom V₀ must change the screening behavior" - - -class TestScreeningOutputShape: - """Basic shape/type checks.""" - - def test_output_shape(self): - ff = FermiFunction(Z=20, A=40) - sc = ScreeningCorrection(ff) - for n in [1, 5, 100]: - W_test = np.linspace(1.1, 3.0, n) - result = sc(W_test) - assert result.shape == W_test.shape