style: canonicalize black formatting on test files

- --normalize-to-one renamed to --separate-branches (no alias; the
  name now describes the output choice: total spectrum with
  intensity-normalized branches vs separate branches, each
  normalized to unit area with intensities left to the data user).
  SpectrumConfig.normalize_to_one -> separate_branches, CSV header
  word to_one -> separate.
- every long option gains a single-letter alias: -N -i -s -m -o -p
  -L -d -c -V -S -b (only -v/-q existed before)
- dry-run drops the top-level Transition: line — the per-branch
  list below already reports each branch's transition type
- version bumped to 0.4.1 (pyproject, __version__, CSV semantics
  docstring, TODO roadmap entry)

.gitignore

@@ -180,3 +180,6 @@ htmlcov/
 .rpiv/
 .pi-lens/
 data/raw/
+
+# ROOT data files produced by notebooks
+*.root

README.md

@@ -60,8 +60,7 @@ beta_shape_pnpi/
 │       │   ├── finite_size.py           ✓ L0 expansion + charge distribution U term
 │       │   ├── screening.py             ✓ Atomic electron screening (ratio method)
 │       │   ├── exchange.py              ✓ Hayen-2018 empirical fit coefficients
-│       │   ├── radiative.py             ✓ Outer radiative corrections, soft-photon resummation
-│       │   └── detector_response.py     ✓ Analytical detector response models
+│       │   └── radiative.py             ✓ Outer radiative corrections, soft-photon resummation
 ├── beta_spectrum/output/    ← output, plotting, visualization, CSV export
 │   ├── __init__.py          # BetaSpectrumAnalyzer, CSVEmitter, Visualizer
 │   ├── core.py              # BetaSpectrumAnalyzer base class
@@ -122,12 +121,19 @@ Defined in `beta_spectrum/spectrum.py`. Orchestrates the calculation:
 
 In `beta_spectrum.output`. Provides debugging and visualization tools:
 
-- `total_spectrum(normalize=True/False)` — computed spectrum array
+- `total_spectrum(normalize=True/False)` — computed spectrum array; raises `ValueError` when every branch is degenerate (non-positive integral)
+- `physical_branch_spectrum(i)` — physical spectrum of branch *i*: the core (phase space × Fermi function) times every enabled higher-order correction, clamped beyond the branch endpoint (unnormalized)
+- `universal_factor()` — product of the components shared across all branches (Fermi function and the universal higher-order corrections)
 - `plot_analysis(save_path=None, show_components=True)` — plot with two modes:
   - `show_components=True` (debug mode): 4-panel figure (spectrum, components, cumulative effect, deviations)
   - `show_components=False` (normal mode): single spectrum plot with nuclear data header
   - All plots include commit hash and UTC timestamp for traceability
-- `export_to_csv(filename)` — export via pandas DataFrame with YAML-style metadata header
+  - Plots are built *from the exported CSV file* (temporary full-verbosity export), so a plot always shows exactly the bytes the CSV export produces
+- `export_to_csv(filename)` — export via pandas DataFrame with YAML-style metadata header:
+  - the basic spectrum is always the core `phase_space × fermi` (the header carries a `# core:` line); finite size, charge distribution, screening, exchange and radiative are higher-order corrections
+  - `branch_{i}_spectrum` columns contain the *physical* branch spectra weighted by intensity and normalized so that the sum of branch columns equals `total_spectrum`
+  - with `separate_branches=True` (CLI `--separate-branches`) each branch column integrates to 1 (intensity ignored, no total_spectrum column)
+- `parse_spectrum_csv(path)` — read an exported CSV and its metadata header (module function in `beta_spectrum.output.visualization`)
 - `get_data()` — return all numerical data for custom analysis
 
 ______________________________________________________________________
@@ -184,9 +190,9 @@ ______________________________________________________________________
 For a complete list of implemented and planned features, see [TODO.md](TODO.md). Key highlights:
 
 - **Physics corrections**: phase space, Fermi function, finite size, screening, exchange, radiative (with delta_cut resummation)
-- **Detector response**: Gaussian, Gaussian+tail, Tikhonov models with tabulated support
+- **Detector response**: not implemented — `SpectrumModel` keeps a duck-typed optional convolution hook; the analytical models were removed (restore-or-delete decision tracked in TODO.md)
 - **Multi-branch decay**: intensity-weighted sum of branch spectra, per-branch calculators with branch-specific endpoint, transition type (from ENSDF), and `--intensity-cutoff` CLI filter
 - **Analysis tools**: χ² fitter (SpectrumFitter + CurveFitter), C(W) extraction pipeline, Kurie plot analysis, g_V/g_A extraction
 - **Experimental data**: ExpSpectrum container, energy calibration, peak fitting, dead-time/pile-up correction
-- **CLI**: `bs_pnpi --nuclide Tc99` — full paceENSDF integration (auto-lookup of Z, A, Q-value, transition type, half-life, branches, forbiddenness from ENSDF)
-- **paceENSDF integration**: `get_decay_info_from_paceENSDF()`, `create_config_from_source("paceENSDF", nuclide="Tc99")`, `_parse_nuclide_symbol()` for Z=1..98, FORBIDDENNESS_MAP, `--nuclide` CLI flag
+- **CLI**: `bs_pnpi --nuclide Tc99` — full paceENSDF integration (auto-lookup of Z, A, Q-value, transition type, half-life, branches, forbiddenness from ENSDF); for JSON input the file's values are used and CLI flags (`--e-step`, `--neutrino-mass`, `--intensity-cutoff`, `--csv-separator`, `--separate-branches`) override them only when explicitly given; every long option has a single-letter alias (`-N`, `-i`, `-s`, `-m`, `-o`, `-p`, `-L`, `-d`, `-c`, `-V`, `-S`, `-b`)
+- **paceENSDF integration**: `get_decay_info_from_paceENSDF()`, `create_config_from_source("paceENSDF", nuclide="Tc99")`, `_parse_nuclide_symbol()` for Z=1..98, FORBIDDENNESS_MAP, `--nuclide` CLI flag; per-branch endpoints are the measured transition endpoint energies from ENSDF (max endpoint across branches defines the spectrum endpoint)

TODO.md

@@ -32,7 +32,7 @@ beta_spectrum/              ← theory-only core
 │   ├── screening.py
 │   ├── exchange.py
 │   ├── radiative.py
-│   └── detector_response.py   ← stays here (part of theory)
+│   └── (detector_response.py — removed; restore-or-delete pending below)
 └── nuclear_data.py         ← paceENSDF integration
 
 exp_data/                   ← experimental data handling
@@ -101,6 +101,11 @@ output/                     ← output, visualization, CSV export
 - [x] 322 tests pass
 - [x] README updated with new architecture diagram and Quick Start examples
 - [x] Version bumped to 0.4.0
+- [ ] Decide: restore or delete the `detector_response` module — it was removed from
+  `beta_spectrum/components/` together with its tests; `fitter/model.py` keeps a
+  duck-typed optional hook only (`DEFAULT_JSON_SCHEMA` detector keys and the CLI
+  dry-run detector line were removed with audit cleanup).
+  Restore from git history (commits up to `5f99566`) if convolution is needed again
 
 ______________________________________________________________________
 
@@ -317,7 +322,7 @@ ______________________________________________________________________
 
 ## Current Status
 
-**Version:** 0.3.1 → **In progress: 0.4.0 — Package restructuring**
+**Version:** 0.4.0 → **v0.4.1 — CLI options refresh & normalization fixes** (released on `dev/core-fermi-csv`)
 
 ### Completed (v0.3.x)
 
@@ -349,6 +354,12 @@ ______________________________________________________________________
   - README updated with new architecture diagram and Quick Start examples
   - Version bumped to 0.4.0
 
+### Completed (v0.4.1)
+
+- Core spectrum made unconditional (`phase_space × fermi`): core toggles removed, physical branch spectra in CSV export, plots built from the exported CSV
+- Normalization audit fixes: lone surviving branch defines the spectrum, CLI options propagate through the JSON flow, positive-integral guard for both export modes, paceENSDF per-branch measured endpoints (Co60 fixed), single delegated paceENSDF config builder
+- CLI surface refresh: `--normalize-to-one` renamed to `--separate-branches` (no alias), single-letter aliases for every option, dry-run drops the redundant top-level `Transition:` line, stale detector_* JSON schema keys removed
+
 ### Next Immediate Step
 
 - **Track A: ⁹⁹Tc spectrum-shape analysis** — start real data processing pipeline

beta_spectrum/__init__.py

@@ -52,4 +52,4 @@ __all__ = [
     "DEFAULT_JSON_SCHEMA",
 ]
 
-__version__ = "0.4.0"
+__version__ = "0.4.1"

beta_spectrum/cli.py

@@ -40,7 +40,6 @@ def _build_parser() -> argparse.ArgumentParser:
             "  bs_pnpi --nuclide Tc99 --output spectrum.csv\n"
             "  bs_pnpi --input custom_input.json --output spectrum.csv\n"
             "  bs_pnpi --nuclide Tc99 --plot spectrum.png --output spectrum.csv\n"
-            "  bs_pnpi --nuclide Tc99 --detector gaussian --output convolved.csv\n"
             "  bs_pnpi --nuclide Tc99 --output spectrum.csv -vv --log-file /tmp/calc.log\n"
             "  bs_pnpi --nuclide Tc99 --dry-run\n"
         ),
