diff --git a/tests/physics/test_core_spectrum.py b/tests/physics/test_core_spectrum.py new file mode 100644 index 0000000..652cdf7 --- /dev/null +++ b/tests/physics/test_core_spectrum.py @@ -0,0 +1,467 @@ +""" +Core-spectrum invariant tests. + +The basic spectrum is ALWAYS core = phase_space × fermi. Fermi is a core +spectrum component on par with phase space (its 0-order Coulomb correction), +never a higher-order correction. Finite size, charge distribution, screening, +exchange and radiative are higher-order corrections applied to the core. + +These tests pin the invariant at every layer: +- SpectrumConfig must not allow disabling the core components +- calculate_components always contains PhaseSpace and Fermi +- total spectrum equals universal × sum(intensity × branch) and every branch + physical spectrum contains the Fermi enhancement +""" + +from __future__ import annotations + +import inspect +from typing import Any, Dict, List, Tuple + +import numpy as np +import pandas as pd +import pytest + +from beta_spectrum import ( + BetaSpectrum, + BranchConfig, + SpectrumConfig, +) +from beta_spectrum.output import BetaSpectrumAnalyzer +from beta_spectrum.output.visualization import Visualizer, parse_spectrum_csv + +ME_MEV = 0.510998950 + + +def multi_branches() -> List[BranchConfig]: + return [ + BranchConfig(endpoint_MeV=0.294, transition_type="F2", intensity=0.98), + BranchConfig(endpoint_MeV=0.200, transition_type="F2U", intensity=0.15), + BranchConfig(endpoint_MeV=0.120, transition_type="A", intensity=0.005), + ] + + +def make_multi_config(**over: Any) -> SpectrumConfig: + return SpectrumConfig( + Z_parent=43, + Z_daughter=44, + A_number=99, + endpoint_MeV=0.294, + transition_type="F2", + e_step_MeV=0.001, + branches=multi_branches(), + **over, + ) + + +def make_single_config(**over: Any) -> SpectrumConfig: + return SpectrumConfig( + Z_parent=15, + Z_daughter=16, + A_number=32, + endpoint_MeV=1.7103, + transition_type="A", + e_step_MeV=0.005, + **over, + ) + + +def make_degenerate_config() -> SpectrumConfig: + """Multi-branch config where every endpoint lies below m_nu.""" + return SpectrumConfig( + Z_parent=43, + Z_daughter=44, + A_number=99, + endpoint_MeV=0.294, + e_step_MeV=0.001, + neutrino_mass_MeV=0.35, + branches=[ + BranchConfig(endpoint_MeV=0.294, transition_type="F2", intensity=0.98), + BranchConfig(endpoint_MeV=0.200, transition_type="F2U", intensity=0.02), + BranchConfig(endpoint_MeV=0.120, transition_type="A", intensity=0.005), + ], + ) + + +def branches_of(config: SpectrumConfig) -> List[BranchConfig]: + branches = config.branches + assert branches is not None + return branches + + +def export_and_read( + analyzer: BetaSpectrumAnalyzer, + tmp_path: Any, + name: str = "out.csv", + csv_verbosity: str = "full", +) -> Tuple[str, pd.DataFrame]: + path = str(tmp_path / name) + analyzer.export_to_csv(path, csv_verbosity=csv_verbosity) + df = pd.read_csv(path, sep=None, engine="python", comment="#") + return path, df + + +def header_lines_of(path: str) -> List[str]: + lines: List[str] = [] + with open(path, "r") as f: + for line in f: + if line.startswith("#"): + lines.append(line.rstrip("\n")) + else: + break + return lines + + +def universal_factor(components: Dict[str, np.ndarray], W: np.ndarray) -> np.ndarray: + """Product of all components shared by all branches. + + Universal components are the bare-named keys except the per-branch core + factors (PhaseSpace) and per-branch radiative corrections. + """ + uni = np.ones_like(np.asarray(W, dtype=float)) + for key, values in components.items(): + if "." in key: + continue + if key in ("PhaseSpace", "Radiative"): + continue + uni = uni * values + return uni + + +class TestCoreComponentsAreUnconditional: + """phase_space × fermi is the core: it can never be disabled.""" + + def test_config_has_no_core_toggles(self) -> None: + """use_fermi / use_phase_space must not exist as config fields. + + A dataclass __init__ rejects unknown keyword arguments with + TypeError, so the fields' absence is exactly the enforcement that + the core components can never be disabled. + """ + params = set(inspect.signature(SpectrumConfig).parameters) + assert "use_fermi" not in params + assert "use_phase_space" not in params + # Corrections remain configurable: + assert "use_radiative" in params + assert "use_screening" in params + + def test_single_branch_components_contain_core(self) -> None: + config = make_single_config() + spectrum = BetaSpectrum.from_config(config) + W, _ = spectrum.get_energy_grid(config) + components = spectrum.calculate_components(W) + assert "PhaseSpace" in components + assert "Fermi" in components + + def test_multi_branch_components_contain_core(self) -> None: + config = make_multi_config() + spectrum = BetaSpectrum.from_config(config) + W, _ = spectrum.get_energy_grid(config) + components = spectrum.calculate_components(W) + # Core phase space lives per-branch; Fermi is shared (bare key). + assert "Fermi" in components + for i in range(len(branches_of(config))): + assert f"branch_{i}.PhaseSpace" in components + + def test_all_corrections_disabled_still_has_core(self) -> None: + config = make_single_config( + use_screening=False, + use_finite_size=False, + use_charge_dist=False, + use_exchange=False, + use_radiative=False, + ) + spectrum = BetaSpectrum.from_config(config) + W, _ = spectrum.get_energy_grid(config) + components = spectrum.calculate_components(W) + assert set(components.keys()) == {"PhaseSpace", "Fermi"} + + def test_single_spectrum_equals_core_times_corrections(self) -> None: + config = make_single_config() + spectrum = BetaSpectrum.from_config(config) + W, _ = spectrum.get_energy_grid(config) + components = spectrum.calculate_components(W) + product = np.ones_like(W) + for values in components.values(): + product = product * values + assert np.allclose(spectrum(W), product, rtol=1e-10) + + +class TestMultiBranchCoreInvariant: + """total = universal × Σ(intensity × physical_branch), Fermi everywhere.""" + + def test_total_equals_universal_times_intensity_sum(self) -> None: + config = make_multi_config() + spectrum = BetaSpectrum.from_config(config) + W, _ = spectrum.get_energy_grid(config) + components = spectrum.calculate_components(W) + uni = universal_factor(components, W) + + branch_sum = np.zeros_like(W) + for i, branch in enumerate(branches_of(config)): + raw = components[f"branch_{i}.PhaseSpace"] + rad_key = f"branch_{i}.Radiative" + raw = raw * components[rad_key] if rad_key in components else raw + W0 = branch.endpoint_MeV / ME_MEV + 1.0 + raw = np.where(W <= W0, raw, 0.0) + branch_sum += branch.intensity * raw + + assert np.allclose(spectrum(W), uni * branch_sum, rtol=1e-10) + + +class TestAnalyzerPhysicalBranchSpectrum: + """analyzer.physical_branch_spectrum(i) = core × corrections, clamped.""" + + def test_multi_physical_includes_fermi(self) -> None: + config = make_multi_config() + spectrum = BetaSpectrum.from_config(config) + analyzer = BetaSpectrumAnalyzer(spectrum, config) + components = analyzer.components + W = analyzer.W + + uni = universal_factor(components, W) + for i, branch in enumerate(branches_of(config)): + raw = components[f"branch_{i}.PhaseSpace"].copy() + rad_key = f"branch_{i}.Radiative" + raw = raw * components[rad_key] if rad_key in components else raw + W0 = branch.endpoint_MeV / ME_MEV + 1.0 + raw = np.where(W <= W0, raw, 0.0) + physical = analyzer.physical_branch_spectrum(i) + expected = uni * raw + assert np.allclose(physical, expected, rtol=1e-10) + + def test_physical_clamped_beyond_endpoint(self) -> None: + config = make_multi_config() + spectrum = BetaSpectrum.from_config(config) + analyzer = BetaSpectrumAnalyzer(spectrum, config) + # Branch 0 (0.294 MeV) must be zero beyond its endpoint, while the + # grid extends to the maximum endpoint. + W = analyzer.W + W0_branch0 = branches_of(config)[0].endpoint_MeV / ME_MEV + 1.0 + W0_max = config.endpoint_MeV / ME_MEV + 1.0 + assert W0_branch0 <= W0_max + phys = analyzer.physical_branch_spectrum(0) + beyond = W > W0_branch0 + assert np.all(phys[beyond] == 0.0) + assert np.all(phys[~beyond] > 0.0) + + def test_sum_intensity_times_physical_equals_spectrum(self) -> None: + config = make_multi_config() + spectrum = BetaSpectrum.from_config(config) + analyzer = BetaSpectrumAnalyzer(spectrum, config) + total_raw = spectrum(analyzer.W) + weighted_sum = np.zeros_like(analyzer.W) + for i, branch in enumerate(branches_of(config)): + weighted_sum += branch.intensity * analyzer.physical_branch_spectrum(i) + assert np.allclose(weighted_sum, total_raw, rtol=1e-10) + + def test_single_branch_physical_equals_full_spectrum(self) -> None: + config = make_single_config() + spectrum = BetaSpectrum.from_config(config) + analyzer = BetaSpectrumAnalyzer(spectrum, config) + physical = analyzer.physical_branch_spectrum(0) + assert np.allclose(physical, spectrum(analyzer.W), rtol=1e-10) + + +class TestCSVExportCoreSemantics: + """CSV branch columns are physical spectra (core x corrections).""" + + def test_branch_columns_sum_to_total(self, tmp_path: Any) -> None: + config = make_multi_config() + spectrum = BetaSpectrum.from_config(config) + analyzer = BetaSpectrumAnalyzer(spectrum, config) + _, df = export_and_read(analyzer, tmp_path) + + total = df["total_spectrum"].to_numpy() + branch_sum = np.zeros_like(total) + for i in range(3): + branch_sum += df[f"branch_{i + 1}_spectrum"].to_numpy() + assert np.allclose(branch_sum, total, rtol=5e-4) + + def test_branch_columns_equal_physical_times_intensity(self, tmp_path: Any) -> None: + config = make_multi_config() + spectrum = BetaSpectrum.from_config(config) + analyzer = BetaSpectrumAnalyzer(spectrum, config) + _, df = export_and_read(analyzer, tmp_path) + + integral = float(np.trapezoid(spectrum(analyzer.W), analyzer.energies_MeV)) + for i, branch in enumerate(branches_of(config)): + physical = analyzer.physical_branch_spectrum(i) + expected = physical * branch.intensity / integral + got = df[f"branch_{i + 1}_spectrum"].to_numpy() + assert np.allclose(got, expected, rtol=5e-4) + + def test_total_column_integrates_to_one(self, tmp_path: Any) -> None: + config = make_multi_config() + spectrum = BetaSpectrum.from_config(config) + analyzer = BetaSpectrumAnalyzer(spectrum, config) + _, df = export_and_read(analyzer, tmp_path) + + area = np.trapezoid(df["total_spectrum"], df["energy_MeV"]) + assert np.isclose(area, 1.0, rtol=1e-3) + + def test_header_contains_core_line(self, tmp_path: Any) -> None: + config = make_multi_config() + spectrum = BetaSpectrum.from_config(config) + analyzer = BetaSpectrumAnalyzer(spectrum, config) + path, _ = export_and_read(analyzer, tmp_path) + + header = header_lines_of(path) + core_line = [h for h in header if h.startswith("# core:")] + assert core_line + assert "phase_space" in core_line[0] + assert "fermi" in core_line[0] + + def test_normalize_to_one_per_branch_unit_area(self, tmp_path: Any) -> None: + config = make_multi_config(normalize_to_one=True) + spectrum = BetaSpectrum.from_config(config) + analyzer = BetaSpectrumAnalyzer(spectrum, config) + _, df = export_and_read(analyzer, tmp_path) + + assert "total_spectrum" not in df.columns + for i in range(3): + col = df[f"branch_{i + 1}_spectrum"] + area = np.trapezoid(col, df["energy_MeV"]) + assert np.isclose(area, 1.0, rtol=1e-3) + + def test_degenerate_export_raises(self, tmp_path: Any) -> None: + config = make_degenerate_config() + spectrum = BetaSpectrum.from_config(config) + analyzer = BetaSpectrumAnalyzer(spectrum, config) + path = str(tmp_path / "degenerate.csv") + with pytest.raises(ValueError): + analyzer.export_to_csv(path) + + def test_concise_values_match_full(self, tmp_path: Any) -> None: + config = make_multi_config() + spectrum = BetaSpectrum.from_config(config) + analyzer = BetaSpectrumAnalyzer(spectrum, config) + _, df_full = export_and_read(analyzer, tmp_path, name="full.csv") + _, df_concise = export_and_read( + analyzer, tmp_path, name="concise.csv", csv_verbosity="concise" + ) + + assert set(df_concise.columns) == { + "energy_MeV", + "total_spectrum", + "branch_1", + "branch_2", + "branch_3", + } + assert np.allclose( + df_concise["branch_1"], df_full["branch_1_spectrum"], rtol=0, atol=0 + ) + assert np.allclose( + df_concise["total_spectrum"], df_full["total_spectrum"], rtol=0, atol=0 + ) + + def test_single_branch_branch_column_equals_total(self, tmp_path: Any) -> None: + config = make_single_config() + spectrum = BetaSpectrum.from_config(config) + analyzer = BetaSpectrumAnalyzer(spectrum, config) + _, df = export_and_read(analyzer, tmp_path) + + assert np.allclose( + df["branch_1_spectrum"], df["total_spectrum"], rtol=0, atol=0 + ) + + +class TestVisualizerBuiltFromCSV: + """Visualizer must plot from the exported CSV file contents only.""" + + def test_parse_spectrum_csv_metadata(self, tmp_path: Any) -> None: + config = make_multi_config() + spectrum = BetaSpectrum.from_config(config) + analyzer = BetaSpectrumAnalyzer(spectrum, config) + path, df = export_and_read(analyzer, tmp_path) + + parsed_df, meta = parse_spectrum_csv(path) + assert list(parsed_df.columns) == list(df.columns) + assert meta["version"] is not None + assert meta["timestamp"] is not None + assert meta["nuclide"] == "Tc99 -> Ru99" + assert meta["endpoint_keV"] == pytest.approx(294.0) + assert meta["e_step_MeV"] == pytest.approx(0.001) + assert meta["normalization"] == "intensity" + assert "screening" in meta["corrections"] + assert len(meta["branches"]) == 3 + assert meta["branches"][0]["intensity"] == pytest.approx(0.98) + assert meta["branches"][0]["excluded"] is None + assert meta["n_points"] == len(df) + + def test_plot_from_csv_curves_match_csv(self, tmp_path: Any) -> None: + config = make_multi_config() + spectrum = BetaSpectrum.from_config(config) + analyzer = BetaSpectrumAnalyzer(spectrum, config) + path, df = export_and_read(analyzer, tmp_path) + save = str(tmp_path / "plot.png") + + fig = Visualizer.plot_from_csv(path, save_path=save, show_components=True) + assert (tmp_path / "plot.png").exists() + + # debug view: 4 panels; panel 0 holds total + active branch curves + assert len(fig.axes) == 4 + panel = fig.axes[0] + curves = [np.asarray(ln.get_xydata()) for ln in panel.lines] + assert len(curves) == 1 + 3 + energies = df["energy_MeV"].to_numpy() + assert np.allclose(curves[0][:, 0], energies, rtol=0, atol=0) + assert np.allclose( + curves[0][:, 1], df["total_spectrum"].to_numpy(), rtol=0, atol=0 + ) + for i in range(3): + col = df[f"branch_{i + 1}_spectrum"].to_numpy() + assert np.allclose(curves[i + 1][:, 1], col, rtol=0, atol=0) + + def test_plot_from_csv_concise_falls_back_to_spectra( + self, tmp_path: Any, capsys: Any + ) -> None: + config = make_multi_config() + spectrum = BetaSpectrum.from_config(config) + analyzer = BetaSpectrumAnalyzer(spectrum, config) + path, _ = export_and_read( + analyzer, tmp_path, name="concise.csv", csv_verbosity="concise" + ) + save = str(tmp_path / "plot.png") + + fig = Visualizer.plot_from_csv(path, save_path=save, show_components=True) + assert len(fig.axes) == 1 + assert "Warning" in capsys.readouterr().out + + def test_plot_analysis_saves_figure(self, tmp_path: Any) -> None: + config = make_single_config() + spectrum = BetaSpectrum.from_config(config) + analyzer = BetaSpectrumAnalyzer(spectrum, config) + save = str(tmp_path / "analysis.png") + + Visualizer.plot_analysis(analyzer, save_path=save) + assert (tmp_path / "analysis.png").exists() + + +class TestAnalyzerUniversalFactor: + """analyzer.universal_factor() = product of components shared by all.""" + + def test_universal_factor_matches_component_product(self) -> None: + config = make_multi_config() + spectrum = BetaSpectrum.from_config(config) + analyzer = BetaSpectrumAnalyzer(spectrum, config) + uni = analyzer.universal_factor() + expected = universal_factor(analyzer.components, analyzer.W) + assert np.allclose(uni, expected, rtol=1e-12) + + def test_universal_factor_drops_disabled_corrections(self) -> None: + config = make_multi_config( + use_screening=False, + use_finite_size=False, + use_charge_dist=False, + use_exchange=False, + use_radiative=False, + ) + spectrum = BetaSpectrum.from_config(config) + analyzer = BetaSpectrumAnalyzer(spectrum, config) + uni = analyzer.universal_factor() + # Only the Fermi function remains universal: it equals the Fermi + # component alone. Fermi enhancement is > 1 for beta-minus Z=44. + assert np.all(uni > 1.0) + fermi_only = analyzer.components["Fermi"] + assert np.allclose(uni, fermi_only, rtol=1e-12) diff --git a/tests/physics/test_detector_response.py b/tests/physics/test_detector_response.py deleted file mode 100644 index c8a5a7a..0000000 --- a/tests/physics/test_detector_response.py +++ /dev/null @@ -1,334 +0,0 @@ -""" -Physics tests for DetectorResponse — unique physical behavior. - -Tests detector response shapes, resolution, and convolution physics. -""" - -import numpy as np -import pytest - -try: - from beta_spectrum.components.detector_response import DetectorResponse -except ImportError: - DetectorResponse = None # type: ignore[misc, assignment] - - -@pytest.mark.skipif( - DetectorResponse is None, reason="DetectorResponse not yet implemented" -) -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" diff --git a/tests/quality/test_api_detector_response.py b/tests/quality/test_api_detector_response.py deleted file mode 100644 index cfe4406..0000000 --- a/tests/quality/test_api_detector_response.py +++ /dev/null @@ -1,136 +0,0 @@ -""" -API / quality tests for DetectorResponse. - -Tests factory methods, repr, error handling. -""" - -import numpy as np -import pytest - -try: - from beta_spectrum.components.detector_response import DetectorResponse -except ImportError: - DetectorResponse = None # type: ignore[misc, assignment] - - -@pytest.mark.skipif( - DetectorResponse is None, reason="DetectorResponse not yet implemented" -) -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 - - -@pytest.mark.skipif( - DetectorResponse is None, reason="DetectorResponse not yet implemented" -) -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