test: pin core-spectrum invariant (phase_space x fermi always included)

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)

tests/physics/test_detector_response.py

@@ -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"