@@ -66,12 +65,14 @@ def _build_parser() -> argparse.ArgumentParser:
     # Source selection flags (backward compatible, mutually exclusive)
     source_group = parser.add_mutually_exclusive_group(required=False)
     source_group.add_argument(
+        "-N",
         "--nuclide",
         type=str,
         help="Nuclide symbol for paceENSDF lookup (e.g. 'Tc99', 'Co60'). "
         "Retrieves decay parameters from the ENSDF database via paceENSDF.",
     )
     source_group.add_argument(
+        "-i",
         "--input",
         type=str,
         metavar="FILE",
@@ -80,20 +81,27 @@ def _build_parser() -> argparse.ArgumentParser:
 
     # Common parameters
     parser.add_argument(
+        "-s",
         "--e-step",
         type=float,
-        default=0.001,
-        help="Energy step size in MeV (default: 0.001). All energies are in MeV throughout the toolkit.",
+        default=None,
+        help="Energy step size in MeV (default: 0.001 for nuclide input; "
+        "for JSON input the file's e_step_MeV is used when the flag is "
+        "omitted). All energies are in MeV throughout the toolkit.",
     )
     parser.add_argument(
+        "-m",
         "--neutrino-mass",
         type=float,
-        default=0.0,
-        help="Neutrino mass in MeV (default: 0.0). Set to non-zero value to include neutrino mass effects in phase space calculation.",
+        default=None,
+        help="Neutrino mass in MeV (default: 0.0, or the JSON file's "
+        "neutrino_mass_MeV when omitted). Set to non-zero value to include "
+        "neutrino mass effects in phase space calculation.",
     )
 
     # Output
     parser.add_argument(
+        "-o",
         "--output",
         type=str,
         metavar="FILE",
@@ -102,6 +110,7 @@ def _build_parser() -> argparse.ArgumentParser:
 
     # Plotting
     parser.add_argument(
+        "-p",
         "--plot",
         type=str,
         metavar="FILE",
@@ -125,6 +134,7 @@ def _build_parser() -> argparse.ArgumentParser:
         help="Suppress all terminal output (only errors shown).",
     )
     log_group.add_argument(
+        "-L",
         "--log-file",
         type=str,
         metavar="PATH",
@@ -132,20 +142,23 @@ def _build_parser() -> argparse.ArgumentParser:
         "Log file is always opened in append mode at DEBUG level.",
     )
     log_group.add_argument(
+        "-d",
         "--dry-run",
         action="store_true",
         help="Validate input and display resolved configuration without computing spectrum.",
     )
     log_group.add_argument(
+        "-c",
         "--intensity-cutoff",
         type=float,
-        default=0.0,
+        default=None,
         metavar="FRAC",
         help="Minimum branch intensity fraction to include in the spectrum "
         "(default: 0.0 = include all branches). "
         "Branches below this threshold are excluded from the calculation and plots.",
     )
     log_group.add_argument(
+        "-V",
         "--csv-verbosity",
         type=str,
         default="full",
@@ -153,15 +166,20 @@ def _build_parser() -> argparse.ArgumentParser:
         help="CSV output verbosity: 'full' (all components) or 'concise' (total spectrum and normalized branch components only) (default: full)",
     )
     log_group.add_argument(
+        "-S",
         "--csv-separator",
         type=str,
-        default=",",
+        default=None,
         help="CSV column separator (default: ','). Use 'tab' for tab separation.",
     )
     log_group.add_argument(
-        "--normalize-to-one",
+        "-b",
+        "--separate-branches",
         action="store_true",
-        help="Normalize each branch to 1, ignoring intensity. In this case, total spectrum will not be produced (only energy and branch-wise spectra).",
+        default=None,
+        help="Emit separate branch spectra, each normalized to integrate to "
+        "1 (intensities ignored; no total_spectrum column). Default: total "
+        "spectrum with intensity-normalized branch columns.",
     )
 
     return parser
@@ -174,15 +192,10 @@ def _run(args: argparse.Namespace) -> None:
     from beta_spectrum.logging_utils import get_logger
     from beta_spectrum.nuclear_data import (
         create_config_from_source,
-        get_decay_info_from_paceENSDF,
         json_to_config,
         load_json_input,
     )
-    from beta_spectrum.spectrum import (
-        BetaSpectrum,
-        BranchConfig,
-        SpectrumConfig,
-    )
+    from beta_spectrum.spectrum import BetaSpectrum
     from beta_spectrum.output import BetaSpectrumAnalyzer
 
     # Determine logging level
@@ -231,66 +244,45 @@ def _run(args: argparse.Namespace) -> None:
     else:
         source_type = "paceENSDF"
 
-    # Build config
+    # Resolve optional CLI flags: None means "not provided". For paceENSDF
+    # they fall back to literal defaults; for JSON they let the file's own
+    # values through unless explicitly overridden.
+    e_step_MeV = 0.001 if args.e_step is None else float(args.e_step)
+    neutrino_mass_MeV = 0.0 if args.neutrino_mass is None else float(args.neutrino_mass)
+    intensity_cutoff = (
+        0.0 if args.intensity_cutoff is None else float(args.intensity_cutoff)
+    )
+    csv_separator = "," if args.csv_separator is None else args.csv_separator
+    separate_branches = bool(args.separate_branches)
+
+    # Build config — a single delegated path for every data source.
     config: SpectrumConfig
     if source_type == "paceENSDF":
         logger.info("Source: paceENSDF nuclide=%s", source_value)
-        decay_info = get_decay_info_from_paceENSDF(source_value, "beta_minus")
-
-        # Normalize branch intensities to sum to 1.0
-        if decay_info.branches:
-            total_intensity = sum(b.intensity for b in decay_info.branches)
-            if total_intensity > 0:
-                branches = [
-                    BranchConfig(
-                        endpoint_MeV=(
-                            (decay_info.endpoint_MeV - b.level_energy_keV / 1000.0)
-                            if b.level_energy_keV > 0
-                            else decay_info.endpoint_MeV
-                        ),
-                        transition_type=b.transition_type,
-                        intensity=b.intensity / total_intensity,
-                    )
-                    for b in decay_info.branches
-                    if b.intensity / total_intensity >= args.intensity_cutoff
-                ]
-                config = SpectrumConfig(
-                    Z_parent=decay_info.Z_parent,
-                    Z_daughter=decay_info.Z_daughter,
-                    A_number=decay_info.A_number,
-                    endpoint_MeV=decay_info.endpoint_MeV,
-                    transition_type=decay_info.transition_type,
-                    e_step_MeV=args.e_step,
-                    neutrino_mass_MeV=args.neutrino_mass,
-                    branches=branches,
-                    intensity_cutoff=args.intensity_cutoff,
-                    csv_separator=args.csv_separator,
-                    normalize_to_one=args.normalize_to_one,
-                )
-            else:
         config = create_config_from_source(
             "paceENSDF",
             nuclide=source_value,
-                    e_step_MeV=args.e_step,
-                    intensity_cutoff=args.intensity_cutoff,
-                    neutrino_mass_MeV=args.neutrino_mass,
-                    csv_separator=args.csv_separator,
-                    normalize_to_one=args.normalize_to_one,
-                )
-        else:
-            config = create_config_from_source(
-                "paceENSDF",
-                nuclide=source_value,
-                e_step_MeV=args.e_step,
-                intensity_cutoff=args.intensity_cutoff,
-                neutrino_mass_MeV=args.neutrino_mass,
-                csv_separator=args.csv_separator,
-                normalize_to_one=args.normalize_to_one,
+            e_step_MeV=e_step_MeV,
+            intensity_cutoff=intensity_cutoff,
+            neutrino_mass_MeV=neutrino_mass_MeV,
+            csv_separator=csv_separator,
+            separate_branches=separate_branches,
         )
     else:
         logger.info("Source: JSON file=%s", source_value)
         data = load_json_input(source_value)
         config = json_to_config(data)
+        # CLI flags override the JSON file values only when explicitly
+        # provided; unset flags keep the file's own values.
+        if args.e_step is not None:
+            config.e_step_MeV = e_step_MeV
+        if args.neutrino_mass is not None:
+            config.neutrino_mass_MeV = neutrino_mass_MeV
+        if args.intensity_cutoff is not None:
+            config.intensity_cutoff = intensity_cutoff
+        if args.csv_separator is not None:
+            config.csv_separator = csv_separator
+        config.separate_branches = separate_branches
 
     # Dry run: display resolved config and exit
     if args.dry_run:
@@ -326,9 +318,20 @@ def _run(args: argparse.Namespace) -> None:
         )
         logger.info("Spectrum exported to %s", args.output)
 
-    # Generate plot
+    # Generate plot — always built from the exported CSV data (single
+    # source of truth): when a full-verbosity CSV was exported it is
+    # plotted directly from that file, otherwise via a temporary export.
     if args.plot:
         analyzer = BetaSpectrumAnalyzer(spectrum, config, logger=logger)
+        if args.output and args.csv_verbosity == "full":
+            from beta_spectrum.output.visualization import Visualizer
+
+            Visualizer.plot_from_csv(
+                args.output,
+                save_path=args.plot,
+                show_components=(args.verbose >= 2),
+            )
+        else:
             analyzer.plot_analysis(
                 save_path=args.plot,
                 show_components=(args.verbose >= 2),
@@ -342,11 +345,8 @@ def _print_dry_run_output(config, source_type: str) -> None:
     """Print resolved configuration for dry-run mode."""
     from beta_spectrum.logging_utils import get_git_short_hash
 
+    print("Core:       phase_space × fermi (basic spectrum, always included)")
     enabled = []
-    if config.use_phase_space:
-        enabled.append("phase_space")
-    if config.use_fermi:
-        enabled.append("fermi")
     if config.use_screening:
         enabled.append("screening")
     if config.use_finite_size:
@@ -366,7 +366,6 @@ def _print_dry_run_output(config, source_type: str) -> None:
         f"Nuclide:    {_decay_notation(config.Z_parent, config.Z_daughter, config.A_number)}"
     )
     print(f"Endpoint:   {config.endpoint_MeV * 1000:.1f} keV")
-    print(f"Transition: {config.transition_type}")
     print(f"e_step:     {config.e_step_MeV:.4f} MeV")
     print(f"neutrino_mass_MeV: {config.neutrino_mass_MeV:.6f}")
     print(f"Corrections: {', '.join(enabled) if enabled else 'none'}")
@@ -386,8 +385,6 @@ def _print_dry_run_output(config, source_type: str) -> None:
     else:
         print("Mode: single-branch")
 
-    # Detector response is a Phase 4 feature — not yet implemented
-    print("Detector:   not configured (Phase 4)")
     print("========================================")
 
 
@@ -509,10 +506,16 @@ def _decay_notation(Z_parent: int, Z_daughter: int, A: int) -> str:
     return f"{_nuclide_symbol(Z_parent, A)} -> {_nuclide_symbol(Z_daughter, A)}"
 
 
-def main() -> None:
-    """Main entry point for the CLI."""
+def main(argv: list[str] | None = None) -> None:
+    """Main entry point for the CLI.
+
+    Parameters
+    ----------
+    argv : list[str], optional
+        Explicit argument vector (for tests/embedding); None reads sys.argv.
+    """
     parser = _build_parser()
-    args = parser.parse_args()
+    args = parser.parse_args(argv)
 
     try:
         _run(args)

beta_spectrum/nuclear_data.py

@@ -96,6 +96,10 @@ class BranchInfo:
         log(ft) value for this branch.
     transition_type : str
         Forbiddenness classification for this branch (A, F1, F1U, F2, ...).
+    endpoint_MeV : float
+        Measured beta endpoint energy of this transition in MeV (ENSDF
+        column 7). 0.0 means "not provided": consumers fall back to
+        base-endpoint minus level-energy subtraction.
     """
 
     level_index: int
@@ -103,6 +107,7 @@ class BranchInfo:
     intensity: float
     log_ft: Optional[float]
     transition_type: str = "A"
+    endpoint_MeV: float = 0.0
 
 
 def _parse_nuclide_symbol(symbol: str) -> Tuple[str, int, int]:
@@ -384,24 +389,19 @@ def get_decay_info_from_paceENSDF(
 
     # Extract Q-value from ground state branch (state[7] = endpoint energy in keV)
     # and half-life from halflife_data
-    endpoint_keV = 0.0
     half_life_str = "unknown"
 
+    # The spectrum endpoint is the maximum measured transition endpoint
+    # across all branches (ENSDF column 7, keV). Selecting the ground
+    # state branch or the first branch is wrong whenever those branches
+    # are absent or non-maximal (e.g. Co60, whose first listed branch
+    # endpoint is below higher-level transition endpoints).
+    endpoint_keV = 0.0
     if beta_data:
-        for k, v in beta_data.items():
+        for _k, v in beta_data.items():
             for state in v:
-                if state[0] == 0:  # ground state branch
-                    endpoint_keV = float(state[7])
-                    break
-            if endpoint_keV > 0:
-                break
-
-        # No ground state branch — use first available branch
-        if endpoint_keV == 0.0:
-            for k, v in beta_data.items():
-                if v and len(v) > 0:
-                    endpoint_keV = float(v[0][7])
-                    break
+                if len(state) > 7 and state[7] is not None:
+                    endpoint_keV = max(endpoint_keV, float(state[7]))
 
     if halflife_data:
         for k, v in halflife_data.items():
@@ -447,6 +447,11 @@ def get_decay_info_from_paceENSDF(
                 intensity=float(state[9]) / 100.0,
                 log_ft=float(state[11]) if state[11] is not None else None,
                 transition_type=branch_transition_type,
+                endpoint_MeV=(
+                    float(state[7]) / 1000.0
+                    if len(state) > 7 and state[7] is not None
+                    else 0.0
+                ),
             )
             branches.append(branch)
 
@@ -472,7 +477,7 @@ def decay_info_to_config(
     intensity_cutoff: float = 0.0,
     neutrino_mass_MeV: float = 0.0,
     csv_separator: str = ",",
-    normalize_to_one: bool = False,
+    separate_branches: bool = False,
 ) -> SpectrumConfig:
     """
     Convert DecayInfo to a SpectrumConfig.
@@ -489,20 +494,28 @@ def decay_info_to_config(
         Neutrino mass in MeV.
     csv_separator : str
         CSV column separator.
-    normalize_to_one : bool
-        If True, normalize each branch to 1, ignoring intensity.
+    separate_branches : bool
+        If True, emit separate branch spectra, each normalized to
+        integrate to 1 (intensity ignored; no total spectrum column).
 
     Returns
     -------
     SpectrumConfig
         Configuration ready for BetaSpectrum.from_config().
     """
-    # Convert BranchInfo -> BranchConfig for SpectrumConfig
+    # Convert BranchInfo -> BranchConfig for SpectrumConfig. Per-branch
+    # endpoints are the measured transition endpoint energies when
+    # available; the base-endpoint minus level-energy subtraction is only
+    # a legacy fallback for hand-built DecayInfo objects without them.
     branches: Optional[List[BranchConfig]] = None
     if info.branches:
         branches = [
             BranchConfig(
-                endpoint_MeV=info.endpoint_MeV - b.level_energy_keV / 1000.0,
+                endpoint_MeV=(
+                    b.endpoint_MeV
+                    if b.endpoint_MeV > 0.0
+                    else info.endpoint_MeV - b.level_energy_keV / 1000.0
+                ),
                 transition_type=b.transition_type,
                 intensity=b.intensity,
             )
@@ -523,7 +536,7 @@ def decay_info_to_config(
         branches=branches,
         intensity_cutoff=intensity_cutoff,
         csv_separator=csv_separator,
-        normalize_to_one=normalize_to_one,
+        separate_branches=separate_branches,
     )
 
 
@@ -542,14 +555,6 @@ DEFAULT_JSON_SCHEMA = {
     "transition_type": "F2",
     "e_step_MeV": 0.001,
     "neutrino_mass_MeV": 0.0,
-    "detector_model": "gaussian",
-    "detector_sigma_a_keV": 1.0,
-    "detector_sigma_b": 0.0,
-    "detector_tail_fraction": 0.0,
-    "detector_tau_keV": 5.0,
-    "detector_fano_factor": 0.12,
-    "detector_n_channels": 4096,
-    "detector_channel_energy_range": [0.0, 0.35],
 }
 
 
@@ -615,8 +620,13 @@ def load_json_input(filepath: str) -> Dict[str, Any]:
     if not path.exists():
         raise FileNotFoundError(f"Input file not found: {filepath}")
 
-    with open(path, "r") as f:
+    try:
+        with open(path, "r", encoding="utf-8") as f:
             data = json.load(f)
+    except OSError as exc:
+        raise FileNotFoundError(f"Input file not readable: {filepath}") from exc
+    except json.JSONDecodeError as exc:
+        raise ValueError(f"Invalid JSON in input file {filepath}: {exc}") from exc
 
     # Fill in defaults for missing optional fields
     for key, default in DEFAULT_JSON_SCHEMA.items():
@@ -699,7 +709,7 @@ def create_config_from_source(
         intensity_cutoff = kwargs.get("intensity_cutoff", 0.0)
         neutrino_mass = kwargs.get("neutrino_mass_MeV", 0.0)
         csv_separator = kwargs.get("csv_separator", ",")
-        normalize_to_one = kwargs.get("normalize_to_one", False)
+        separate_branches = kwargs.get("separate_branches", False)
         info = get_decay_info_from_paceENSDF(nuclide, decay_type, decay_index)
         return decay_info_to_config(
             info,
@@ -707,7 +717,7 @@ def create_config_from_source(
             intensity_cutoff=intensity_cutoff,
             neutrino_mass_MeV=neutrino_mass,
             csv_separator=csv_separator,
-            normalize_to_one=normalize_to_one,
+            separate_branches=separate_branches,
         )
     elif source == "json":
         if not json_path:

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

beta_spectrum/output/core.py

@@ -14,6 +14,8 @@ from typing import TYPE_CHECKING, Any, Dict, List, Optional
 
 import numpy as np
 
+from beta_spectrum.components.phase_space import PhaseSpace
+from beta_spectrum.components.radiative import RadiativeCorrection
 from beta_spectrum.utils import T_to_W
 
 if TYPE_CHECKING:
@@ -47,17 +49,26 @@ class BetaSpectrumAnalyzer:
         self.config = config
         self.W, self.energies_MeV = spectrum.get_energy_grid(config)
         self._components_cache: Optional[Dict[str, np.ndarray]] = None
-        self._branch_spectra_cache: Optional[List[np.ndarray]] = None
+        self._universal_cache: Optional[np.ndarray] = None
+        self._physical_branch_cache: Dict[int, np.ndarray] = {}
         self._logger = logger
         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,
                 )
             ]
@@ -82,15 +93,76 @@ class BetaSpectrumAnalyzer:
         assert self._components_cache is not None
         return self._components_cache
 
-    @property
-    def branch_spectra(self) -> List[np.ndarray]:
-        """Intensity-weighted branch spectra."""
-        if self._branch_spectra_cache is None:
-            self._branch_spectra_cache = self.spectrum.get_branch_normalized_spectra(
-                self.W
+    def universal_factor(self) -> np.ndarray:
+        """Product of the components shared by all branches.
+
+        The universal factor always contains the Fermi function (the core
+        0-order Coulomb correction) and every enabled higher-order
+        correction that does not depend on the branch endpoint
+        (screening, finite size, charge distribution, exchange).
+
+        In single-branch mode this equals the product of all components
+        except the per-branch phase space and radiative corrections.
+        """
+        if self._universal_cache is None:
+            prod: Optional[np.ndarray] = None
+            for comp in self.spectrum.components:
+                if not self._is_multi and isinstance(
+                    comp, (PhaseSpace, RadiativeCorrection)
+                ):
+                    continue
+                values = comp(self.W)
+                prod = (
+                    np.array(values, dtype=np.float64)
+                    if prod is None
+                    else prod * values
                 )
-        assert self._branch_spectra_cache is not None
-        return self._branch_spectra_cache
+            if prod is None:
+                prod = np.ones_like(self.W, dtype=np.float64)
+            self._universal_cache = prod
+        assert self._universal_cache is not None
+        return self._universal_cache
+
+    def physical_branch_spectrum(self, i: int) -> np.ndarray:
+        """Physical spectrum of branch i: core × corrections, unnormalized.
+
+        The physical branch spectrum is
+            core(phase_space × fermi) × all enabled higher-order corrections,
+        masked to zero beyond the branch endpoint. This is the canonical
+        branch accessor: every consumer (CSV export, plots) must use it so
+        that the Fermi function is always present in branch-level spectra.
+        """
+        if i in self._physical_branch_cache:
+            return self._physical_branch_cache[i]
+        if i < len(self.spectrum.branch_spectra):
+            raw = self.spectrum.branch_spectra[i](self.W)
+            W0 = T_to_W(self._unified_branches[i].endpoint_MeV)
+            raw = np.where(self.W <= W0, raw, 0.0)
+            physical = self.universal_factor() * raw
+        else:
+            # Single-branch mode: the spectrum itself is the full physical branch
+            physical = np.asarray(self.spectrum(self.W), dtype=np.float64)
+        self._physical_branch_cache[i] = physical
+        return physical
+
+    def total_integral(self) -> float:
+        """Integral of the raw total spectrum over kinetic energy.
+
+        Raises
+        ------
+        ValueError
+            If the integral is not positive (e.g. all branches are
+            degenerate because the neutrino mass exceeds every endpoint).
+        """
+        total_raw = self.spectrum(self.W)
+        integral = float(np.trapezoid(total_raw, self.energies_MeV))
+        if integral <= 0.0:
+            raise ValueError(
+                "Total spectrum integral is non-positive: every branch is "
+                "degenerate (check neutrino_mass vs branch endpoints and "
+                "intensity_cutoff)"
+            )
+        return integral
 
     def total_spectrum(self, normalize: bool = True) -> np.ndarray:
         """Calculate the total spectrum."""
@@ -98,6 +170,12 @@ class BetaSpectrumAnalyzer:
 
         if normalize:
             integral = np.trapezoid(total, self.energies_MeV)
+            if integral <= 0.0:
+                raise ValueError(
+                    "Total spectrum integral is non-positive: every branch "
+                    "is degenerate (check neutrino_mass vs branch endpoints "
+                    "and intensity_cutoff)"
+                )
             total = total / integral
 
         return total
@@ -111,36 +189,12 @@ class BetaSpectrumAnalyzer:
         """Number of branches (unified view)."""
         return len(self._unified_branches)
 
-    def _get_branch_spectrum(self, i: int) -> np.ndarray:
-        """Get the spectrum array for branch i."""
-        if i < len(self.spectrum.branch_spectra):
-            bs_obj = self.spectrum.branch_spectra[i]
-            # branch_spectra stores BetaSpectrum objects — call them to get arrays
-            if hasattr(bs_obj, "__call__"):
-                val = bs_obj(self.W)
-                # Clamp to zero beyond endpoint (phase space goes negative)
-                W0 = T_to_W(self._unified_branches[i].endpoint_MeV)
-                return np.where(self.W <= W0, val, 0.0)
-            return np.asarray(bs_obj)
-        # Single-branch: the spectrum itself is the only branch
-        return self.spectrum(self.W)
-
     def _get_component_key(self, i: int, name: str) -> str:
         """Get the canonical key for a per-branch component."""
         if self._is_multi:
             return f"branch_{i}.{name}"
         return name
 
-    def _has_component(self, i: int, name: str) -> bool:
-        """Check if branch i has a given component."""
-        key = self._get_component_key(i, name)
-        return key in self.components
-
-    def _get_component(self, i: int, name: str) -> np.ndarray:
-        """Get a per-branch component array."""
-        key = self._get_component_key(i, name)
-        return self.components[key]
-
     def _element_symbol(self, Z: int) -> str:
         """Return chemical symbol for atomic number Z."""
         symbols = [

beta_spectrum/output/csv_export.py

@@ -1,14 +1,28 @@
 """
 CSV export module — CSV export functionality for BetaSpectrumAnalyzer.
+
+Column semantics (v0.4.1 onward):
+- The basic spectrum is ALWAYS core = phase_space x fermi; the Fermi function
+  is a core spectrum component, never an optional correction. Higher-order
+  corrections (finite size, charge distribution, screening, exchange,
+  radiative) are multiplied on top of the core.
+- ``branch_{i}_spectrum`` columns carry the PHYSICAL branch spectrum
+  (core x all enabled corrections), weighted by intensity and
+  area-normalized so that the sum of branch columns equals total_spectrum
+  (the total column itself is normalized to integrate to 1). In the legacy
+  export branch columns contained only the per-branch phase-space x
+  radiative product WITHOUT the Fermi function — physically wrong.
+- ``separate_branches`` mode: each branch column is normalized to
+  integrate to 1 individually, intensity is ignored, no total column.
 """
 
 from __future__ import annotations
 
-import pandas as pd
 from datetime import datetime, timezone
 from typing import TYPE_CHECKING, Dict
 
 import numpy as np
+import pandas as pd
 
 from beta_spectrum.logging_utils import get_git_short_hash
 
@@ -37,64 +51,92 @@ class CSVEmitter:
         source_type : str
             Source type (paceENSDF, json, etc.).
         csv_verbosity : str
-            Verbosity level: 'full' (all components) or 'concise' (total spectrum and normalized branch components only).
-        """
-        total = analyzer.total_spectrum(normalize=True)
-        components = analyzer.components
+            Verbosity level: 'full' (all components) or 'concise' (total
+            spectrum and physical branch spectra only).
 
-        # Only include enabled corrections that are not basic components
-        # Basic components (phase_space and fermi) are always used and should not be listed
+        Raises
+        ------
+        ValueError
+            If every branch spectrum integrates to a non-positive area
+            (e.g. all branch endpoints lie below the neutrino mass).
+        """
+        components = analyzer.components
+        config = analyzer.config
+        separate = config.separate_branches
+        n_branches = analyzer._n_branches
+
+        # Physical branch spectra: core(phase_space x fermi) x corrections.
+        physical = [analyzer.physical_branch_spectrum(i) for i in range(n_branches)]
+
+        # Only higher-order corrections are toggleable; list them enabled.
         enabled = []
-        if analyzer.config.use_screening:
+        if config.use_screening:
             enabled.append("screening")
-        if analyzer.config.use_finite_size:
+        if config.use_finite_size:
             enabled.append("finite_size")
-        if analyzer.config.use_charge_dist:
+        if config.use_charge_dist:
             enabled.append("charge_dist")
-        if analyzer.config.use_radiative:
+        if config.use_radiative:
             enabled.append("radiative")
-        if analyzer.config.use_exchange:
+        if config.use_exchange:
             enabled.append("exchange")
 
-        parent_symbol = analyzer._element_symbol(analyzer.config.Z_parent)
-        daughter_symbol = analyzer._element_symbol(analyzer.config.Z_daughter)
+        # 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)
 
-        # Get git commit hash
         git_commit = get_git_short_hash()
 
-        # Format source line with paceENSDF URL if applicable
         if source_type == "paceENSDF":
             source_line = "# source: paceENSDF (https://pypi.org/project/paceENSDF/)"
         else:
             source_line = f"# source: {source_type}"
 
         header_lines = [
-            f"# beta-spectrum v{__import__('beta_spectrum').__version__}(build.{git_commit})",
+            f"# beta-spectrum v{__import__('beta_spectrum').__version__}"
+            f"(build.{git_commit})",
             f"# timestamp: {datetime.now(timezone.utc).isoformat()}",
             source_line,
-            f"# nuclide: {parent_symbol}{analyzer.config.A_number} -> {daughter_symbol}{analyzer.config.A_number} (Z={analyzer.config.Z_parent}->{analyzer.config.Z_daughter})",
-            f"# basic_endpoint: {analyzer.config.endpoint_MeV * 1000:.1f} keV",
+            f"# nuclide: {parent_symbol}{config.A_number} "
+            f"-> {daughter_symbol}{config.A_number} "
+            f"(Z={config.Z_parent}->{config.Z_daughter})",
+            f"# basic_endpoint: {basic_endpoint_MeV * 1000:.1f} keV",
         ]
 
-        # Add neutrino mass info if non-zero
-        if analyzer.config.neutrino_mass_MeV > 0:
+        if config.neutrino_mass_MeV > 0:
             header_lines.append(
-                f"# neutrino_mass_MeV: {analyzer.config.neutrino_mass_MeV:.6f} (non-zero neutrino mass effects included)"
+                f"# neutrino_mass_MeV: {config.neutrino_mass_MeV:.6f} "
+                f"(non-zero neutrino mass effects included)"
             )
         else:
-            header_lines.append(
-                f"# neutrino_mass_MeV: {analyzer.config.neutrino_mass_MeV:.6f}"
-            )
+            header_lines.append(f"# neutrino_mass_MeV: {config.neutrino_mass_MeV:.6f}")
 
+        header_lines.append("# core: phase_space x fermi (basic spectrum)")
         header_lines.append(f"# corrections: {', '.join(enabled)}")
-        header_lines.append(f"# e_step: {analyzer.config.e_step_MeV:.4f} MeV")
+        header_lines.append(f"# e_step: {config.e_step_MeV:.4f} MeV")
         header_lines.append(f"# n_points: {len(analyzer.energies_MeV)}")
 
-        # Add normalization mode
-        if analyzer.config.normalize_to_one:
-            header_lines.append("# branch normalization: to_one (ignoring intensity)")
+        if separate:
+            header_lines.append(
+                "# branch normalization: separate "
+                "(each branch integrates to 1, intensity ignored)"
+            )
         else:
-            header_lines.append("# branch normalization: intensity")
+            header_lines.append(
+                "# branch normalization: intensity "
+                "(sum of branch columns == total_spectrum)"
+            )
 
         # Separate branches into active and excluded lists
         active_branches = []
@@ -102,118 +144,99 @@ class CSVEmitter:
 
         for i, branch in enumerate(analyzer._unified_branches):
             branch_endpoint_keV = branch.endpoint_MeV * 1000
-            neutrino_mass_keV = analyzer.config.neutrino_mass_MeV * 1000
+            neutrino_mass_keV = config.neutrino_mass_MeV * 1000
 
-            # Check if branch is excluded due to neutrino mass
             if branch_endpoint_keV <= neutrino_mass_keV:
                 excluded_branches.append(
-                    f"  branch_{i + 1}: E₀={branch_endpoint_keV:.1f} keV, transition={branch.transition_type}, intensity={branch.intensity:.4f} - excluded: endpoint < m_nu"
+                    f"  branch_{i + 1}: E₀={branch_endpoint_keV:.1f} keV, "
+                    f"transition={branch.transition_type}, "
+                    f"intensity={branch.intensity:.4f} - excluded: endpoint < m_nu"
                 )
-            # Check if branch is excluded due to intensity cutoff
-            elif branch.intensity < analyzer.config.intensity_cutoff:
+            elif branch.intensity < config.intensity_cutoff:
                 excluded_branches.append(
-                    f"  branch_{i + 1}: E₀={branch_endpoint_keV:.1f} keV, transition={branch.transition_type}, intensity={branch.intensity:.4f} - excluded: intensity < cutoff"
+                    f"  branch_{i + 1}: E₀={branch_endpoint_keV:.1f} keV, "
+                    f"transition={branch.transition_type}, "
+                    f"intensity={branch.intensity:.4f} "
+                    "- excluded: intensity < cutoff"
                 )
             else:
                 active_branches.append(
-                    f"  branch_{i + 1}: E₀={branch_endpoint_keV:.1f} keV, transition={branch.transition_type}, intensity={branch.intensity:.4f}"
+                    f"  branch_{i + 1}: E₀={branch_endpoint_keV:.1f} keV, "
+                    f"transition={branch.transition_type}, "
+                    f"intensity={branch.intensity:.4f}"
                 )
 
-        # Add active branches list
         if active_branches:
             header_lines.append("# active branches:")
-            for branch_info in active_branches:
-                header_lines.append(f"# {branch_info}")
+            header_lines.extend(f"# {b}" for b in active_branches)
 
-        # Add excluded branches list if there are any
         if excluded_branches:
             header_lines.append("# excluded branches:")
-            for branch_info in excluded_branches:
-                header_lines.append(f"# {branch_info}")
+            header_lines.extend(f"# {b}" for b in excluded_branches)
 
-        # Create data dictionary with concise column headers
+        # ----- data columns -----
         data: Dict[str, np.ndarray] = {
             "energy_MeV": analyzer.energies_MeV,
         }
 
-        # Only add total_spectrum if not normalizing to one
-        if not analyzer.config.normalize_to_one:
-            data["total_spectrum"] = total
+        # Positive-integral guard applies to BOTH normalization modes: a
+        # degenerate spectrum must never be written (not even as zeros).
+        total_raw = analyzer.spectrum(analyzer.W)
+        integral = float(np.trapezoid(total_raw, analyzer.energies_MeV))
+        if integral <= 0.0:
+            raise ValueError(
+                "Total spectrum integral is non-positive: every branch "
+                "is degenerate (check neutrino_mass vs branch endpoints "
+                "and intensity_cutoff)"
+            )
+        if not separate:
+            data["total_spectrum"] = total_raw / integral
+
+        for i in range(n_branches):
+            key_i = i + 1
+            if csv_verbosity == "concise":
+                col_name = f"branch_{key_i}"
+            else:
+                col_name = f"branch_{key_i}_spectrum"
+
+            if separate:
+                area = float(np.trapezoid(physical[i], analyzer.energies_MeV))
+                if area > 0.0:
+                    data[col_name] = physical[i] / area
+                else:
+                    data[col_name] = np.zeros_like(physical[i])
+            else:
+                intensity_i = analyzer._unified_branches[i].intensity
+                data[col_name] = physical[i] * intensity_i / integral
 
         if csv_verbosity == "full":
-            # Per-branch spectra and components with concise branch naming
-            for i in range(analyzer._n_branches):
-                # Get branch spectrum
-                branch_spectrum = analyzer._get_branch_spectrum(i)
-
-                if analyzer.config.normalize_to_one:
-                    # Normalize to 1, ignoring intensity
-                    sum_spectrum = np.sum(branch_spectrum)
-                    if sum_spectrum > 0:
-                        normalized_spectrum = branch_spectrum / sum_spectrum
-                    else:
-                        normalized_spectrum = branch_spectrum
-                else:
-                    # Use original spectrum with intensity
-                    normalized_spectrum = (
-                        branch_spectrum * analyzer._unified_branches[i].intensity
-                    )
-
-                data[f"branch_{i + 1}_spectrum"] = normalized_spectrum
-
+            # Per-branch component columns (factorized ingredients).
+            # Only multi-branch mode carries per-branch factors: in single
+            # mode the bare PhaseSpace/Radiative columns below already
+            # describe the single spectrum, and prefixed duplicates would
+            # silently disagree with them once renamed.
+            if analyzer._is_multi:
+                for i in range(n_branches):
                     for comp_name in ["PhaseSpace", "Radiative"]:
                         key = analyzer._get_component_key(i, comp_name)
                         if key in components:
                             data[f"branch_{i + 1}_{comp_name}"] = components[key]
 
-            # Universal components (no prefix)
+            # Universal components (no prefix): Fermi plus universal
+            # higher-order corrections. In single mode the core factors
+            # (PhaseSpace, Radiative) also appear here unprefixed.
             for name, values in components.items():
                 if not name.startswith("branch_"):
                     data[name] = values
-        elif csv_verbosity == "concise":
-            # Only total spectrum and normalized branch components with concise naming
-            for i in range(analyzer._n_branches):
-                branch_spectrum = analyzer._get_branch_spectrum(i)
-                # Check if branch is degenerate (neutrino mass >= branch endpoint)
-                branch_endpoint_keV = analyzer._unified_branches[i].endpoint_MeV * 1000
-                neutrino_mass_keV = analyzer.config.neutrino_mass_MeV * 1000
-                if branch_endpoint_keV <= neutrino_mass_keV:
-                    # Branch is degenerate, spectrum should be all zeros
-                    data[f"branch_{i + 1}"] = np.zeros_like(branch_spectrum)
-                else:
-                    # Branch is valid, normalize as usual
-                    if analyzer.config.normalize_to_one:
-                        # Normalize to 1, ignoring intensity
-                        sum_spectrum = np.sum(branch_spectrum)
-                        if sum_spectrum > 0:
-                            data[f"branch_{i + 1}"] = branch_spectrum / sum_spectrum
-                        else:
-                            data[f"branch_{i + 1}"] = branch_spectrum
-                    else:
-                        # Normalize as usual with intensity
-                        sum_spectrum = np.sum(branch_spectrum)
-                        if sum_spectrum > 0:
-                            data[f"branch_{i + 1}"] = (
-                                branch_spectrum
-                                / sum_spectrum
-                                * analyzer._unified_branches[i].intensity
-                            )
-                        else:
-                            data[f"branch_{i + 1}"] = branch_spectrum
 
         df = pd.DataFrame(data)
-        # Determine separator
-        separator = (
-            "\t"
-            if analyzer.config.csv_separator == "tab"
-            else analyzer.config.csv_separator
-        )
-        df.to_csv(filename, index=False, float_format="%.4e", sep=separator)
-
-        with open(filename, "r") as f:
-            content = f.read()
+        separator = "\t" if config.csv_separator == "tab" else config.csv_separator
+        try:
             with open(filename, "w") as f:
-            f.write("\n".join(header_lines) + "\n" + content)
+                f.write("\n".join(header_lines) + "\n")
+                df.to_csv(f, index=False, float_format="%.4e", sep=separator)
+        except OSError as exc:
+            raise ValueError(f"Cannot write spectrum CSV {filename}: {exc}") from exc
 
         if analyzer._logger:
             analyzer._logger.info("CSV exported to %s", filename)

beta_spectrum/output/visualization.py

@@ -1,68 +1,195 @@
 """
-Visualization module — Plotting functionality for BetaSpectrumAnalyzer.
+Visualization module — plotting built ON TOP OF the CSV export.
+
+All graphical representations are produced from the exported CSV file
+(data columns + metadata header) and never recomputed from the analyzer in
+parallel. The CSV file is the single source of truth: a plot is produced
+exactly from the bytes a researcher would download, which turns every plot
+into an additional quality check of the export routine.
 """
 
 from __future__ import annotations
 
-from datetime import datetime, timezone
-from typing import TYPE_CHECKING, Dict, Optional
+import re
+import shutil
+import tempfile
+from typing import TYPE_CHECKING, Any, Dict, List, Optional, Tuple
 
 import numpy as np
-import matplotlib.pyplot as plt
+import pandas as pd
 
-from beta_spectrum.logging_utils import get_git_short_hash
-from beta_spectrum.utils import T_to_W
+import matplotlib.pyplot as plt
+from matplotlib.figure import Figure
+
+from .csv_export import CSVEmitter
 
 if TYPE_CHECKING:
-    from matplotlib.figure import Figure
-
     from beta_spectrum.output.core import BetaSpectrumAnalyzer
 
+_BRANCH_LINE_RE = re.compile(
+    r"^\s*branch_(\d+):\s*E₀=([0-9.]+)\s*keV,\s*transition=(\w+),\s*"
+    r"intensity=([0-9.]+)(.*)$"
+)
+_EXCLUDED_SUFFIX_RE = re.compile(r"-\s*excluded:\s*(.*)$")
 
-class Visualizer:
-    """Visualization functionality for BetaSpectrumAnalyzer."""
+_BRANCH_COLORS = [
+    "#1f77b4",
+    "#ff7f0e",
+    "#2ca02c",
+    "#d62728",
+    "#9467bd",
+    "#8c564b",
+    "#e377c2",
+    "#7f7f7f",
+    "#bcbd22",
+    "#17becf",
+]
 
-    @staticmethod
-    def plot_analysis(
-        analyzer: BetaSpectrumAnalyzer,
-        save_path: Optional[str] = None,
-        show_components: bool = True,
-    ) -> None:
-        """
-        Create visualization of the spectrum and all correction factors.
+
+def parse_spectrum_csv(
+    csv_path: str,
+) -> Tuple[pd.DataFrame, Dict[str, Any]]:
+    """Read an exported spectrum CSV and parse its metadata header.
 
     Parameters
     ----------
-        analyzer : BetaSpectrumAnalyzer
-            The analyzer instance to visualize.
-        save_path : str, optional
-            Path to save the figure.
-        show_components : bool
-            If True, show full debug view.
-            If False, show only the spectrum plot.
+    csv_path : str
+        Path to a CSV file produced by CSVEmitter.export_to_csv.
+
+    Returns
+    -------
+    tuple[pandas.DataFrame, dict]
+        The data frame (comment lines skipped, separator auto-detected)
+        and a metadata dict with keys: version, build, timestamp, source,
+        nuclide, endpoint_keV, neutrino_mass_MeV, corrections, e_step_MeV,
+        n_points, normalization, branches.
+
+    Raises
+    ------
+    ValueError
+        If the file cannot be read or parsed.
     """
-        total = analyzer.total_spectrum(normalize=True)
-        components = analyzer.components
-        analyzer._timestamp = datetime.now(timezone.utc).strftime(
-            "%Y-%m-%d %H:%M:%S UTC"
-        )
-
-        if show_components:
-            Visualizer._plot_debug_view(analyzer, total, components, save_path)
+    try:
+        df = pd.read_csv(csv_path, sep=None, engine="python", comment="#")
+        raw_lines = []
+        with open(csv_path, "r") as f:
+            for line in f:
+                if line.startswith("#"):
+                    raw_lines.append(line.rstrip("\n"))
                 else:
-            Visualizer._plot_spectra_view(analyzer, total, components, save_path)
+                    break  # header block ended
+    except OSError as exc:
+        raise ValueError(f"Cannot read spectrum CSV {csv_path}: {exc}") from exc
 
-    @staticmethod
-    def _build_figure_header(analyzer: BetaSpectrumAnalyzer, fig: "Figure") -> None:
-        """Add a compact header with nuclide info, e_step, commit and timestamp outside axes."""
-        parent = analyzer._element_symbol(analyzer.config.Z_parent)
-        daughter = analyzer._element_symbol(analyzer.config.Z_daughter)
-        e_step_keV = analyzer.config.e_step_MeV * 1000
-        header = (
-            f"{parent}{analyzer.config.A_number} -> {daughter}{analyzer.config.A_number}  "
-            f"|  E₀={analyzer.config.endpoint_MeV * 1000:.1f} keV  "
-            f"|  ΔE={e_step_keV:.1f} keV"
+    meta: Dict[str, Any] = {
+        "version": None,
+        "build": None,
+        "timestamp": None,
+        "source": None,
+        "nuclide": None,
+        "endpoint_keV": None,
+        "neutrino_mass_MeV": None,
+        "corrections": [],
+        "e_step_MeV": None,
+        "n_points": None,
+        "normalization": None,
+        "branches": [],
+    }
+
+    for line in raw_lines:
+        body = line.lstrip("#").strip()
+        m = re.match(r"beta-spectrum v([\d.]+)\(build\.(\w+)\)", body)
+        if m:
+            meta["version"] = m.group(1)
+            meta["build"] = m.group(2)
+            continue
+        m = re.match(r"timestamp:\s*(.+)", body)
+        if m:
+            meta["timestamp"] = m.group(1).strip()
+            continue
+        m = re.match(r"source:\s*(\S+)", body)
+        if m:
+            meta["source"] = m.group(1)
+            continue
+        m = re.match(r"nuclide:\s*(.+?)\s+\(Z=", body)
+        if m:
+            meta["nuclide"] = m.group(1).strip()
+            continue
+        m = re.match(r"basic_endpoint:\s*([0-9.]+)\s*keV", body)
+        if m:
+            meta["endpoint_keV"] = float(m.group(1))
+            continue
+        m = re.match(r"neutrino_mass_MeV:\s*([0-9.]+)", body)
+        if m:
+            meta["neutrino_mass_MeV"] = float(m.group(1))
+            continue
+        m = re.match(r"corrections:\s*(.*)", body)
+        if m:
+            meta["corrections"] = [
+                c.strip() for c in m.group(1).split(",") if c.strip()
+            ]
+            continue
+        m = re.match(r"e_step:\s*([0-9.]+)\s*MeV", body)
+        if m:
+            meta["e_step_MeV"] = float(m.group(1))
+            continue
+        m = re.match(r"n_points:\s*(\d+)", body)
+        if m:
+            meta["n_points"] = int(m.group(1))
+            continue
+        m = re.match(r"branch normalization:\s*(\S+)", body)
+        if m:
+            meta["normalization"] = m.group(1)
+            continue
+        m = _BRANCH_LINE_RE.match(body)
+        if m:
+            excluded = _EXCLUDED_SUFFIX_RE.search(m.group(5))
+            meta["branches"].append(
+                {
+                    "index": int(m.group(1)),
+                    "endpoint_keV": float(m.group(2)),
+                    "transition": m.group(3),
+                    "intensity": float(m.group(4)),
+                    "excluded": None if excluded is None else excluded.group(1).strip(),
+                }
             )
+            continue
+
+    if meta["n_points"] is None:
+        meta["n_points"] = len(df)
+
+    return df, meta
+
+
+def _branch_column(df: pd.DataFrame, index: int) -> Optional[str]:
+    """Find the branch spectrum column name for a 1-based branch index."""
+    for candidate in (f"branch_{index}_spectrum", f"branch_{index}"):
+        if candidate in df.columns:
+            return candidate
+    return None
+
+
+def _active_branches(meta: Dict[str, Any]) -> List[Dict[str, Any]]:
+    return [b for b in meta["branches"] if b["excluded"] is None]
+
+
+def _spectrum_title(meta: Dict[str, Any]) -> str:
+    n_branches = len(meta["branches"])
+    active = len(_active_branches(meta))
+    title = f"Beta-decay: {meta['nuclide']}"
+    if n_branches > 1:
+        if active < n_branches:
+            title += f" ({active} of {n_branches} branches)"
+        else:
+            title += f" ({n_branches} branches)"
+    return title
+
+
+def _apply_figure_header(fig: Figure, meta: Dict[str, Any]) -> None:
+    """Write the figure header/footer lines from the CSV metadata."""
+    e_step_keV = (meta["e_step_MeV"] or 0.0) * 1000.0
+    endpoint = meta["endpoint_keV"]
+    header = f"{meta['nuclide']}  |  E₀={endpoint:.1f} keV  |  ΔE={e_step_keV:.1f} keV"
     fig.text(
         0.01,
         0.98,
@@ -74,13 +201,11 @@ class Visualizer:
         fontfamily="monospace",
     )
 
-        # Add neutrino mass info to header if non-zero
-        if analyzer.config.neutrino_mass_MeV > 0:
-            nu_mass_text = f"|  m_ν={analyzer.config.neutrino_mass_MeV:.6f} MeV"
+    if meta["neutrino_mass_MeV"]:
         fig.text(
             0.01,
             0.96,
-                nu_mass_text,
+            f"|  m_ν={meta['neutrino_mass_MeV']:.6f} MeV",
             transform=fig.transFigure,
             fontsize=8,
             verticalalignment="top",
@@ -89,11 +214,11 @@ class Visualizer:
             color="red",
         )
 
-        id_text = f"{get_git_short_hash(6)}  |  {analyzer._timestamp}"
+    if meta["build"]:
         fig.text(
             0.01,
             0.01,
-            id_text,
+            f"{meta['build']}  |  {meta['timestamp']}",
             transform=fig.transFigure,
             fontsize=8,
             verticalalignment="bottom",
@@ -102,309 +227,141 @@ class Visualizer:
             color="gray",
         )
 
-    @staticmethod
-    def _build_figure_title(analyzer: BetaSpectrumAnalyzer, fig: "Figure") -> str:
-        """Build a unified figure title."""
-        parent = analyzer._element_symbol(analyzer.config.Z_parent)
-        daughter = analyzer._element_symbol(analyzer.config.Z_daughter)
-        n = analyzer.config.A_number
-        title = f"Beta-decay: {n}{parent} -> {n}{daughter}"
-        if analyzer._n_branches > 1:
-            # Count active branches (not excluded by neutrino mass or intensity cutoff)
-            active_branches = 0
-            for i, branch in enumerate(analyzer._unified_branches):
-                branch_endpoint_keV = branch.endpoint_MeV * 1000
-                neutrino_mass_keV = analyzer.config.neutrino_mass_MeV * 1000
-                is_excluded_by_mass = branch_endpoint_keV <= neutrino_mass_keV
-                is_excluded_by_intensity = (
-                    branch.intensity < analyzer.config.intensity_cutoff
-                )
-                if not is_excluded_by_mass and not is_excluded_by_intensity:
-                    active_branches += 1
 
-            if active_branches < analyzer._n_branches:
-                title += f" ({active_branches} of {analyzer._n_branches} branches)"
-            else:
-                title += f" ({analyzer._n_branches} branches)"
-        return title
-
-    @staticmethod
-    def _format_branch_label(analyzer: BetaSpectrumAnalyzer, i: int) -> str:
-        """Abbreviated branch label for legend."""
-        branch = analyzer._unified_branches[i]
-        e0 = branch.endpoint_MeV * 1000
-        intensity_pct = branch.intensity * 100
-
-        # Check if branch is excluded
-        branch_endpoint_keV = branch.endpoint_MeV * 1000
-        neutrino_mass_keV = analyzer.config.neutrino_mass_MeV * 1000
-        is_excluded_by_mass = branch_endpoint_keV <= neutrino_mass_keV
-        is_excluded_by_intensity = branch.intensity < analyzer.config.intensity_cutoff
-
-        if is_excluded_by_mass:
-            return f"Br. {i + 1}: {e0:.1f} keV ({intensity_pct:.1f}%) - excluded (endpoint < m_nu)"
-        elif is_excluded_by_intensity:
-            return f"Br. {i + 1}: {e0:.1f} keV ({intensity_pct:.1f}%) - excluded (intensity < cutoff)"
-        else:
-            return f"Br. {i + 1}: {e0:.1f} keV ({intensity_pct:.1f}%)"
-
-    @staticmethod
-    def _plot_spectra_view(
-        analyzer: BetaSpectrumAnalyzer,
-        total: np.ndarray,
-        components: Dict[str, np.ndarray],
+def _plot_spectra_view_from_df(
+    df: pd.DataFrame,
+    meta: Dict[str, Any],
     save_path: Optional[str],
-    ) -> None:
-        """Plot total spectrum with individual branch decomposition."""
+) -> Figure:
+    """Plot total spectrum and per-branch spectra from CSV columns."""
+    energies = df["energy_MeV"].to_numpy()
+
     fig, ax = plt.subplots(figsize=(12, 7))
 
-        branch_colors = [
-            "#1f77b4",
-            "#ff7f0e",
-            "#2ca02c",
-            "#d62728",
-            "#9467bd",
-            "#8c564b",
-            "#e377c2",
-            "#7f7f7f",
-            "#bcbd22",
-            "#17becf",
-        ]
-
-        e_step_keV = analyzer.config.e_step_MeV * 1000
-
-        # Compute branch contributions (universal * intensity * branch_spectrum)
-        branch_contribs = []
-        # Universal components are shared (no per-branch key prefix)
-        universal_names = [
-            "Fermi",
-            "Screening",
-            "FiniteSizeL0",
-            "ChargeDistributionU",
-            "Exchange",
-        ]
-        for i in range(analyzer._n_branches):
-            universal_product = np.ones_like(analyzer.W)
-            for name in universal_names:
-                key = analyzer._get_component_key(i, name)
-                if key in components:
-                    universal_product *= components[key]
-                elif name in components:  # fallback: universal components
-                    universal_product *= components[name]
-            bs = analyzer._get_branch_spectrum(i)
-            contrib = universal_product * analyzer._unified_branches[i].intensity * bs
-            branch_contribs.append(contrib)
-
-        total_raw = sum(branch_contribs)
-
+    if "total_spectrum" in df.columns:
         ax.plot(
-            analyzer.energies_MeV,
-            total_raw,
+            energies,
+            df["total_spectrum"].to_numpy(),
             "k-",
             lw=2.5,
             label="Total spectrum",
             zorder=10,
         )
 
-        # Only plot active branches (not excluded by neutrino mass or intensity cutoff)
-        active_branch_indices = []
-        for i in range(analyzer._n_branches):
-            branch_endpoint_keV = analyzer._unified_branches[i].endpoint_MeV * 1000
-            neutrino_mass_keV = analyzer.config.neutrino_mass_MeV * 1000
-            is_excluded_by_mass = branch_endpoint_keV <= neutrino_mass_keV
-            is_excluded_by_intensity = (
-                analyzer._unified_branches[i].intensity
-                < analyzer.config.intensity_cutoff
-            )
-            if not is_excluded_by_mass and not is_excluded_by_intensity:
-                active_branch_indices.append(i)
-
-        for i in active_branch_indices:
-            contrib = branch_contribs[i]
-            color = branch_colors[i % len(branch_colors)]
+    for color_i, b in enumerate(_active_branches(meta)):
+        col = _branch_column(df, int(b["index"]))
+        if col is None:
+            continue
         ax.plot(
-                analyzer.energies_MeV,
-                contrib,
-                color=color,
+            energies,
+            df[col].to_numpy(),
+            color=_BRANCH_COLORS[color_i % len(_BRANCH_COLORS)],
             lw=1.5,
             alpha=0.7,
-                label=Visualizer._format_branch_label(analyzer, i),
+            label=(
+                f"Br. {b['index']}: {b['endpoint_keV']:.1f} keV "
+                f"({b['intensity'] * 100:.1f}%)"
+            ),
         )
 
     ax.set_xlabel(r"Electron kinetic energy $E$ [MeV]", fontsize=12)
-        ax.set_ylabel(f"Normalized Counts / {e_step_keV:.1f} keV", fontsize=12)
-        fig.suptitle(
-            Visualizer._build_figure_title(analyzer, fig),
-            fontsize=14,
-            fontweight="bold",
-        )
+    ax.set_ylabel(r"Normalized spectrum [MeV$^{-1}$]", fontsize=12)
+    fig.suptitle(_spectrum_title(meta), fontsize=14, fontweight="bold")
     ax.set_yscale("log")
     ax.grid(True, alpha=0.3)
-        ax.legend(
-            loc="best",
-            fontsize=9,
-            ncol=min(2, len(active_branch_indices) + 1),
-        )
+    ax.legend(loc="best", fontsize=9, ncol=2)
 
-        Visualizer._build_figure_header(analyzer, fig)
-
-        plt.tight_layout()
+    _apply_figure_header(fig, meta)
+    fig.tight_layout()
 
     if save_path:
         fig.savefig(save_path, dpi=150, bbox_inches="tight")
         print(f"Figure saved to {save_path}")
+    return fig
 
-    @staticmethod
-    def _plot_debug_view(
-        analyzer: BetaSpectrumAnalyzer,
-        total: np.ndarray,
-        components: Dict[str, np.ndarray],
+
+def _plot_debug_view_from_df(
+    df: pd.DataFrame,
+    meta: Dict[str, Any],
     save_path: Optional[str],
-    ) -> None:
-        """Debug view: vertical layout with 4 panels."""
-        n_branches = analyzer._n_branches
-        e_step_keV = analyzer.config.e_step_MeV * 1000
+) -> Figure:
+    """Debug view (4 panels) built entirely from CSV columns."""
+    energies = df["energy_MeV"].to_numpy()
+    active = _active_branches(meta)
 
     fig, axes = plt.subplots(4, 1, figsize=(12, 16), sharex=True)
-        fig.suptitle(
-            Visualizer._build_figure_title(analyzer, fig),
-            fontsize=14,
-            fontweight="bold",
-        )
-
-        branch_colors = [
-            "#1f77b4",
-            "#ff7f0e",
-            "#2ca02c",
-            "#d62728",
-            "#9467bd",
-            "#8c564b",
-            "#e377c2",
-            "#7f7f7f",
-            "#bcbd22",
-            "#17becf",
-        ]
-
-        # Panel 1: Total spectrum + per-branch spectra
-        branch_contribs = []
-        universal_names = [
-            "Fermi",
-            "Screening",
-            "FiniteSizeL0",
-            "ChargeDistributionU",
-            "Exchange",
-        ]
-        for i in range(n_branches):
-            universal_product = np.ones_like(analyzer.W)
-            for name in universal_names:
-                key = analyzer._get_component_key(i, name)
-                if key in components:
-                    universal_product *= components[key]
-                elif name in components:
-                    universal_product *= components[name]
-            bs = analyzer._get_branch_spectrum(i)
-            contrib = universal_product * analyzer._unified_branches[i].intensity * bs
-            branch_contribs.append(contrib)
-
-        total_raw = sum(branch_contribs)
+    fig.suptitle(_spectrum_title(meta), fontsize=14, fontweight="bold")
 
+    # Panel 1: total spectrum + per-branch spectra
     ax = axes[0]
+    if "total_spectrum" in df.columns:
         ax.plot(
-            analyzer.energies_MeV,
-            total_raw,
+            energies,
+            df["total_spectrum"].to_numpy(),
             "k-",
             lw=2,
             label="Total spectrum",
             zorder=10,
         )
-
-        # Only plot active branches (not excluded by neutrino mass or intensity cutoff)
-        active_branch_indices = []
-        for i in range(n_branches):
-            branch_endpoint_keV = analyzer._unified_branches[i].endpoint_MeV * 1000
-            neutrino_mass_keV = analyzer.config.neutrino_mass_MeV * 1000
-            is_excluded_by_mass = branch_endpoint_keV <= neutrino_mass_keV
-            is_excluded_by_intensity = (
-                analyzer._unified_branches[i].intensity
-                < analyzer.config.intensity_cutoff
-            )
-            if not is_excluded_by_mass and not is_excluded_by_intensity:
-                active_branch_indices.append(i)
-
-        for i in active_branch_indices:
-            contrib = branch_contribs[i]
-            color = branch_colors[i % len(branch_colors)]
+    for color_i, b in enumerate(active):
+        col = _branch_column(df, int(b["index"]))
+        if col is None:
+            continue
         ax.plot(
-                analyzer.energies_MeV,
-                contrib,
-                color=color,
+            energies,
+            df[col].to_numpy(),
+            color=_BRANCH_COLORS[color_i % len(_BRANCH_COLORS)],
             lw=1.5,
             alpha=0.7,
-                label=Visualizer._format_branch_label(analyzer, i),
+            label=f"Br. {b['index']} ({b['intensity'] * 100:.1f}%)",
         )
     ax.set_yscale("log")
-        ax.set_ylabel(f"Normalized Counts / {e_step_keV:.1f} keV", fontsize=10)
+    ax.set_ylabel(r"Normalized spectrum [MeV$^{-1}$]", fontsize=10)
     ax.set_title("Spectrum", fontsize=11)
     ax.grid(True, alpha=0.3)
-        ax.legend(fontsize=8, loc="best", ncol=min(2, len(active_branch_indices) + 1))
+    ax.legend(fontsize=8, loc="best", ncol=2)
 
     # Panel 2: Fermi function + all phase spaces
     ax = axes[1]
-        if "Fermi" in components:
-            fermi_vals = components["Fermi"]
-            fermi_scaled = fermi_vals / np.max(fermi_vals)
+    if "Fermi" in df.columns:
+        fermi = df["Fermi"].to_numpy()
+        fermi_scaled = fermi / np.max(fermi)
+        ax.plot(energies, fermi_scaled, "r-", lw=1.5, label="Fermi", alpha=0.8)
+    for color_i, b in enumerate(active):
+        col = f"branch_{b['index']}_PhaseSpace"
+        if col not in df.columns:
+            # Single-branch CSVs carry the core factors unprefixed.
+            col = "PhaseSpace"
+        if col not in df.columns:
+            continue
+        ps_vals = df[col].to_numpy()
+        ps_max = np.max(ps_vals) if np.any(ps_vals > 0) else 1.0
         ax.plot(
-                analyzer.energies_MeV,
-                fermi_scaled,
-                "r-",
-                lw=1.5,
-                label="Fermi",
-                alpha=0.8,
-            )
-        for i in active_branch_indices:
-            color = branch_colors[i % len(branch_colors)]
-            branch = analyzer._unified_branches[i]
-            key = analyzer._get_component_key(i, "PhaseSpace")
-            if key in components:
-                ps_vals = components[key]
-                W0_branch = T_to_W(branch.endpoint_MeV)
-                mask = analyzer.W <= W0_branch
-                ps_masked = np.where(mask, ps_vals, 0.0)
-                ps_max = np.max(ps_masked)
-                if ps_max > 0:
-                    ps_scaled = ps_masked / ps_max
-                else:
-                    ps_scaled = ps_masked
-                ax.plot(
-                    analyzer.energies_MeV,
-                    ps_scaled,
-                    color=color,
+            energies,
+            ps_vals / ps_max if ps_max > 0 else ps_vals,
+            color=_BRANCH_COLORS[color_i % len(_BRANCH_COLORS)],
             lw=1.5,
             alpha=0.7,
-                    label=f"PS Br. {i + 1}",
+            label=f"PS Br. {b['index']}",
         )
     ax.set_ylabel("Normalized Factor", fontsize=10)
     ax.set_title("Fermi & Phase Space", fontsize=11)
     ax.grid(True, alpha=0.3)
     ax.legend(fontsize=8, loc="best")
 
-        # Panel 3: Universal W0-independent corrections
+    # Panel 3: universal higher-order corrections
     ax = axes[2]
-        universal_components = [
-            "Screening",
-            "FiniteSizeL0",
-            "ChargeDistributionU",
-            "Exchange",
-        ]
-        colors = ["green", "orange", "purple", "brown"]
-        for name, color in zip(universal_components, colors):
-            if name in components:
-                values = components[name]
+    correction_colors = ["green", "orange", "purple", "brown"]
+    for corr_i, name in enumerate(
+        ["Screening", "FiniteSizeL0", "ChargeDistributionU", "Exchange"]
+    ):
+        if name in df.columns:
+            values = df[name].to_numpy()
             values_scaled = values / np.max(values)
             ax.plot(
-                    analyzer.energies_MeV,
+                energies,
                 values_scaled,
-                    color=color,
+                color=correction_colors[corr_i],
                 lw=1.5,
                 label=name,
                 alpha=0.8,
@@ -414,31 +371,98 @@ class Visualizer:
     ax.grid(True, alpha=0.3)
     ax.legend(fontsize=8, loc="best")
 
-        # Panel 4: W0-dependent corrections (Radiative for each branch)
+    # Panel 4: W0-dependent corrections (Radiative per branch)
     ax = axes[3]
-        for i in active_branch_indices:
-            color = branch_colors[i % len(branch_colors)]
-            key = analyzer._get_component_key(i, "Radiative")
-            if key in components:
-                rad_vals = components[key]
+    for color_i, b in enumerate(active):
+        col = f"branch_{b['index']}_Radiative"
+        if col not in df.columns:
+            # Single-branch CSVs carry the core factors unprefixed.
+            col = "Radiative"
+        if col not in df.columns:
+            continue
         ax.plot(
-                    analyzer.energies_MeV,
-                    rad_vals,
-                    color=color,
+            energies,
+            df[col].to_numpy(),
+            color=_BRANCH_COLORS[color_i % len(_BRANCH_COLORS)],
             lw=1.5,
             alpha=0.7,
-                    label=f"Br. {i + 1}",
+            label=f"Br. {b['index']}",
         )
     ax.set_ylabel("Radiative Correction", fontsize=10)
     ax.set_xlabel(r"Electron kinetic energy $E$ [MeV]", fontsize=10)
     ax.set_title("Radiative Corrections", fontsize=11)
     ax.grid(True, alpha=0.3)
-        ax.legend(fontsize=8, loc="best", ncol=min(2, len(active_branch_indices) + 1))
+    ax.legend(fontsize=8, loc="best", ncol=2)
 
-        Visualizer._build_figure_header(analyzer, fig)
-
-        plt.tight_layout()
+    _apply_figure_header(fig, meta)
+    fig.tight_layout()
 
     if save_path:
         fig.savefig(save_path, dpi=150, bbox_inches="tight")
         print(f"Figure saved to {save_path}")
+    return fig
+
+
+class Visualizer:
+    """Visualization built from the exported CSV data (single source of truth)."""
+
+    @staticmethod
+    def plot_from_csv(
+        csv_path: str,
+        save_path: Optional[str] = None,
+        show_components: bool = True,
+    ) -> Figure:
+        """Create a visualization from an exported CSV spectrum file.
+
+        Parameters
+        ----------
+        csv_path : str
+            Path to a CSV exported by CSVEmitter.export_to_csv. The plots
+            are built exclusively from the file contents (data columns and
+            metadata header) — nothing is recomputed from the model.
+        save_path : str, optional
+            Path to save the figure.
+        show_components : bool
+            If True, produce the full 4-panel debug view; requires the
+            'full' CSV verbosity. If the file lacks component columns, a
+            warning is printed and the simple spectrum view is produced.
+        """
+        df, meta = parse_spectrum_csv(csv_path)
+        if show_components and "Fermi" in df.columns:
+            return _plot_debug_view_from_df(df, meta, save_path)
+        if show_components:
+            print(
+                "Warning: requested component panels but CSV lacks component "
+                "columns (export with --csv-verbosity full); falling back to "
+                "simple spectrum view."
+            )
+        return _plot_spectra_view_from_df(df, meta, save_path)
+
+    @staticmethod
+    def plot_analysis(
+        analyzer: BetaSpectrumAnalyzer,
+        save_path: Optional[str] = None,
+        show_components: bool = True,
+    ) -> Figure:
+        """Create a visualization by exporting CSV and plotting from it.
+
+        The export goes to a temporary file using the 'full' verbosity so
+        that the plots exercise the same data pipeline (header + columns)
+        a researcher receives.
+        """
+        tmp_dir = tempfile.mkdtemp(prefix="bs_plot_")
+        try:
+            csv_path = f"{tmp_dir}/spectrum.csv"
+            CSVEmitter.export_to_csv(
+                analyzer,
+                csv_path,
+                source_type="plot",
+                csv_verbosity="full",
+            )
+            return Visualizer.plot_from_csv(
+                csv_path,
+                save_path=save_path,
+                show_components=show_components,
+            )
+        finally:
+            shutil.rmtree(tmp_dir, ignore_errors=True)

beta_spectrum/spectrum.py

@@ -57,9 +57,6 @@ class SpectrumConfig:
     """
     Configuration for beta spectrum calculation.
 
-    Supports declarative detector response specification:
-    to enable automatic detector smearing via convolve_with_detector().
-
     Supports multi-branch decay via the `branches` field. When branches are
     provided, the total spectrum is the intensity-weighted sum of individual
     branch spectra. Each branch has its own endpoint energy and transition type.
@@ -73,9 +70,9 @@ class SpectrumConfig:
     e_step_MeV: float = 0.001
     neutrino_mass_MeV: float = 0.0
 
-    # Toggle components and corrections
-    use_phase_space: bool = True
-    use_fermi: bool = True
+    # Higher-order corrections. The core spectrum components (phase space
+    # and the Fermi function, both 0-order) are unconditional: there are
+    # deliberately no toggles for them.
     use_screening: bool = True
     use_finite_size: bool = True
     use_charge_dist: bool = True
@@ -90,9 +87,25 @@ class SpectrumConfig:
 
     # CSV output options
     csv_separator: str = ","  # Default separator for CSV output
-    normalize_to_one: bool = (
-        False  # If True, normalize each branch to 1, ignoring intensity
+    separate_branches: bool = (
+        False  # If True: emit separate branch spectra, each normalized to
     )
+    # integrate to 1 (intensity ignored); no total spectrum column
+
+
+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:
@@ -252,23 +265,27 @@ class BetaSpectrum:
     def _from_config_single(
         cls, config: SpectrumConfig, logger: Optional[logging.Logger] = None
     ) -> BetaSpectrum:
-        """Create single-branch spectrum calculator."""
-        components: List[SpectrumComponent] = []
-        W0 = float(T_to_W(config.endpoint_MeV))
+        """Create single-branch spectrum calculator.
 
-        if config.use_phase_space:
+        The core spectrum components (phase space x Fermi function, both
+        0-order) are always built; only higher-order corrections are
+        conditional.
+        """
+        components: List[SpectrumComponent] = []
+        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
         m_nu_natural = config.neutrino_mass_MeV / ME_MEV
         components.append(
             PhaseSpace(
                 W0=W0,
-                    transition_type=config.transition_type,
+                transition_type=eff_transition,
                 m_nu=m_nu_natural,
                 logger=logger,
             )
         )
-
-        if config.use_fermi:
         components.append(
             FermiFunction(Z=config.Z_daughter, A=config.A_number, logger=logger)
         )
@@ -322,13 +339,13 @@ class BetaSpectrum:
         Note: Radiative correction is branch-dependent (depends on W0), so it goes
         into each branch's calculator, not the universal components.
         """
+        if config.branches is None:
+            raise ValueError("branches must be provided for multi-branch mode")
         if logger:
             logger.info("Multi-branch mode: %d branches", len(config.branches))
 
         # Create per-branch calculators
         branch_spectra: List[BetaSpectrum] = []
-        if config.branches is None:
-            raise ValueError("branches must be provided for multi-branch mode")
         for i, branch in enumerate(config.branches):
             if logger:
                 logger.info(
@@ -339,37 +356,42 @@ class BetaSpectrum:
                     branch.intensity,
                 )
 
-            # Create a per-branch config with the branch's endpoint and transition type
-            branch_config = SpectrumConfig(
-                Z_parent=config.Z_parent,
-                Z_daughter=config.Z_daughter,
-                A_number=config.A_number,
-                endpoint_MeV=branch.endpoint_MeV,
+            # Per-branch calculator: the per-branch core phase space
+            # (unconditional) plus optional per-branch radiative
+            # corrections. The universal Fermi function and universal
+            # higher-order corrections are applied outside, so they must
+            # NOT be built here (a Fermi function per branch would
+            # double-apply it).
+            W0_branch = float(T_to_W(branch.endpoint_MeV))
+            m_nu_natural = config.neutrino_mass_MeV / ME_MEV
+            branch_components: List[SpectrumComponent] = [
+                PhaseSpace(
+                    W0=W0_branch,
                     transition_type=branch.transition_type,
-                e_step_MeV=config.e_step_MeV,
-                neutrino_mass_MeV=config.neutrino_mass_MeV,
-                # Branch-dependent components
-                use_phase_space=config.use_phase_space,
-                use_radiative=config.use_radiative,
-                # Universal components disabled (handled separately)
-                use_fermi=False,
-                use_screening=False,
-                use_finite_size=False,
-                use_charge_dist=False,
-                use_exchange=False,
+                    m_nu=m_nu_natural,
+                    logger=logger,
+                )
+            ]
+            if config.use_radiative:
+                branch_components.append(
+                    RadiativeCorrection(
+                        W0=W0_branch,
+                        Z=config.Z_daughter,
+                        A=config.A_number,
+                        use_endpoint_resummation=True,
+                        logger=logger,
+                    )
                 )
 
-            branch_spectrum = cls._from_config_single(branch_config, logger)
-            branch_spectra.append(branch_spectrum)
+            branch_spectra.append(cls(branch_components, logger=logger))
 
         # Create universal components (shared across all branches)
         # Note: Radiative is NOT included here because it depends on W0
-        universal_components: List[SpectrumComponent] = []
-
-        if config.use_fermi:
-            universal_components.append(
+        # The Fermi function is a core spectrum component (0-order Coulomb
+        # correction): always included, not a higher-order correction.
+        universal_components: List[SpectrumComponent] = [
             FermiFunction(Z=config.Z_daughter, A=config.A_number, logger=logger)
-            )
+        ]
 
         if config.use_finite_size:
             universal_components.append(
@@ -411,6 +433,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
 

data/custom_input_example.json

@@ -5,14 +5,5 @@
     "A_number": 99,
     "endpoint_MeV": 0.294,
     "transition_type": "A",
-    "e_step_MeV": 0.001,
-    "use_detector_response": false,
-    "detector_model": "gaussian",
-    "detector_sigma_a_keV": 1.0,
-    "detector_sigma_b": 0.0,
-    "detector_tail_fraction": 0.0,
-    "detector_tau_keV": 5.0,
-    "detector_fano_factor": 0.12,
-    "detector_n_channels": 4096,
-    "detector_channel_energy_range": [0.0, 0.35]
+    "e_step_MeV": 0.001
 }

fitter/model.py

@@ -11,13 +11,12 @@ Usage:
 
 from __future__ import annotations
 
-from typing import TYPE_CHECKING, Optional
+from typing import TYPE_CHECKING, Any, Optional
 
 import numpy as np
 
 if TYPE_CHECKING:
     from beta_spectrum.spectrum import BetaSpectrum
-    from beta_spectrum.components.detector_response import DetectorResponse
     from exp_data.spectrum import ExpSpectrum
 
 
@@ -39,8 +38,11 @@ class SpectrumModel:
     ----------
     spectrum : BetaSpectrum
         Theoretical spectrum calculator.
-    detector_response : DetectorResponse, optional
-        Detector response function. If None, no convolution is applied.
+    detector_response : optional
+        Detector response object exposing ``channel_energies`` and
+        ``convolve`` (duck-typed: the DetectorResponse module was removed
+        from beta_spectrum and is pending a restore-or-delete decision).
+        If None, no convolution is applied.
     channel_energies : np.ndarray, optional
         Energy values for detector channels (keV). Used to set up
         the response if detector_response is None.
@@ -49,7 +51,7 @@ class SpectrumModel:
     def __init__(
         self,
         spectrum: "BetaSpectrum",
-        detector_response: Optional["DetectorResponse"] = None,
+        detector_response: Optional[Any] = None,
         channel_energies: Optional[np.ndarray] = None,
     ):
         self.spectrum = spectrum
@@ -79,11 +81,10 @@ class SpectrumModel:
         if energies_keV is not None:
             W = self._keV_to_W(energies_keV)
             return self.spectrum(W)
-        return self.spectrum(
-            self.spectrum._W_grid
-            if hasattr(self.spectrum, "_W_grid") and self.spectrum._W_grid is not None
-            else np.array([])
-        )
+        w_grid = getattr(self.spectrum, "_W_grid", None)
+        if w_grid is None:
+            return np.array([], dtype=np.float64)
+        return self.spectrum(w_grid)
 
     def _keV_to_W(self, energies_keV: np.ndarray) -> np.ndarray:
         """Convert kinetic energy in keV to total energy in m_e units."""

notebooks/01-beta-spectrum-calculator.ipynb

@@ -94,8 +94,6 @@
     "    endpoint_MeV=0.294,\n",
     "    transition_type=\"F2\",\n",
     "    e_step_MeV=0.001,\n",
-    "    use_phase_space=True,\n",
-    "    use_fermi=True,\n",
     "    use_screening=True,\n",
     "    use_finite_size=True,\n",
     "    use_charge_dist=True,\n",
@@ -410,13 +408,16 @@
     }
    ],
    "source": [
+    "multi_analyzer = BetaSpectrumAnalyzer(multi_spectrum, multi_config)\n",
+    "\n",
     "fig, ax = plt.subplots(figsize=(12, 6))\n",
     "ax.plot(T_keV, total, \"k-\", linewidth=2.5, label=\"Total spectrum\")\n",
     "\n",
-    "branch_spectra = multi_spectrum.get_branch_normalized_spectra(W)\n",
-    "for i, bs in enumerate(branch_spectra):\n",
-    "    ax.plot(T_keV, bs, alpha=0.6, linewidth=1.5,\n",
-    "            label=\"Branch {} (E0={:.0f} keV)\".format(i+1, branches[i].endpoint_MeV*1000))\n",
+    "# Physical branch spectra: core (phase_space x fermi) x all enabled corrections\n",
+    "for i, b in enumerate(branches):\n",
+    "    physical = multi_analyzer.physical_branch_spectrum(i)\n",
+    "    ax.plot(T_keV, physical, alpha=0.6, linewidth=1.5,\n",
+    "            label=\"Branch {} (E0={:.0f} keV)\".format(i+1, b.endpoint_MeV*1000))\n",
     "\n",
     "for i, b in enumerate(branches):\n",
     "    br_end_keV = W_to_T(T_to_W(b.endpoint_MeV))\n",
@@ -430,7 +431,7 @@
     "ax.set_xlim(0, W0_max_keV)\n",
     "ax.grid(True, alpha=0.3)\n",
     "plt.tight_layout()\n",
-    "plt.show()"
+    "plt.show()\n"
    ]
   },
   {
@@ -541,8 +542,8 @@
     }
    ],
    "source": [
-    "analyzer = BetaSpectrumAnalyzer(multi_spectrum, multi_config)\n",
-    "analyzer.plot_analysis(show_components=True)"
+    "analyzer = multi_analyzer\n",
+    "analyzer.plot_analysis(show_components=True)\n"
    ]
   },
   {

notebooks/exp_calibration.ipynb

@@ -81,11 +81,7 @@
     "\n",
     "# Detect header size from file\n",
     "def detect_header_size(filepath):\n",
-    "    \"\"\"Detect MCA header size by checking for binary pattern.\"\"\"\n",
-    "    with open(filepath, 'rb') as f:\n",
-    "        header = f.read(128)\n",
-    "    # Check if first 32 bytes look like a header (small ints)\n",
-    "    # and bytes 32-128 also look like header\n",
+    "    \"\"\"Detect MCA header size from the file size.\"\"\"\n",
     "    # Simple heuristic: if file is exactly 32 + 4096*4 = 16412, header=32\n",
     "    # If file is 128 + 4096*4 = 16512, header=128\n",
     "    import os\n",
@@ -97,10 +93,6 @@
     "    elif size == expected_128:\n",
     "        return 128\n",
     "    else:\n",
-    "        # Try to detect by checking if byte 32 looks like an int\n",
-    "        with open(filepath, 'rb') as f:\n",
-    "            f.seek(32)\n",
-    "            test = f.read(4)\n",
     "        return 32  # default\n",
     "\n",
     "# --- Raw run spectrum ---\n",

pyproject.toml

@@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta"
 
 [project]
 name = "beta-spectrum"
-version = "0.4.0"
+version = "0.4.1"
 description = "High-precision beta spectrum calculation for nuclear physics"
 authors = [
     { name = "Evgeniy Unzhakov" }

scripts/analyze_raw.py

@@ -17,6 +17,7 @@ from __future__ import annotations
 
 import sys
 from pathlib import Path
+from typing import cast
 
 import numpy as np
 from scipy.stats import linregress
@@ -101,7 +102,10 @@ def estimate_endpoint_from_counts(
 
     # Linear fit: K = a + b * E, endpoint = -a/b
     try:
-        slope, intercept, r_val, _, _ = linregress(E_sel, K)[:5]
+        lin = cast("tuple[float, float, float, float, float]", linregress(E_sel, K))
+        slope = lin[0]
+        intercept = lin[1]
+        r_val = lin[2]
         if slope < 0 and abs(r_val) > 0.3:
             endpoint = -intercept / slope
             # Clamp to reasonable range (Q ± 100 keV)
@@ -129,7 +133,7 @@ def compute_regional_rates(counts: np.ndarray, energies_keV: np.ndarray) -> dict
     }
 
 
-def main():
+def main() -> None:
     """Main analysis pipeline."""
     data_dir = Path(__file__).parent.parent / "data" / "raw"
     output_dir = Path(__file__).parent.parent / "output"
@@ -168,8 +172,6 @@ def main():
     print("\nSubtracting consecutive runs...")
     individual_runs = []
     for i in range(1, len(cumulative_specs)):
-        parent_a = cumulative_specs[i].run_id
-        parent_b = cumulative_specs[i - 1].run_id
         diff = subtract_spectra(
             cumulative_a=cumulative_specs[i],
             cumulative_b=cumulative_specs[i - 1],
@@ -215,33 +217,50 @@ def main():
     print(f"\n{'=' * 60}")
     print(f"Analysis of {len(run_ids)} individual Tc-99 runs")
     print(f"{'=' * 60}")
-    print(f"Run range: {run_ids.min()} - {run_ids.max()}")
-    print(f"Total counts: {total_rates.min():.0f} - {total_rates.max():.0f}")
-    print(f"Low-energy fraction: {low_fractions.min():.3f} - {low_fractions.max():.3f}")
-    print(f"Endpoint estimate: {endpoints.min():.1f} - {endpoints.max():.1f} keV")
+    print(f"Run range: {min(run_ids)} - {max(run_ids)}")
+    print(f"Total counts: {min(total_rates):.0f} - {max(total_rates):.0f}")
+    print(f"Low-energy fraction: {min(low_fractions):.3f} - {max(low_fractions):.3f}")
+    print(f"Endpoint estimate: {min(endpoints):.1f} - {max(endpoints):.1f} keV")
 
     # Endpoint drift
     if len(run_ids) > 10:
-        slope, intercept, r_val, _, _ = linregress(run_ids, endpoints)
-        print(f"\nEndpoint drift: {slope:.3f} keV/run (R²={r_val**2:.4f})")
-        print(f"  Start (run 1): {intercept + slope * run_ids[0]:.1f} keV")
-        print(f"  End (run {run_ids[-1]}): {intercept + slope * run_ids[-1]:.1f} keV")
+        lin_e = cast(
+            "tuple[float, float, float, float, float]",
+            linregress(run_ids, endpoints),
+        )
+        slope_e = lin_e[0]
+        intercept_e = lin_e[1]
+        r_e = lin_e[2]
+        print(f"\nEndpoint drift: {slope_e:.3f} keV/run (R²={r_e**2:.4f})")
+        print(f"  Start (run 1): {intercept_e + slope_e * run_ids[0]:.1f} keV")
+        print(
+            f"  End (run {run_ids[-1]}): {intercept_e + slope_e * run_ids[-1]:.1f} keV"
+        )
 
     # Rate drift
-    slope, intercept, r_val, _, _ = linregress(run_ids, total_rates)
-    print(f"\nTotal count rate drift: {slope:.0f} counts/run (R²={r_val**2:.4f})")
+    lin_r = cast(
+        "tuple[float, float, float, float, float]", linregress(run_ids, total_rates)
+    )
+    slope_r = lin_r[0]
+    r_r = lin_r[2]
+    print(f"\nTotal count rate drift: {slope_r:.0f} counts/run (R²={r_r**2:.4f})")
 
     # Save results
     results_file = output_dir / "rate_analysis.csv"
+    try:
         with open(results_file, "w") as f:
             f.write(
-            "run,total_counts,low_energy_counts,high_energy_counts,low_fraction,endpoint_keV\n"
+                "run,total_counts,low_energy_counts,high_energy_counts,"
+                "low_fraction,endpoint_keV\n"
             )
             for i in range(len(run_ids)):
                 f.write(
                     f"{run_ids[i]},{total_rates[i]:.0f},{low_rates[i]:.0f},"
                     f"{high_rates[i]:.0f},{low_fractions[i]:.4f},{endpoints[i]:.1f}\n"
                 )
+    except OSError as exc:
+        print(f"Cannot write {results_file}: {exc}")
+        raise SystemExit(1) from exc
     print(f"\nResults saved to {results_file}")
 
     # Plot

scripts/docs_check.py

@@ -471,7 +471,7 @@ CHECKS = [
 ]
 
 
-def main():
+def main() -> int:
     parser = argparse.ArgumentParser(description="Check documentation consistency")
     parser.add_argument(
         "--strict", action="store_true", help="Treat warnings as errors (exit 1)"

tests/common/test_hypothesis_tests.py

@@ -10,7 +10,6 @@ seconds to complete. They are intended for CI/full test runs, not for
 fast feedback during development.
 """
 
-
 import numpy as np
 from hypothesis import given, settings, strategies as st
 

tests/common/test_property_tests.py

@@ -116,14 +116,13 @@ def _test_non_negative_values(
     ), f"{component_name}: Inf values found in output."
 
     if allow_zero:
-        assert np.all(result >= 0), (
-            f"{component_name}: negative values found. " f"Min = {np.min(result):.10f}"
-        )
+        assert np.all(
+            result >= 0
+        ), f"{component_name}: negative values found. Min = {np.min(result):.10f}"
     else:
-        assert np.all(result > 0), (
-            f"{component_name}: non-positive values found. "
-            f"Min = {np.min(result):.10f}"
-        )
+        assert np.all(
+            result > 0
+        ), f"{component_name}: non-positive values found. Min = {np.min(result):.10f}"
 
 
 # ---------------------------------------------------------------------------
@@ -140,10 +139,9 @@ def _test_output_shape_preserved(component, W, component_name="Component"):
     assert isinstance(
         result, np.ndarray
     ), f"{component_name}: output is {type(result)}, expected np.ndarray"
-    assert result.shape == W_arr.shape, (
-        f"{component_name}: shape mismatch. "
-        f"Input {W_arr.shape} → Output {result.shape}"
-    )
+    assert (
+        result.shape == W_arr.shape
+    ), f"{component_name}: shape mismatch. Input {W_arr.shape} → Output {result.shape}"
 
 
 @pytest.mark.quality
@@ -216,10 +214,9 @@ def _test_output_is_ndarray(component, W, component_name="Component"):
     W_arr = np.asarray(W, dtype=float)
     result = component(W_arr)
 
-    assert isinstance(result, np.ndarray), (
-        f"{component_name}: output type is {type(result).__name__}, "
-        f"expected np.ndarray"
-    )
+    assert isinstance(
+        result, np.ndarray
+    ), f"{component_name}: output type is {type(result).__name__}, expected np.ndarray"
 
 
 # ---------------------------------------------------------------------------

tests/fitter/test_fit_engine.py

@@ -9,7 +9,7 @@ from fitter.model import SpectrumModel
 class TestFitConfig:
     """Test FitConfig defaults."""
 
-    def test_defaults(self):
+    def test_defaults(self) -> None:
         config = FitConfig()
         assert config.method == "trf"
         assert config.max_nfev == 1000
@@ -19,11 +19,13 @@ class TestFitConfig:
 class TestSpectrumFitter:
     """Test SpectrumFitter."""
 
-    def _make_fitter(self):
+    def _make_fitter(self) -> "SpectrumFitter":
         """Create a simple fitter with synthetic data."""
 
         # Create a simple linear model
-        def linear_model(energies_keV, norm, bg):
+        def linear_model(
+            energies_keV: np.ndarray, norm: float, bg: float
+        ) -> np.ndarray:
             return norm * (1.0 - energies_keV / 100.0) + bg
 
         # Synthetic data
@@ -33,7 +35,9 @@ class TestSpectrumFitter:
 
         # Mock model
         class MockModel:
-            def evaluate(self, energies_keV, norm, bg):
+            def evaluate(
+                self, energies_keV: np.ndarray, norm: float, bg: float
+            ) -> np.ndarray:
                 return linear_model(energies_keV, norm, bg)
 
         model = MockModel()
@@ -42,7 +46,7 @@ class TestSpectrumFitter:
         )
         return fitter
 
-    def test_fit_basic(self):
+    def test_fit_basic(self) -> None:
         fitter = self._make_fitter()
         result = fitter.fit(x0=[10.0, 2.0], param_names=["norm", "bg"])
 
@@ -53,7 +57,7 @@ class TestSpectrumFitter:
         assert result.chi2_per_dof > 0
         assert len(result.parameters) == 2
 
-    def test_fit_with_bounds(self):
+    def test_fit_with_bounds(self) -> None:
         fitter = self._make_fitter()
         result = fitter.fit(
             x0=[10.0, 2.0],
@@ -61,36 +65,36 @@ class TestSpectrumFitter:
         )
         assert result.success
 
-    def test_fit_returns_residuals(self):
+    def test_fit_returns_residuals(self) -> None:
         fitter = self._make_fitter()
         result = fitter.fit(x0=[10.0, 2.0])
 
         assert len(result.residuals) == 200
         assert len(result.model_values) == 200
 
-    def test_parameters_property(self):
+    def test_parameters_property(self) -> None:
         fitter = self._make_fitter()
         fitter.fit(x0=[10.0, 2.0])
         params = fitter.parameters
         assert len(params) == 2
 
-    def test_covariance_property(self):
+    def test_covariance_property(self) -> None:
         fitter = self._make_fitter()
         fitter.fit(x0=[10.0, 2.0])
         cov = fitter.covariance
         assert cov.shape == (2, 2)
 
-    def test_chi2_per_dof(self):
+    def test_chi2_per_dof(self) -> None:
         fitter = self._make_fitter()
         result = fitter.fit(x0=[10.0, 2.0])
         assert result.chi2_per_dof == result.chi2 / result.n_free
 
-    def test_p_value(self):
+    def test_p_value(self) -> None:
         fitter = self._make_fitter()
         result = fitter.fit(x0=[10.0, 2.0])
         assert 0.0 <= result.p_value <= 1.0
 
-    def test_parameters_with_errors(self):
+    def test_parameters_with_errors(self) -> None:
         fitter = self._make_fitter()
         result = fitter.fit(x0=[10.0, 2.0])
         params = result.parameters_with_errors
@@ -100,7 +104,7 @@ class TestSpectrumFitter:
             assert isinstance(err, float)
             assert err >= 0
 
-    def test_correlation_matrix(self):
+    def test_correlation_matrix(self) -> None:
         fitter = self._make_fitter()
         result = fitter.fit(x0=[10.0, 2.0])
         corr = result.correlation_matrix
@@ -108,7 +112,7 @@ class TestSpectrumFitter:
         # Diagonal should be 1
         np.testing.assert_allclose(np.diag(corr), 1.0)
 
-    def test_summary(self):
+    def test_summary(self) -> None:
         fitter = self._make_fitter()
         result = fitter.fit(x0=[10.0, 2.0], param_names=["norm", "bg"])
         summary = result.summary(param_names=["norm", "bg"])
@@ -116,7 +120,7 @@ class TestSpectrumFitter:
         assert "bg" in summary
         assert "Fit Results" in summary
 
-    def test_profile_likelihood(self):
+    def test_profile_likelihood(self) -> None:
         fitter = self._make_fitter()
         result = fitter.fit(x0=[10.0, 2.0])
         param_vals, chi2_vals, norm_vals = result.profile_likelihood(0, n_points=10)
@@ -125,7 +129,7 @@ class TestSpectrumFitter:
         assert len(norm_vals) == 10
         assert np.min(norm_vals) >= 0  # Normalized χ² should be >= 0
 
-    def test_confidence_interval(self):
+    def test_confidence_interval(self) -> None:
         fitter = self._make_fitter()
         result = fitter.fit(x0=[10.0, 2.0])
         ci = result.confidence_interval(0, confidence=0.95)
@@ -136,7 +140,7 @@ class TestSpectrumFitter:
 class TestSpectrumFitterWithRealModel:
     """Test SpectrumFitter with a real SpectrumModel."""
 
-    def test_fit_with_real_model(self):
+    def test_fit_with_real_model(self) -> None:
         from beta_spectrum import BetaSpectrum, SpectrumConfig
 
         config = SpectrumConfig(
@@ -144,16 +148,12 @@ class TestSpectrumFitterWithRealModel:
             Z_daughter=44,
             A_number=99,
             endpoint_MeV=0.294,
-            detector_model="gaussian",
-            detector_sigma_a_keV=0.5,
-            detector_n_channels=200,
         )
         spectrum = BetaSpectrum.from_config(config)
-        detector = BetaSpectrum.create_detector_from_config(config)
-        model = SpectrumModel(spectrum, detector)
+        model = SpectrumModel(spectrum)
 
-        # Generate synthetic data from model
-        energies_keV = detector.channel_energies
+        # Generate synthetic data from the unconvolved model on a keV grid
+        energies_keV = np.linspace(1.0, 293.0, 200)
         true_norm = 5.0
         true_bg = 0.1
         counts = model.evaluate(energies_keV, true_norm, true_bg)

tests/fitter/test_model.py

@@ -10,8 +10,8 @@ from fitter.model import SpectrumModel
 class TestSpectrumModel:
     """Test SpectrumModel."""
 
-    def _make_model(self, with_detector: bool = False):
-        """Create a test SpectrumModel."""
+    def _make_model(self) -> tuple[SpectrumModel, SpectrumConfig]:
+        """Create a test SpectrumModel (no detector response)."""
         config = SpectrumConfig(
             Z_parent=43,
             Z_daughter=44,
@@ -20,24 +20,10 @@ class TestSpectrumModel:
         )
         spectrum = BetaSpectrum.from_config(config)
 
-        detector = None
-        if with_detector:
-            detector = BetaSpectrum.create_detector_from_config(
-                SpectrumConfig(
-                    Z_parent=43,
-                    Z_daughter=44,
-                    A_number=99,
-                    endpoint_MeV=0.294,
-                    detector_model="gaussian",
-                    detector_sigma_a_keV=1.0,
-                    detector_n_channels=100,
-                )
-            )
-
-        model = SpectrumModel(spectrum, detector)
+        model = SpectrumModel(spectrum)
         return model, config
 
-    def test_evaluate_unconvolved(self):
+    def test_evaluate_unconvolved(self) -> None:
         model, config = self._make_model()
         energies_keV = np.linspace(0.01, 0.29, 50)
         values = model.evaluate_unconvolved(energies_keV)
@@ -45,44 +31,34 @@ class TestSpectrumModel:
         assert np.all(np.isfinite(values))
         assert np.all(values >= 0)
 
-    def test_evaluate_no_detector(self):
+    def test_evaluate_no_detector(self) -> None:
         model, config = self._make_model()
         energies_keV = np.linspace(0.01, 0.29, 50)
         values = model.evaluate(energies_keV, normalization=1.0, background=0.0)
         assert len(values) == 50
         assert np.all(np.isfinite(values))
 
-    def test_evaluate_with_normalization(self):
+    def test_evaluate_with_normalization(self) -> None:
         model, config = self._make_model()
         energies_keV = np.linspace(0.01, 0.29, 50)
         values1 = model.evaluate(energies_keV, normalization=1.0, background=0.0)
         values2 = model.evaluate(energies_keV, normalization=2.0, background=0.0)
         np.testing.assert_allclose(values2, values1 * 2.0)
 
-    def test_evaluate_with_background(self):
+    def test_evaluate_with_background(self) -> None:
         model, config = self._make_model()
         energies_keV = np.linspace(0.01, 0.29, 50)
         values1 = model.evaluate(energies_keV, normalization=1.0, background=0.0)
         values2 = model.evaluate(energies_keV, normalization=1.0, background=5.0)
         np.testing.assert_allclose(values2, values1 + 5.0)
 
-    def test_convolve_requires_detector(self):
+    def test_convolve_requires_detector(self) -> None:
         model, _ = self._make_model()
         energies_keV = np.linspace(0.01, 0.29, 50)
         with pytest.raises(RuntimeError):
             model.convolve(energies_keV)
 
-    def test_convolve_with_detector(self):
-        model, _ = self._make_model(with_detector=True)
-        energies_keV = model.detector_response.channel_energies
-        values = model.convolve(energies_keV)
-        assert len(values) == model.detector_response.n_channels
-        assert np.all(np.isfinite(values))
-        # Convolution can produce small negative values due to numerical effects
-        # — check that most values are non-negative
-        assert np.sum(values >= 0) > len(values) * 0.9
-
-    def test_default_params(self):
+    def test_default_params(self) -> None:
         from exp_data.spectrum import ExpSpectrum
 
         energies = np.linspace(0, 100, 200)
@@ -93,7 +69,7 @@ class TestSpectrumModel:
         assert norm > 0
         assert bg >= 0
 
-    def test_default_bounds(self):
+    def test_default_bounds(self) -> None:
         lower, upper = SpectrumModel.default_bounds(100)
         assert len(lower) == 2
         assert len(upper) == 2
@@ -104,8 +80,11 @@ class TestSpectrumModel:
 class TestSpectrumModelKeVToW:
     """Test energy unit conversion in SpectrumModel."""
 
-    def test_keV_to_W(self):
-        model = SpectrumModel(None)  # type: ignore
+    def test_keV_to_W(self) -> None:
+        config = SpectrumConfig(
+            Z_parent=43, Z_daughter=44, A_number=99, endpoint_MeV=0.294
+        )
+        model = SpectrumModel(BetaSpectrum.from_config(config))
         # 0 keV → W = 1.0
         W = model._keV_to_W(np.array([0.0]))
         assert np.isclose(W[0], 1.0)

tests/integration/test_cli_flows.py

@@ -0,0 +1,208 @@
+"""CLI flow tests: command-line options must reach the config in every
+data-source flow (JSON file and paceENSDF), and degenerate exports must
+exit non-zero without writing a file (regression D2/D4).
+
+The CLI is invoked in-process via beta_spectrum.cli.main(argv) so the
+paceENSDF dataset cache (loaded once per process) is shared with the
+rest of the test session.
+"""
+
+from __future__ import annotations
+
+import json
+from typing import Any, Dict, List, Optional
+
+import numpy as np
+import pandas as pd
+import pytest
+
+from beta_spectrum.cli import main
+
+pytestmark = pytest.mark.integration
+
+
+def _base_json() -> Dict[str, Any]:
+    return {
+        "parent_symbol": "Tc99",
+        "Z_parent": 43,
+        "Z_daughter": 44,
+        "A_number": 99,
+        "endpoint_MeV": 0.294,
+        "transition_type": "A",
+    }
+
+
+def write_json(tmp_path: Any, extra: Optional[Dict[str, Any]] = None) -> str:
+    data = _base_json()
+    if extra:
+        data.update(extra)
+    path = tmp_path / "input.json"
+    path.write_text(json.dumps(data))
+    return str(path)
+
+
+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 read_export(path: str) -> pd.DataFrame:
+    return pd.read_csv(path, sep=None, engine="python", comment="#")
+
+
+def n_points_of(header: List[str]) -> int:
+    for h in header:
+        if h.startswith("# n_points:"):
+            return int(h.split(":")[1])
+    raise AssertionError(f"n_points header missing: {header}")
+
+
+class TestCLIJsonFlows:
+    """Options must propagate through the JSON config flow (D2)."""
+
+    def test_json_honors_separate_branches_and_tab_separator(
+        self, tmp_path: Any
+    ) -> None:
+        src = write_json(tmp_path)
+        out = str(tmp_path / "to_one.csv")
+        main([src, "--output", out, "--separate-branches", "--csv-separator", "tab"])
+
+        with open(out) as f:
+            body = "".join(ln for ln in f if not ln.startswith("#"))
+        assert "\t" in body, "CSV must be tab-separated"
+
+        df = read_export(out)
+        assert "total_spectrum" not in df.columns
+        area = np.trapezoid(df["branch_1_spectrum"], df["energy_MeV"])
+        assert np.isclose(area, 1.0, rtol=1e-3)
+
+    def test_json_defaults_keep_total_column(self, tmp_path: Any) -> None:
+        src = write_json(tmp_path)
+        out = str(tmp_path / "default.csv")
+        main([src, "--output", out])
+
+        df = read_export(out)
+        assert "total_spectrum" in df.columns
+        area = np.trapezoid(df["total_spectrum"], df["energy_MeV"])
+        assert np.isclose(area, 1.0, rtol=1e-3)
+
+    def test_json_cli_overrides_file_physics_values(self, tmp_path: Any) -> None:
+        src = write_json(tmp_path, {"e_step_MeV": 0.001, "neutrino_mass_MeV": 0.0})
+        out = str(tmp_path / "overrides.csv")
+        main([src, "--output", out, "--e-step", "0.02", "--neutrino-mass", "0.15"])
+
+        header = header_lines_of(out)
+        assert any("neutrino_mass_MeV: 0.150000" in h for h in header)
+        expected_points = int(np.arange(0.02, 0.294, 0.02).size)
+        assert n_points_of(header) == expected_points
+
+    def test_json_file_values_used_when_flags_absent(self, tmp_path: Any) -> None:
+        src = write_json(tmp_path, {"e_step_MeV": 0.005})
+        out = str(tmp_path / "file_values.csv")
+        main([src, "--output", out])
+
+        header = header_lines_of(out)
+        expected_points = int(np.arange(0.005, 0.294, 0.005).size)
+        assert n_points_of(header) == expected_points
+
+
+class TestCLIDegenerateFlows:
+    """Degenerate exports (m_nu above every endpoint) must fail cleanly."""
+
+    def test_degenerate_intensity_mode_exits_nonzero_no_file(
+        self, tmp_path: Any
+    ) -> None:
+        src = write_json(tmp_path)
+        out = str(tmp_path / "deg.csv")
+        with pytest.raises(SystemExit) as exc:
+            main([src, "--output", out, "--neutrino-mass", "0.35"])
+        assert exc.value.code == 1
+        assert not (tmp_path / "deg.csv").exists()
+
+    def test_degenerate_to_one_mode_exits_nonzero_no_file(self, tmp_path: Any) -> None:
+        src = write_json(tmp_path)
+        out = str(tmp_path / "deg_to_one.csv")
+        with pytest.raises(SystemExit) as exc:
+            main(
+                [
+                    src,
+                    "--output",
+                    out,
+                    "--neutrino-mass",
+                    "0.35",
+                    "--separate-branches",
+                ]
+            )
+        assert exc.value.code == 1
+        assert not (tmp_path / "deg_to_one.csv").exists()
+
+
+@pytest.mark.requires_ensdf
+class TestCLIPaceENSDF:
+    """paceENSDF flow honors the same export options."""
+
+    def test_paceensdf_honors_separate_branches(self, tmp_path: Any) -> None:
+        out = str(tmp_path / "p32.csv")
+        main(["P32", "--output", out, "--separate-branches"])
+
+        df = read_export(out)
+        assert "total_spectrum" not in df.columns
+        area = np.trapezoid(df["branch_1_spectrum"], df["energy_MeV"])
+        assert np.isclose(area, 1.0, rtol=1e-3)
+
+    def test_paceensdf_degenerate_exits_nonzero_no_file(self, tmp_path: Any) -> None:
+        out = str(tmp_path / "tc99deg.csv")
+        with pytest.raises(SystemExit) as exc:
+            main(["Tc99", "--output", out, "--neutrino-mass", "0.35"])
+        assert exc.value.code == 1
+        assert not (tmp_path / "tc99deg.csv").exists()
+
+    def test_paceensdf_co60_dry_run_multi_branch(self, capsys: Any) -> None:
+        """Co60 parses into a healthy multi-branch config (D3): no crash,
+        all branch endpoints positive."""
+        main(["Co60", "--dry-run"])
+        out = capsys.readouterr().out
+        assert "Multi-branch mode: 3 branches" in out
+
+    def test_paceensdf_co60_export(self, tmp_path: Any) -> None:
+        out = str(tmp_path / "co60.csv")
+        main(["Co60", "--output", out])
+        df = read_export(out)
+        area = np.trapezoid(df["total_spectrum"], df["energy_MeV"])
+        assert np.isclose(area, 1.0, rtol=1e-3)
+        eps = header_lines_of(out)
+        assert any("E₀=1492.0" in h for h in eps)
+        assert not any("excluded" in h for h in eps)
+
+    def test_cli_builds_paceensdf_config_via_library(
+        self, monkeypatch: Any, capsys: Any
+    ) -> None:
+        """CLI must delegate paceENSDF config building to the library
+        (create_config_from_source), not duplicate it inline (D5)."""
+        import beta_spectrum.nuclear_data as nuclear_data
+
+        real_builder = nuclear_data.create_config_from_source
+        calls: List[Dict[str, Any]] = []
+
+        def spy(source: str, **kwargs: Any) -> Any:
+            calls.append({"source": source, **kwargs})
+            return real_builder(source, **kwargs)
+
+        monkeypatch.setattr(nuclear_data, "create_config_from_source", spy)
+
+        main(["Cs137", "--dry-run", "--intensity-cutoff", "0.9"])
+        capsys.readouterr()
+
+        assert len(calls) == 1, (
+            "CLI must delegate paceENSDF config building to the library "
+            "builder exactly once"
+        )
+        assert calls[0]["source"] == "paceENSDF"
+        assert calls[0]["nuclide"] == "Cs137"
+        assert calls[0]["intensity_cutoff"] == 0.9

tests/integration/test_spectrum_pipeline.py

@@ -7,7 +7,6 @@ work together correctly.
 """
 
 import numpy as np
-import pytest
 
 from beta_spectrum import BetaSpectrum, SpectrumConfig
 
@@ -61,13 +60,12 @@ class TestBetaSpectrumToggles:
     """Test that component toggling works correctly."""
 
     def test_all_disabled(self):
-        """With all corrections disabled, only phase space remains — still valid output."""
+        """With all corrections disabled, the core remains — still valid output."""
         config = SpectrumConfig(
             Z_parent=19,
             Z_daughter=20,
             A_number=40,
             endpoint_MeV=2.5,
-            use_fermi=False,
             use_screening=False,
             use_finite_size=False,
             use_charge_dist=False,
@@ -82,14 +80,14 @@ class TestBetaSpectrumToggles:
         assert len(values) > 0, "Should return non-empty array"
         assert np.all(np.isfinite(values)), "All values must be finite"
 
-    def test_fermi_only(self):
-        """With only fermi function enabled: should give F₀(Z,W) — not just phase space."""
+    def test_core_equals_phase_space_times_fermi(self):
+        """Core components are unconditional: with corrections disabled the
+        spectrum equals phase space × Fermi function exactly."""
         config_minimal = SpectrumConfig(
             Z_parent=19,
             Z_daughter=20,
             A_number=40,
             endpoint_MeV=2.5,
-            use_phase_space=False,
             use_screening=False,
             use_finite_size=False,
             use_charge_dist=False,
@@ -101,8 +99,11 @@ class TestBetaSpectrumToggles:
         W, _ = spectrum.get_energy_grid(config_minimal)
         values = spectrum(W)
 
-        # Fermi function alone is ~1 for Z=20 at mid-energy
-        assert np.all(np.isfinite(values)), "Fermi-only output must be finite"
+        components = spectrum.calculate_components(W)
+        assert set(components.keys()) == {"PhaseSpace", "Fermi"}
+        expected = components["PhaseSpace"] * components["Fermi"]
+        assert np.allclose(values, expected, rtol=1e-10)
+        assert np.all(np.isfinite(values)), "Core-only output must be finite"
 
 
 class TestBetaSpectrumComponentExtraction:
@@ -158,178 +159,3 @@ class TestBetaSpectrumRealisticConfig:
         assert np.all(np.isfinite(values)), "C-14 must not produce NaN/inf"
         # Phase space should peak in the bulk for low-Q decay
         assert np.any(values > 0), "Should have non-zero spectrum values"
-
-
-@pytest.mark.skip(reason="DetectorResponse not yet implemented")
-class TestDeclarativeDetectorResponse:
-    """Test declarative detector response workflow via SpectrumConfig."""
-
-    def test_create_detector_from_config(self):
-        """SpectrumConfig detector params should create valid DetectorResponse."""
-        config = SpectrumConfig(
-            Z_parent=43,
-            Z_daughter=44,
-            A_number=99,
-            endpoint_MeV=0.294,
-            detector_model="gaussian",
-            detector_sigma_a_keV=1.0,
-            detector_n_channels=512,
-        )
-
-        detector = BetaSpectrum.create_detector_from_config(config)
-
-        from beta_spectrum.components.detector_response import DetectorResponse
-        assert isinstance(detector, DetectorResponse)
-        assert detector.n_channels == 512
-        assert detector._mode == "analytical"
-        assert detector.model == "gaussian"
-
-    @pytest.mark.skip(reason="DetectorResponse not yet implemented")
-    def test_detector_sigma_conversion(self):
-        """Sigma in keV should be correctly converted to m_e units."""
-        from beta_spectrum.utils import T_to_W
-
-        config = SpectrumConfig(
-            Z_parent=43,
-            Z_daughter=44,
-            A_number=99,
-            endpoint_MeV=0.294,
-            detector_sigma_a_keV=1.0,
-            detector_fano_factor=0.0,
-        )
-
-        detector = BetaSpectrum.create_detector_from_config(config)
-
-        # 1.0 keV = 1.0/511.0 m_e ≈ 0.00196
-        expected_sigma = 1.0 / 0.510998950
-        sigma_at_endpoint = detector._resolution_sigma(T_to_W(0.294))
-        assert np.isclose(sigma_at_endpoint, expected_sigma, rtol=0.01)
-
-    @pytest.mark.skip(reason="DetectorResponse not yet implemented")
-    def test_convolve_detector_single_call(self):
-        """convolve_detector() should work in a single call with config params."""
-        config = SpectrumConfig(
-            Z_parent=43,
-            Z_daughter=44,
-            A_number=99,
-            endpoint_MeV=0.294,
-            detector_model="gaussian",
-            detector_sigma_a_keV=1.0,
-            detector_n_channels=512,
-            detector_fano_factor=0.0,
-        )
-
-        spectrum = BetaSpectrum.from_config(config)
-        W, _ = spectrum.get_energy_grid(config)
-
-        convolved = spectrum.convolve_detector(config, W=W)
-
-        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"
-
-    @pytest.mark.skip(reason="DetectorResponse not yet implemented")
-    def test_convolve_detector_with_tail(self):
-        """Convolution with tail model should produce low-energy tailing."""
-        config = SpectrumConfig(
-            Z_parent=43,
-            Z_daughter=44,
-            A_number=99,
-            endpoint_MeV=0.294,
-            detector_model="gaussian_tail",
-            detector_sigma_a_keV=1.0,
-            detector_tail_fraction=0.2,
-            detector_tau_keV=5.0,
-            detector_n_channels=512,
-            detector_fano_factor=0.0,
-        )
-
-        spectrum = BetaSpectrum.from_config(config)
-        W, _ = spectrum.get_energy_grid(config)
-
-        convolved = spectrum.convolve_detector(config, W=W)
-
-        assert convolved.shape == (512,)
-        assert np.all(convolved >= 0)
-
-    @pytest.mark.skip(reason="DetectorResponse not yet implemented")
-    def test_analyzer_convolved_from_config(self):
-        """BetaSpectrumAnalyzer should use config detector params automatically."""
-        config = SpectrumConfig(
-            Z_parent=43,
-            Z_daughter=44,
-            A_number=99,
-            endpoint_MeV=0.294,
-            use_detector_response=True,
-            detector_model="gaussian",
-            detector_sigma_a_keV=1.0,
-            detector_n_channels=256,
-            detector_fano_factor=0.0,
-        )
-
-        spectrum = BetaSpectrum.from_config(config)
-        from beta_spectrum import BetaSpectrumAnalyzer
-
-        analyzer = BetaSpectrumAnalyzer(spectrum, config)
-        detector = BetaSpectrum.create_detector_from_config(config)
-
-        convolved = analyzer.convolved_spectrum(detector, normalize=True)
-
-        assert convolved.shape == (256,)
-        # Normalized spectrum should integrate to ~1 over channel energies
-        integral = np.trapezoid(convolved, detector.channel_energies)
-        assert np.isclose(
-            integral, 1.0, rtol=0.05
-        ), f"Integral should be ~1.0, got {integral}"
-
-    @pytest.mark.skip(reason="DetectorResponse not yet implemented")
-    def test_detector_range_extends_to_endpoint(self):
-        """Detector channel range should extend to cover spectrum endpoint."""
-        config = SpectrumConfig(
-            Z_parent=43,
-            Z_daughter=44,
-            A_number=99,
-            endpoint_MeV=0.294,
-            detector_channel_energy_range=(0.0, 0.2),  # Less than endpoint
-            detector_n_channels=256,
-            detector_fano_factor=0.0,
-        )
-
-        detector = BetaSpectrum.create_detector_from_config(config)
-
-        # Range should extend beyond endpoint in m_e units
-        from beta_spectrum.utils import T_to_W
-
-        W0 = T_to_W(0.294)
-        assert detector.channel_energies[-1] >= W0, (
-            f"Detector range {detector.channel_energies[-1]:.4f} should cover "
-            f"endpoint {W0:.4f}"
-        )
-
-    @pytest.mark.skip(reason="DetectorResponse not yet implemented")
-    def test_different_detector_models(self):
-        """All detector models should produce valid convolved spectra."""
-        models = ["gaussian", "gaussian_tail", "tikhonov"]
-
-        for model in models:
-            config = SpectrumConfig(
-                Z_parent=43,
-                Z_daughter=44,
-                A_number=99,
-                endpoint_MeV=0.294,
-                detector_model=model,
-                detector_sigma_a_keV=1.0,
-                detector_n_channels=256,
-                detector_fano_factor=0.0,
-            )
-            if model == "gaussian_tail":
-                config.detector_tail_fraction = 0.1
-                config.detector_tau_keV = 5.0
-
-            spectrum = BetaSpectrum.from_config(config)
-            W, _ = spectrum.get_energy_grid(config)
-
-            convolved = spectrum.convolve_detector(config, W=W)
-
-            assert convolved.shape == (256,)
-            assert np.all(convolved >= 0), f"Model {model} produced negative values"

tests/physics/test_core_spectrum.py

@@ -0,0 +1,731 @@
+"""
+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(**over: Any) -> 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,
+        **over,
+        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_separate_branches_per_branch_unit_area(self, tmp_path: Any) -> None:
+        config = make_multi_config(separate_branches=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, match="non-positive"):
+            analyzer.export_to_csv(path)
+        assert not (tmp_path / "degenerate.csv").exists()
+
+    def test_degenerate_export_raises_in_separate_branches_mode(
+        self, tmp_path: Any
+    ) -> None:
+        """The positive-integral guard must apply in separate-branches mode too:
+        no all-zero CSV may be written for a degenerate spectrum."""
+        config = make_degenerate_config(separate_branches=True)
+        spectrum = BetaSpectrum.from_config(config)
+        analyzer = BetaSpectrumAnalyzer(spectrum, config)
+        path = str(tmp_path / "degenerate_to_one.csv")
+        with pytest.raises(ValueError, match="non-positive"):
+            analyzer.export_to_csv(path)
+        assert not (tmp_path / "degenerate_to_one.csv").exists()
+
+    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()
+
+
+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."""
+
+    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)
+
+
+class TestNormalizationModeMatrix:
+    """separate_branches and intensity modes must behave consistently across CSV
+    verbosity, branch count and the CSV-driven visualizer."""
+
+    def test_separate_branches_concise_multi_unit_area_no_total_no_components(
+        self, tmp_path: Any
+    ) -> None:
+        config = make_multi_config(separate_branches=True)
+        spectrum = BetaSpectrum.from_config(config)
+        analyzer = BetaSpectrumAnalyzer(spectrum, config)
+        path, df = export_and_read(
+            analyzer, tmp_path, "to_one_concise.csv", csv_verbosity="concise"
+        )
+
+        assert set(df.columns) == {
+            "energy_MeV",
+            "branch_1",
+            "branch_2",
+            "branch_3",
+        }
+        header = header_lines_of(path)
+        assert any("branch normalization: separate" in h for h in header)
+        energies = df["energy_MeV"].to_numpy()
+        for col in ("branch_1", "branch_2", "branch_3"):
+            area = np.trapezoid(df[col].to_numpy(), energies)
+            assert np.isclose(area, 1.0, rtol=5e-4)
+
+    def test_separate_branches_single_mode_unit_area(self, tmp_path: Any) -> None:
+        config = make_single_config(separate_branches=True)
+        spectrum = BetaSpectrum.from_config(config)
+        analyzer = BetaSpectrumAnalyzer(spectrum, config)
+        _path, df = export_and_read(analyzer, tmp_path, "to_one_single.csv")
+
+        assert "total_spectrum" not in df.columns
+        area = np.trapezoid(
+            df["branch_1_spectrum"].to_numpy(), df["energy_MeV"].to_numpy()
+        )
+        assert np.isclose(area, 1.0, rtol=5e-4)
+
+    def test_intensity_concise_branch_columns_sum_to_total(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, "intensity_concise.csv", csv_verbosity="concise"
+        )
+
+        summed = (
+            df["branch_1"].to_numpy()
+            + df["branch_2"].to_numpy()
+            + df["branch_3"].to_numpy()
+        )
+        total = df["total_spectrum"].to_numpy()
+        nonzero = total > 1e-12
+        assert np.all(np.abs(summed[nonzero] - total[nonzero]) <= total[nonzero] * 5e-4)
+
+    def test_plot_from_csv_separate_branches_plots_branches_without_total(
+        self, tmp_path: Any
+    ) -> None:
+        config = make_multi_config(separate_branches=True)
+        spectrum = BetaSpectrum.from_config(config)
+        analyzer = BetaSpectrumAnalyzer(spectrum, config)
+        path, df = export_and_read(analyzer, tmp_path, "to_one_plot.csv")
+
+        fig = Visualizer.plot_from_csv(path, save_path=None, show_components=True)
+        assert len(fig.axes) == 4
+        # Panel 1 has exactly the 3 branch curves — no total spectrum
+        # exists in separate-branches mode.
+        panel = fig.axes[0]
+        curves = [np.asarray(ln.get_xydata()) for ln in panel.lines]
+        assert len(curves) == 3
+        energies = df["energy_MeV"].to_numpy()
+        assert np.allclose(curves[0][:, 0], energies)
+        assert np.allclose(curves[0][:, 1], df["branch_1_spectrum"].to_numpy())
+
+    def test_plot_analysis_separate_branches_returns_figure(
+        self, tmp_path: Any
+    ) -> None:
+        config = make_multi_config(separate_branches=True)
+        spectrum = BetaSpectrum.from_config(config)
+        analyzer = BetaSpectrumAnalyzer(spectrum, config)
+        save = str(tmp_path / "analysis.png")
+        fig = Visualizer.plot_analysis(analyzer, save_path=save)
+        assert len(fig.axes) == 4
+        assert (tmp_path / "analysis.png").exists()
+
+    def test_plot_from_csv_separate_branches_concise_warns_and_falls_back(
+        self, tmp_path: Any, capsys: Any
+    ) -> None:
+        config = make_multi_config(separate_branches=True)
+        spectrum = BetaSpectrum.from_config(config)
+        analyzer = BetaSpectrumAnalyzer(spectrum, config)
+        path, df = export_and_read(
+            analyzer, tmp_path, "to_one_concise2.csv", csv_verbosity="concise"
+        )
+
+        fig = Visualizer.plot_from_csv(path, save_path=None, show_components=True)
+        captured = capsys.readouterr()
+        assert "Warning" in captured.out
+        # Simple spectra view: one axis, 3 branch curves, no total curve.
+        assert len(fig.axes) == 1
+        assert len(fig.axes[0].lines) == 3
+
+
+class TestSingleModeColumnCleanup:
+    """Single-mode CSVs carry the factorized components once, unprefixed;
+    multi-mode CSVs carry them per branch; the visualizer reads whichever
+    naming the file actually uses."""
+
+    def test_single_mode_full_export_has_unprefixed_components_only(
+        self, tmp_path: Any
+    ) -> None:
+        config = make_single_config(use_radiative=True)
+        spectrum = BetaSpectrum.from_config(config)
+        analyzer = BetaSpectrumAnalyzer(spectrum, config)
+        path, df = export_and_read(analyzer, tmp_path, "single_components.csv")
+
+        assert "PhaseSpace" in df.columns
+        assert "Radiative" in df.columns
+        assert "Fermi" in df.columns
+        assert "branch_1_spectrum" in df.columns
+        assert "branch_1_PhaseSpace" not in df.columns
+        assert "branch_1_Radiative" not in df.columns
+
+    def test_multi_mode_full_export_has_prefixed_core_factors(
+        self, tmp_path: Any
+    ) -> None:
+        config = make_multi_config(use_radiative=True)
+        spectrum = BetaSpectrum.from_config(config)
+        analyzer = BetaSpectrumAnalyzer(spectrum, config)
+        path, df = export_and_read(analyzer, tmp_path, "multi_components.csv")
+
+        assert "branch_1_PhaseSpace" in df.columns
+        assert "branch_2_Radiative" in df.columns
+        assert "Fermi" in df.columns
+        # The bare core-factor names belong to single mode only.
+        assert "PhaseSpace" not in df.columns
+        assert "Radiative" not in df.columns
+
+    def test_debug_view_single_csv_includes_core_panels(self, tmp_path: Any) -> None:
+        config = make_single_config(use_radiative=True)
+        spectrum = BetaSpectrum.from_config(config)
+        analyzer = BetaSpectrumAnalyzer(spectrum, config)
+        path, _df = export_and_read(analyzer, tmp_path, "single_plot.csv")
+
+        fig = Visualizer.plot_from_csv(path, None, show_components=True)
+        assert len(fig.axes) == 4
+        # Panel 2: Fermi + the single branch's phase space.
+        assert len(fig.axes[1].lines) == 2
+        labels = [ln.get_label() for ln in fig.axes[1].lines]
+        assert "Fermi" in labels
+        assert any("PS Br" in str(lbl) for lbl in labels)
+        # Panel 4: the radiative correction for the single branch.
+        assert len(fig.axes[3].lines) == 1
+
+    def test_ylabel_uses_per_MeV_units(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, "ylabel.csv")
+
+        fig = Visualizer.plot_from_csv(path, None, show_components=True)
+        assert "MeV$^{-1}$" in fig.axes[0].get_ylabel()

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"

tests/physics/test_multi_branch.py

@@ -21,8 +21,7 @@ We test:
  6. **Component extraction**: universal + per-branch components correctly separated.
  7. **CSV export**: multi-branch CSV has correct column structure.
  8. **Energy grid**: uses max endpoint across all branches.
- 9. **Detector response**: multi-branch spectra are convolvable.
- 10. **BranchConfig dataclass**: basic construction and validation.
+ 9. **BranchConfig dataclass**: basic construction and validation.
 """
 
 import numpy as np
@@ -249,8 +248,8 @@ class TestMultiBranchSpectrum:
             Z_daughter=44,
             A_number=99,
             endpoint_MeV=0.350,
-            # Disable universal components to simplify
-            use_fermi=False,
+            # Disable higher-order universal corrections to simplify:
+            # the universal factor reduces to the core Fermi function.
             use_screening=False,
             use_finite_size=False,
             use_charge_dist=False,
@@ -266,8 +265,11 @@ class TestMultiBranchSpectrum:
         weighted = spectrum.get_branch_normalized_spectra(W)
         assert len(weighted) == 2
 
-        # With no universal components, total = sum of weighted branch spectra
-        total_from_branches = weighted[0] + weighted[1]
+        # Universal factor is the core Fermi function only (corrections off):
+        # total = fermi × sum(weighted branch spectra)
+        components = spectrum.calculate_components(W)
+        fermi = components["Fermi"]
+        total_from_branches = fermi * (weighted[0] + weighted[1])
         total = spectrum(W)
         assert np.allclose(total, total_from_branches, rtol=1e-10)
 
@@ -278,8 +280,6 @@ class TestMultiBranchSpectrum:
             Z_daughter=44,
             A_number=99,
             endpoint_MeV=0.350,
-            use_phase_space=True,
-            use_fermi=False,
             use_screening=False,
             use_finite_size=False,
             use_charge_dist=False,
@@ -293,13 +293,16 @@ class TestMultiBranchSpectrum:
         spectrum = BetaSpectrum.from_config(config)
         W, _ = spectrum.get_energy_grid(config)
 
-        # Both branches have same endpoint and transition type,
-        # so phase space should be identical. Total = 0.8*PS + 0.2*PS = PS
+        # Both branches have same endpoint and transition type, so their
+        # per-branch cores are identical: sum = 0.8*core + 0.2*core = core.
+        # Total = universal (Fermi) × sum.
         total = spectrum(W)
         branch0 = spectrum.branch_spectra[0](W)
         branch1 = spectrum.branch_spectra[1](W)
+        components = spectrum.calculate_components(W)
+        fermi = components["Fermi"]
 
-        expected = 0.8 * branch0 + 0.2 * branch1
+        expected = fermi * (0.8 * branch0 + 0.2 * branch1)
         assert np.allclose(total, expected, rtol=1e-10)
 
 
@@ -371,8 +374,8 @@ class TestMultiBranchComponents:
         assert "Fermi" in components
         assert "branch_0.PhaseSpace" in components
 
-    def test_analyzer_branch_spectra_property(self):
-        """Analyzer should expose branch spectra property."""
+    def test_analyzer_physical_branch_spectra(self):
+        """Analyzer physical branch spectra are core × corrections arrays."""
         config = SpectrumConfig(
             Z_parent=43,
             Z_daughter=44,
@@ -386,8 +389,11 @@ class TestMultiBranchComponents:
         spectrum = BetaSpectrum.from_config(config)
         analyzer = BetaSpectrumAnalyzer(spectrum, config)
 
-        branch_spectra = analyzer.branch_spectra
-        assert len(branch_spectra) == 2
+        physical = [
+            analyzer.physical_branch_spectrum(i) for i in range(analyzer._n_branches)
+        ]
+        assert len(physical) == 2
+        assert physical[0].shape == analyzer.W.shape
 
     def test_total_equals_universal_times_branch_sum(self):
         """Total = universal_components * sum(intensity * branch_spectra)."""
@@ -481,6 +487,13 @@ class TestMultiBranchCSVExport:
         # Should have universal columns
         assert "Fermi" in df.columns
 
+        # Branch columns are physical spectra: their sum equals total_spectrum
+        total = df["total_spectrum"].to_numpy()
+        branch_sum = (
+            df["branch_1_spectrum"].to_numpy() + df["branch_2_spectrum"].to_numpy()
+        )
+        assert np.allclose(branch_sum, total, rtol=1e-3)
+
     def test_single_branch_csv_columns(self, tmp_path):
         """Single-branch CSV should have standard columns."""
         config = SpectrumConfig(
@@ -582,42 +595,6 @@ class TestMultiBranchTransitionTypes:
         assert np.all(values >= 0)
 
 
-# ---------------------------------------------------------------------------
-# Multi-branch with detector response
-# ---------------------------------------------------------------------------
-
-
-class TestMultiBranchDetectorResponse:
-    """Test multi-branch with detector response convolution."""
-
-    def test_multi_branch_with_detector(self):
-        """Multi-branch spectrum should be convolvable with detector response."""
-        config = SpectrumConfig(
-            Z_parent=43,
-            Z_daughter=44,
-            A_number=99,
-            endpoint_MeV=0.350,
-            use_detector_response=True,
-            detector_model="gaussian",
-            detector_sigma_a_keV=1.0,
-            detector_n_channels=256,
-            detector_fano_factor=0.0,
-            branches=[
-                BranchConfig(endpoint_MeV=0.294, transition_type="F2", intensity=0.7),
-                BranchConfig(endpoint_MeV=0.150, transition_type="A", intensity=0.3),
-            ],
-        )
-        spectrum = BetaSpectrum.from_config(config)
-        W, _ = spectrum.get_energy_grid(config)
-
-        detector = BetaSpectrum.create_detector_from_config(config)
-        convolved = spectrum.convolve_with_detector(detector, W=W, config=config)
-
-        assert convolved.shape == (256,)
-        assert np.all(convolved >= 0), "Convolved spectrum should be non-negative"
-        assert np.sum(convolved) > 0, "Convolved spectrum should have non-zero area"
-
-
 # ---------------------------------------------------------------------------
 # Edge cases
 # ---------------------------------------------------------------------------
@@ -705,6 +682,7 @@ class TestMultiBranchEdgeCases:
                 ),
             ],
         )
+        assert config.branches is not None
         spectrum = BetaSpectrum.from_config(config)
         W, energies = spectrum.get_energy_grid(config)
 
@@ -748,6 +726,7 @@ class TestMultiBranchEdgeCases:
                 BranchConfig(endpoint_MeV=0.3, transition_type="F1", intensity=0.2),
             ],
         )
+        assert config.branches is not None
         spectrum = BetaSpectrum.from_config(config)
         W, energies = spectrum.get_energy_grid(config)
 
@@ -775,6 +754,7 @@ class TestMultiBranchEdgeCases:
         # Each branch integral should be proportional to its intensity
         # (within numerical tolerance, since different endpoints give different
         # phase space integrals)
+        assert config.branches is not None
         for i, (integral, intensity) in enumerate(zip(integrals, config.branches)):
             assert integral > 0, f"Branch {i} has zero integral"
 

tests/physics/test_neutrino_mass_effects.py

@@ -118,10 +118,15 @@ class TestNeutrinoMassCLIIntegration:
         assert args.neutrino_mass == 0.01, "CLI should parse neutrino_mass correctly"
 
     def test_cli_neutrino_mass_default_value(self):
-        """Test that CLI neutrino mass parameter defaults to 0.0."""
+        """Unset flag parses as None ("not provided"); the effective
+        default (0.0 MeV) is applied when the config is built."""
         from beta_spectrum.cli import _build_parser
 
         parser = _build_parser()
         args = parser.parse_args(["--nuclide", "Tc99"])
 
-        assert args.neutrino_mass == 0.0, "CLI neutrino_mass default should be 0.0"
+        assert args.neutrino_mass is None, (
+            "unset --neutrino-mass must parse as None so JSON file values "
+            "are not overwritten"
+        )
+        assert hasattr(args, "neutrino_mass")

tests/quality/test_api_detector_response.py

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

tests/quality/test_cli_options.py

@@ -0,0 +1,137 @@
+"""CLI option surface tests: renamed branch-separation flag, short
+aliases for every long option, and dry-run output shape (audit follow-up
+by the maintainer: rename --normalize-to-one, add short options, drop the
+redundant top-level Transition line)."""
+
+from __future__ import annotations
+
+import dataclasses
+from typing import Any, Optional
+
+import pytest
+
+from beta_spectrum import BranchConfig, SpectrumConfig
+from beta_spectrum.cli import _build_parser, _print_dry_run_output
+
+pytestmark = pytest.mark.quality
+
+
+def make_config() -> SpectrumConfig:
+    return SpectrumConfig(
+        Z_parent=59,
+        Z_daughter=60,
+        A_number=144,
+        endpoint_MeV=2.98,
+        transition_type="F1",
+        e_step_MeV=0.001,
+        branches=[
+            BranchConfig(endpoint_MeV=2.98, transition_type="F1", intensity=0.9),
+            BranchConfig(endpoint_MeV=2.1, transition_type="A", intensity=0.1),
+        ],
+    )
+
+
+class TestSeparateBranchesRename:
+    """--normalize-to-one was renamed to --separate-branches (no alias)."""
+
+    def test_flag_parses_true(self) -> None:
+        parser = _build_parser()
+        args = parser.parse_args(["input.json", "--separate-branches"])
+        assert args.separate_branches is True
+
+    def test_flag_default_is_none(self) -> None:
+        parser = _build_parser()
+        args = parser.parse_args(["input.json"])
+        assert args.separate_branches is None
+
+    def test_old_flag_is_rejected(self) -> None:
+        parser = _build_parser()
+        with pytest.raises(SystemExit):
+            parser.parse_args(["input.json", "--normalize-to-one"])
+
+    def test_config_field_renamed(self) -> None:
+        names = {f.name for f in dataclasses.fields(SpectrumConfig)}
+        assert "separate_branches" in names
+        assert "normalize_to_one" not in names
+
+    def test_config_old_kwarg_rejected(self) -> None:
+        kwargs: dict[str, Any] = {
+            "Z_parent": 15,
+            "Z_daughter": 16,
+            "A_number": 32,
+            "endpoint_MeV": 1.7103,
+            "transition_type": "A",
+            "normalize_to_one": True,
+        }
+        with pytest.raises(TypeError):
+            SpectrumConfig(**kwargs)
+
+
+# (short, long, value-or-None-for-store_true)
+_SHORT_EQUIV = [
+    ("-N", "--nuclide", "Cs137"),
+    ("-i", "--input", "in.json"),
+    ("-s", "--e-step", "0.01"),
+    ("-m", "--neutrino-mass", "0.15"),
+    ("-o", "--output", "out.csv"),
+    ("-p", "--plot", "plot.png"),
+    ("-L", "--log-file", "calc.log"),
+    ("-d", "--dry-run", None),
+    ("-c", "--intensity-cutoff", "0.5"),
+    ("-V", "--csv-verbosity", "concise"),
+    ("-S", "--csv-separator", "tab"),
+    ("-b", "--separate-branches", None),
+]
+
+
+class TestShortOptions:
+    """Every long option gains a single-letter alias (only -v/-q existed)."""
+
+    @pytest.mark.parametrize(
+        "short,long,value",
+        _SHORT_EQUIV,
+        ids=[s.lstrip("-") for s, _, _ in _SHORT_EQUIV],
+    )
+    def test_short_matches_long(
+        self, short: str, long: str, value: Optional[str]
+    ) -> None:
+        parser = _build_parser()
+        argv_long = [long] if value is None else [long, value]
+        argv_short = [short] if value is None else [short, value]
+        args_long = parser.parse_args(["Tc99"] + argv_long)
+        args_short = parser.parse_args(["Tc99"] + argv_short)
+        dest = long.lstrip("-").replace("-", "_")
+        assert getattr(args_short, dest) == getattr(args_long, dest)
+        assert getattr(args_long, dest) is not None or value is None
+
+    def test_existing_shorts_unchanged(self) -> None:
+        parser = _build_parser()
+        args = parser.parse_args(["Tc99", "-v", "-q"])
+        assert args.verbose == 1
+        assert args.quiet is True
+
+
+class TestDryRunOutputShape:
+    """The top-level Transition line is redundant with the per-branch
+    list (for multi-branch runs) and was removed."""
+
+    def test_multi_branch_dry_run_no_top_level_transition(self, capsys: Any) -> None:
+        _print_dry_run_output(make_config(), "paceENSDF")
+        out = capsys.readouterr().out
+        assert "Transition:" not in out
+        assert "transition=F1" in out
+        assert "transition=A" in out
+        assert "Multi-branch mode: 2 branches" in out
+
+    def test_single_branch_dry_run_reports_endpoint(self, capsys: Any) -> None:
+        config = SpectrumConfig(
+            Z_parent=15,
+            Z_daughter=16,
+            A_number=32,
+            endpoint_MeV=1.7103,
+            transition_type="A",
+        )
+        _print_dry_run_output(config, "json")
+        out = capsys.readouterr().out
+        assert "Transition:" not in out
+        assert "Mode: single-branch" in out

tests/quality/test_logging_utils.py

@@ -154,7 +154,7 @@ class TestCSVMetadataHeader:
             lines = f.readlines()
 
         assert any("# beta-spectrum v" in line for line in lines)
-        assert any("0.4.0" in line for line in lines)
+        assert any("0.4.1" in line for line in lines)
         os.unlink(filename)
 
     def test_header_contains_timestamp(self):

tests/quality/test_neutrino_mass_cli.py

@@ -76,10 +76,14 @@ class TestCLIArguments:
         assert args.neutrino_mass == 0.01, "CLI neutrino_mass should be 0.01"
 
     def test_neutrino_mass_cli_default_value(self):
-        """Test that --neutrino-mass defaults to 0.0 in CLI."""
+        """Unset flag parses as None ("not provided") so JSON file
+        values are not overwritten; the effective default (0.0 MeV) is
+        applied when the config is built."""
         from beta_spectrum.cli import _build_parser
 
         parser = _build_parser()
 
         args = parser.parse_args(["--nuclide", "Tc99"])
-        assert args.neutrino_mass == 0.0, "CLI neutrino_mass default should be 0.0"
+        assert (
+            args.neutrino_mass is None
+        ), "unset --neutrino-mass must parse as None (not a value)"

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
@@ -160,9 +161,12 @@ class TestJsonInput:
         path = self._write_json(data)
         try:
             loaded = load_json_input(str(path))
-            # Defaults should be filled in
-            assert loaded["detector_model"] == "gaussian"
-            assert loaded["detector_sigma_a_keV"] == 1.0
+            # Defaults should be filled in. The detector_* schema keys were
+            # removed together with the detector_response module; JSON
+            # defaults carry spectrum physics only.
+            assert "detector_model" not in loaded
+            assert loaded["neutrino_mass_MeV"] == 0.0
+            assert loaded["transition_type"] == "F2"
         finally:
             path.unlink()
 
@@ -241,12 +245,96 @@ 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
 # ---------------------------------------------------------------------------
 
 
+@pytest.mark.slow
+@pytest.mark.requires_ensdf
+@pytest.mark.skipif(not _has_paceENSDF(), reason="paceENSDF not installed")
+class TestPaceENSDFBranchEndpoints:
+    """Per-branch endpoints are the measured transition endpoint energies
+    (ENSDF column 7), not base-endpoint minus level energy subtractions
+    with a base taken from an arbitrary first branch (D3)."""
+
+    def _branch_endpoints(self, nuclide: str) -> list[float]:
+        info = get_decay_info_from_paceENSDF(nuclide, "beta_minus")
+        cfg = decay_info_to_config(info)
+        assert cfg.branches is not None
+        return sorted(
+            (b.endpoint_MeV for b in cfg.branches),
+            reverse=True,
+        )
+
+    def test_co60_branch_endpoints_positive(self) -> None:
+        # Measured endpoints: 1492.0, 670.0, 317.88 keV (levels above the
+        # first branch's endpoint previously produced negative values and
+        # a hard ValueError from BranchConfig validation).
+        eps = self._branch_endpoints("Co60")
+        assert eps == pytest.approx([1.492, 0.670, 0.31788], abs=1e-6)
+        assert all(e > 0 for e in eps)
+
+    def test_tc99_branch_endpoints(self) -> None:
+        eps = self._branch_endpoints("Tc99")
+        assert eps == pytest.approx([0.2975, 0.2080], abs=1e-6)
+
+    def test_cs137_branch_endpoints(self) -> None:
+        eps = self._branch_endpoints("Cs137")
+        assert eps == pytest.approx([1.176, 0.89213, 0.51403], abs=1e-6)
+
+    def test_overall_endpoint_is_max_branch_endpoint(self) -> None:
+        info = get_decay_info_from_paceENSDF("Co60", "beta_minus")
+        assert info.endpoint_MeV == pytest.approx(1.492, abs=1e-6)
+
+
 @pytest.mark.requires_ensdf
 @pytest.mark.skipif(not _has_paceENSDF(), reason="paceENSDF not installed")
 class TestPaceENSDFIntegration: