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).
This commit is contained in:
2026-10-01 22:04:17 +03:00
parent 6bb7b950c7
commit ee7aedc80a
6 changed files with 179 additions and 10 deletions
-3
View File
@@ -12,9 +12,6 @@ from .core import BetaSpectrumAnalyzer, _UnifiedBranch
from .csv_export import CSVEmitter from .csv_export import CSVEmitter
from .visualization import Visualizer from .visualization import Visualizer
# Re-export the main class for backward compatibility
BetaSpectrumAnalyzer = BetaSpectrumAnalyzer
__all__ = [ __all__ = [
"BetaSpectrumAnalyzer", "BetaSpectrumAnalyzer",
"CSVEmitter", "CSVEmitter",
+11 -3
View File
@@ -55,12 +55,20 @@ class BetaSpectrumAnalyzer:
self._is_multi = spectrum._is_multi self._is_multi = spectrum._is_multi
self._timestamp: Optional[str] = None 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 == []: 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] = [ self._unified_branches: List[_UnifiedBranch] = [
_UnifiedBranch( _UnifiedBranch(
endpoint_MeV=config.endpoint_MeV, endpoint_MeV=endpoint_MeV,
transition_type=config.transition_type, transition_type=transition_type,
intensity=1.0, intensity=1.0,
) )
] ]
+13 -1
View File
@@ -81,6 +81,18 @@ class CSVEmitter:
if config.use_exchange: if config.use_exchange:
enabled.append("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) parent_symbol = analyzer._element_symbol(config.Z_parent)
daughter_symbol = analyzer._element_symbol(config.Z_daughter) daughter_symbol = analyzer._element_symbol(config.Z_daughter)
@@ -99,7 +111,7 @@ class CSVEmitter:
f"# nuclide: {parent_symbol}{config.A_number} " f"# nuclide: {parent_symbol}{config.A_number} "
f"-> {daughter_symbol}{config.A_number} " f"-> {daughter_symbol}{config.A_number} "
f"(Z={config.Z_parent}->{config.Z_daughter})", 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: if config.neutrino_mass_MeV > 0:
+23 -2
View File
@@ -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: class BetaSpectrum:
""" """
Precise beta spectrum calculator. Precise beta spectrum calculator.
@@ -259,7 +274,8 @@ class BetaSpectrum:
conditional. conditional.
""" """
components: List[SpectrumComponent] = [] 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. # Core spectrum components (0-order): always included.
# Convert neutrino mass from MeV to natural units # Convert neutrino mass from MeV to natural units
@@ -267,7 +283,7 @@ class BetaSpectrum:
components.append( components.append(
PhaseSpace( PhaseSpace(
W0=W0, W0=W0,
transition_type=config.transition_type, transition_type=eff_transition,
m_nu=m_nu_natural, m_nu=m_nu_natural,
logger=logger, logger=logger,
) )
@@ -419,6 +435,11 @@ class BetaSpectrum:
""" """
if self._is_multi and self.branches: if self._is_multi and self.branches:
max_endpoint = max(b.endpoint_MeV for b in 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: else:
max_endpoint = config.endpoint_MeV max_endpoint = config.endpoint_MeV
+84
View File
@@ -438,6 +438,90 @@ class TestVisualizerBuiltFromCSV:
assert (tmp_path / "analysis.png").exists() 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: class TestAnalyzerUniversalFactor:
"""analyzer.universal_factor() = product of components shared by all.""" """analyzer.universal_factor() = product of components shared by all."""
+48 -1
View File
@@ -5,6 +5,7 @@ import tempfile
from pathlib import Path from pathlib import Path
from typing import Any, Dict from typing import Any, Dict
import numpy as np
import pytest import pytest
# ruff: noqa: E402, ANN201, ANN001 # ruff: noqa: E402, ANN201, ANN001
@@ -19,7 +20,7 @@ from beta_spectrum.nuclear_data import ( # noqa: E402
validate_json_input, validate_json_input,
_parse_nuclide_symbol, _parse_nuclide_symbol,
) )
from beta_spectrum.spectrum import SpectrumConfig from beta_spectrum.spectrum import BetaSpectrum, SpectrumConfig
def _has_paceENSDF() -> bool: # noqa: ANN201 def _has_paceENSDF() -> bool: # noqa: ANN201
@@ -241,6 +242,52 @@ class TestDecayInfoToConfig:
assert len(config.branches) == 1 assert len(config.branches) == 1
assert config.branches[0].intensity == 0.5 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 # get_decay_info_from_paceENSDF