From ee7aedc80abdea2534609e3bc368aab513e41fe9 Mon Sep 17 00:00:00 2001 From: Evgeniy Unzhakov Date: Thu, 1 Oct 2026 22:04:17 +0300 Subject: [PATCH] fix(spectrum): lone branch defines endpoint in single-branch mode When config.branches carries exactly one branch (e.g. after --intensity-cutoff filtering), the branch endpoint/transition define the spectrum, energy grid, analyzer unified branch view and CSV header. The stale top-level endpoint inherited from the parent decay no longer leaks into physics (proven wrong on Cs137: spectrum extended to 1176 keV while the lone survivor ends at 514 keV). --- beta_spectrum/output/__init__.py | 3 -- beta_spectrum/output/core.py | 14 +++-- beta_spectrum/output/csv_export.py | 14 ++++- beta_spectrum/spectrum.py | 25 ++++++++- tests/physics/test_core_spectrum.py | 84 +++++++++++++++++++++++++++++ tests/quality/test_nuclear_data.py | 49 ++++++++++++++++- 6 files changed, 179 insertions(+), 10 deletions(-) diff --git a/beta_spectrum/output/__init__.py b/beta_spectrum/output/__init__.py index 0cc56d8..58e5b35 100644 --- a/beta_spectrum/output/__init__.py +++ b/beta_spectrum/output/__init__.py @@ -12,9 +12,6 @@ from .core import BetaSpectrumAnalyzer, _UnifiedBranch from .csv_export import CSVEmitter from .visualization import Visualizer -# Re-export the main class for backward compatibility -BetaSpectrumAnalyzer = BetaSpectrumAnalyzer - __all__ = [ "BetaSpectrumAnalyzer", "CSVEmitter", diff --git a/beta_spectrum/output/core.py b/beta_spectrum/output/core.py index fe88b0d..ded7f87 100644 --- a/beta_spectrum/output/core.py +++ b/beta_spectrum/output/core.py @@ -55,12 +55,20 @@ class BetaSpectrumAnalyzer: self._is_multi = spectrum._is_multi self._timestamp: Optional[str] = None - # Unified branch view: single-branch treated as one-branch multi-branch + # Unified branch view: single-branch treated as one-branch multi-branch. + # When config.branches carries exactly one branch, that branch (not the + # top-level decay fields) defines the endpoint and transition type. if not self._is_multi and self.spectrum.branch_spectra == []: + endpoint_MeV = config.endpoint_MeV + transition_type = config.transition_type + if config.branches is not None and len(config.branches) == 1: + lone = config.branches[0] + endpoint_MeV = float(lone.endpoint_MeV) + transition_type = str(lone.transition_type) self._unified_branches: List[_UnifiedBranch] = [ _UnifiedBranch( - endpoint_MeV=config.endpoint_MeV, - transition_type=config.transition_type, + endpoint_MeV=endpoint_MeV, + transition_type=transition_type, intensity=1.0, ) ] diff --git a/beta_spectrum/output/csv_export.py b/beta_spectrum/output/csv_export.py index 33bbd31..ec92c6b 100644 --- a/beta_spectrum/output/csv_export.py +++ b/beta_spectrum/output/csv_export.py @@ -81,6 +81,18 @@ class CSVEmitter: if config.use_exchange: enabled.append("exchange") + # The header endpoint must describe the spectrum that is actually + # computed. Single mode with exactly one configured branch follows + # the branch (the top-level endpoint may be stale from the parent + # decay, e.g. after intensity-cutoff filtering). + basic_endpoint_MeV = config.endpoint_MeV + if ( + not analyzer._is_multi + and config.branches is not None + and len(config.branches) == 1 + ): + basic_endpoint_MeV = analyzer._unified_branches[0].endpoint_MeV + parent_symbol = analyzer._element_symbol(config.Z_parent) daughter_symbol = analyzer._element_symbol(config.Z_daughter) @@ -99,7 +111,7 @@ class CSVEmitter: f"# nuclide: {parent_symbol}{config.A_number} " f"-> {daughter_symbol}{config.A_number} " f"(Z={config.Z_parent}->{config.Z_daughter})", - f"# basic_endpoint: {config.endpoint_MeV * 1000:.1f} keV", + f"# basic_endpoint: {basic_endpoint_MeV * 1000:.1f} keV", ] if config.neutrino_mass_MeV > 0: diff --git a/beta_spectrum/spectrum.py b/beta_spectrum/spectrum.py index f966b06..00286e9 100644 --- a/beta_spectrum/spectrum.py +++ b/beta_spectrum/spectrum.py @@ -95,6 +95,21 @@ class SpectrumConfig: ) +def _effective_single_fields(config: SpectrumConfig) -> Tuple[float, str]: + """Resolve (endpoint_MeV, transition_type) for single-branch mode. + + When config.branches carries exactly one branch (e.g. after intensity + cutoff filtering removed all but one), that branch defines the spectrum: + its endpoint and transition type win over the (possibly stale) top-level + fields inherited from the parent decay. With branches=None the top-level + fields are used as before. + """ + if config.branches is not None and len(config.branches) == 1: + lone = config.branches[0] + return float(lone.endpoint_MeV), str(lone.transition_type) + return float(config.endpoint_MeV), str(config.transition_type) + + class BetaSpectrum: """ Precise beta spectrum calculator. @@ -259,7 +274,8 @@ class BetaSpectrum: conditional. """ components: List[SpectrumComponent] = [] - W0 = float(T_to_W(config.endpoint_MeV)) + eff_endpoint, eff_transition = _effective_single_fields(config) + W0 = float(T_to_W(eff_endpoint)) # Core spectrum components (0-order): always included. # Convert neutrino mass from MeV to natural units @@ -267,7 +283,7 @@ class BetaSpectrum: components.append( PhaseSpace( W0=W0, - transition_type=config.transition_type, + transition_type=eff_transition, m_nu=m_nu_natural, logger=logger, ) @@ -419,6 +435,11 @@ class BetaSpectrum: """ if self._is_multi and self.branches: max_endpoint = max(b.endpoint_MeV for b in self.branches) + elif config.branches is not None and len(config.branches) == 1: + # Single-mode with exactly one configured branch: the branch + # endpoint defines the grid, not the (possibly stale) top-level + # endpoint from the parent decay. + max_endpoint = float(config.branches[0].endpoint_MeV) else: max_endpoint = config.endpoint_MeV diff --git a/tests/physics/test_core_spectrum.py b/tests/physics/test_core_spectrum.py index 652cdf7..f9f8ed9 100644 --- a/tests/physics/test_core_spectrum.py +++ b/tests/physics/test_core_spectrum.py @@ -438,6 +438,90 @@ class TestVisualizerBuiltFromCSV: assert (tmp_path / "analysis.png").exists() +def make_single_survivor_config(**over: Any) -> SpectrumConfig: + """A config whose top-level endpoint exceeds a lone branch's endpoint. + + This mirrors the intensity_cutoff / decay_info_to_config path where the + parent endpoint_MeV keeps the strongest-decay value while the only + surviving branch sits far below it. + """ + over.setdefault("endpoint_MeV", 1.176) + over.setdefault("transition_type", "F2") + over.setdefault( + "branches", + [BranchConfig(endpoint_MeV=0.514341, transition_type="A", intensity=1.0)], + ) + return SpectrumConfig( + Z_parent=55, + Z_daughter=56, + A_number=137, + e_step_MeV=0.001, + **over, + ) + + +class TestSurvivingBranchEndpointFollows: + """When config.branches has exactly 1 entry, single-mode spectrum, + energy grid and analyzer must follow the branch endpoint/transition, + not the top-level endpoint fields.""" + + def test_single_mode_builds_from_branch_endpoint(self) -> None: + config = make_single_survivor_config() + spectrum = BetaSpectrum.from_config(config) + assert not spectrum._is_multi + W, E = spectrum.get_energy_grid(config) + # The energy grid must be capped at the lone branch endpoint + # (~514 keV), not the stale top-level endpoint (1176 keV), and + # stay non-negative everywhere it covers. + assert float(E.max()) < 0.514341 + 1e-12 + values = spectrum(W) + assert np.all(np.isfinite(values)) + assert np.all(values >= 0.0) + + def test_single_mode_matches_reference_bare_config(self) -> None: + """Branch-driven single mode == a bare config set to the branch endpoint.""" + branch_ep = 0.514341 + cfg_branch = make_single_survivor_config() + cfg_ref = SpectrumConfig( + Z_parent=55, + Z_daughter=56, + A_number=137, + endpoint_MeV=branch_ep, + transition_type="A", + e_step_MeV=0.001, + ) + spec_b = BetaSpectrum.from_config(cfg_branch) + spec_ref = BetaSpectrum.from_config(cfg_ref) + Wb, Eb = spec_b.get_energy_grid(cfg_branch) + Wr, Er = spec_ref.get_energy_grid(cfg_ref) + vb = spec_b(Wb) + vr = spec_ref(Wr) + # The branch-driven grid and values must be identical to a bare + # config built directly at the branch endpoint/transition. + assert len(Eb) == len(Er) + assert np.allclose(Eb, Er, rtol=0, atol=1e-12) + assert np.allclose(vb, vr, rtol=1e-10) + + def test_analyzer_unified_view_follows_branch(self) -> None: + config = make_single_survivor_config() + spectrum = BetaSpectrum.from_config(config) + analyzer = BetaSpectrumAnalyzer(spectrum, config) + ub = analyzer._unified_branches[0] + assert ub.endpoint_MeV == pytest.approx(0.514341) + assert ub.transition_type == "A" + + def test_csv_header_reports_branch_endpoint(self, tmp_path: Any) -> None: + config = make_single_survivor_config() + spectrum = BetaSpectrum.from_config(config) + analyzer = BetaSpectrumAnalyzer(spectrum, config) + path, _ = export_and_read(analyzer, tmp_path) + header = header_lines_of(path) + # The lone surviving branch must be reported at its own endpoint, + # not the 1176 keV top-level endpoint. + assert any("514.3 keV" in h for h in header) + assert not any("1176.0 keV" in h for h in header) + + class TestAnalyzerUniversalFactor: """analyzer.universal_factor() = product of components shared by all.""" diff --git a/tests/quality/test_nuclear_data.py b/tests/quality/test_nuclear_data.py index 0c95aee..2e54089 100644 --- a/tests/quality/test_nuclear_data.py +++ b/tests/quality/test_nuclear_data.py @@ -5,6 +5,7 @@ import tempfile from pathlib import Path from typing import Any, Dict +import numpy as np import pytest # ruff: noqa: E402, ANN201, ANN001 @@ -19,7 +20,7 @@ from beta_spectrum.nuclear_data import ( # noqa: E402 validate_json_input, _parse_nuclide_symbol, ) -from beta_spectrum.spectrum import SpectrumConfig +from beta_spectrum.spectrum import BetaSpectrum, SpectrumConfig def _has_paceENSDF() -> bool: # noqa: ANN201 @@ -241,6 +242,52 @@ class TestDecayInfoToConfig: assert len(config.branches) == 1 assert config.branches[0].intensity == 0.5 + def test_cutoff_single_survivor_endpoint_follows_branch(self): + """A lone surviving branch must define the spectrum, not the stale + top-level endpoint inherited from the parent decay (D1).""" + info = DecayInfo( + parent_symbol="Cs137", + Z_parent=55, + Z_daughter=56, + A_number=137, + endpoint_keV=1176.0, + endpoint_MeV=1.176, + transition_type="F2", + forbiddenness_code="2F", + parent_J=7.0, + parent_pi=1, + half_life="y", + branches=[ + BranchInfo( + level_index=0, + level_energy_keV=0.0, + intensity=0.053, + log_ft=5.0, + transition_type="F2", + ), + BranchInfo( + level_index=1, + level_energy_keV=661.66, + intensity=0.947, + log_ft=5.0, + transition_type="A", + ), + ], + ) + cfg = decay_info_to_config(info, intensity_cutoff=0.9) + assert cfg.branches is not None + assert len(cfg.branches) == 1 + lone_ep = cfg.branches[0].endpoint_MeV + assert lone_ep == pytest.approx(1.176 - 0.66166) + spectrum = BetaSpectrum.from_config(cfg) + assert not spectrum._is_multi + W, E = spectrum.get_energy_grid(cfg) + # Grid capped at the surviving branch endpoint, not 1.176 MeV. + assert float(E.max()) < lone_ep + 1e-12 + values = spectrum(W) + assert np.all(np.isfinite(values)) + assert np.all(values >= 0.0) + # --------------------------------------------------------------------------- # get_decay_info_from_paceENSDF