style: apply black formatting after cleanup

Remove gA-extraction notes, nuclear-structure notes, paper-summary
refs, detector-response correction note and the finished test
migration checklist. Rewrite the vault index for the single-package
theory calculator; repoint dangling wiki-links in surviving notes to
inline citations; refresh docs_check.py (drop bs_exp CLI check,
register removed files).

Literature and nuclear-structure research resumes in dedicated
projects; the removed notes are recoverable from history.

.gitignore

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

AGENTS.md

@@ -1,6 +1,6 @@
 # Agent Instructions — beta-spectrum
 
-High-precision beta decay spectrum calculator. Extracts shape factor C(W) from ⁹⁹Tc beta spectrum to parametrize g_V coupling.
+High-precision beta decay spectrum shape calculator. Since v0.5.0 the repository is dedicated solely to spectrum shape calculation: detector response, calibration, fitting and C(W)/g_V/g_A extraction live in separate, dedicated projects.
 
 ## Stack & Commands
 
@@ -68,7 +68,6 @@ After quality gates pass, run:
 
 ```bash
 bs_pnpi --nuclide Tc99 -vv --output ./output/verification.csv --plot ./output/verification.png --log-file ./output/debug.log 2>&1 | tee ./output/debug_stdout.log
-bs_exp --help                    # verify experimental data CLI
 ```
 
 Verify parameter consistency across all components. Save artifacts to `./output/`.

CONVENTIONS.md

@@ -18,7 +18,7 @@ m_e = ℏ = c = 1
 
 **External API:** human-readable MeV/keV. Check the unit suffix in parameter names:
 - `_MeV` suffix → value is in MeV
-- `_keV` suffix → value is in keV (detector params only)
+- `_keV` suffix → value is in keV
 - No suffix on energy → value is in m_e units (internal)
 
 ### Conversion

README.md

@@ -2,6 +2,8 @@
 
 A modular Python package for calculating high-precision **beta decay energy spectra**. It models the spectral shape of nuclear transitions by combining well-established physical corrections into a coherent, configurable pipeline.
 
+Since v0.5.0 the project scope is the spectrum-shape **calculation** alone: detector response modeling, data calibration and fitting/extraction live in separate, dedicated projects.
+
 **Features:**
 
 - Numerical evaluation from near-zero kinetic energies up to the endpoint
@@ -11,6 +13,8 @@ A modular Python package for calculating high-precision **beta decay energy spec
   - Normal mode: single spectrum plot with nuclear data header
   - Debug mode (`-vv`): 4-panel view with all correction factors
   - All plots include commit hash and UTC timestamp for traceability
+- Multi-branch decays with intensity weighting and per-branch endpoints
+- Nuclear data retrieval from ENSDF via paceENSDF (`bs_pnpi --nuclide Tc99`)
 - Comprehensive scientific reference in [`docs/`](docs/)
 - Debug verification workflow: `bs_pnpi -vv` validates all parameters end-to-end
 
@@ -22,7 +26,7 @@ The beta spectrum is constructed as a product of independent correction factors:
 
 $$N(W) \\propto p W (W_0 - W)^2 \\times \\prod_i C_i(W)$$
 
-Each $C_i(W)$ is an independently implemented and toggleable component.
+Each $C_i(W)$ is an independently implemented and toggleable component. The phase space and the Fermi function are the unconditional core: every spectrum carries their product, and the higher-order corrections multiply on top of it.
 
 ______________________________________________________________________
 
@@ -31,68 +35,46 @@ ______________________________________________________________________
 ```
 beta_shape_pnpi/
 ├── pyproject.toml
-├── docs/                          # Obsidian knowledge base (physics theory & references)
-│   ├── README.md                  # Vault overview + navigation index
-│   ├── 00-beta-spectrum-overview.md    # Master equation, all corrections table
-│   ├── 01-fermi-function.md           # Fermi function F₀(Z,W)
-│   ├── 02-phase-space.md              # Phase space p·W·(W₀−W)²
-│   ├── 03-finite-size.md              # L₀, U, DFS corrections
-│   ├── 04-screening-correction.md     # Atomic screening S(Z,W)
-│   ├── 05-exchange-correction.md      # Electron exchange X(Z,W)
-│   ├── 06-radiative-corrections.md    # Outer radiative δᵣ(W,W₀)
-│   ├── 07-recoil-effects.md           # Weak magnetism, recoil kinematics
-│   ├── 08-detector-response.md        # Detector response models
-│   ├── 10-nuclear-structure.md        # Shape factor C(Z,W), BB vs HS formalisms
-│   ├── 12-atomic-overlap.md           # Bahcall correction r(Z,W)
-│   ├── 13-chemical-effects.md         # Molecular environment effects
-├── beta_spectrum/           ← pure theory core
-│   ├── __init__.py          # Public API: constants, utilities, classes
-│   ├── base.py              # Abstract SpectrumComponent base class
-│   ├── constants.py         # Physical constants (natural units, m_e = 1)
-│   ├── utils.py             # Helpers: T<->W conversion, momentum, nuclear radius
-│   ├── spectrum.py          # BetaSpectrum + SpectrumConfig + BranchConfig
-│   ├── nuclear_data.py      # paceENSDF integration + JSON input support
-│   ├── cli.py               # Command-line interface (bs_pnpi)
-│   ├── logging_utils.py     # Structured logging helpers
-│   └── components/
-│       │   ├── phase_space.py           ✓ Phase space shape (p·W·(W₀−W)²)
-│       │   ├── fermi.py                 ✓ Coulomb correction (loggamma for stability)
-│       │   ├── 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
-├── beta_spectrum/output/    ← output, plotting, visualization, CSV export
-│   ├── __init__.py          # BetaSpectrumAnalyzer, CSVEmitter, Visualizer
-│   ├── core.py              # BetaSpectrumAnalyzer base class
-│   ├── csv_export.py        # CSV export functionality
-│   └── visualization.py     # Plotting and visualization functionality
-├── exp_data/                ← experimental data handling
-│   ├── __init__.py          # ExpSpectrum, EnergyCalibrator, GaussianFitter, etc.
-│   ├── spectrum.py          # ExpSpectrum dataclass (energies, counts, errors, metadata)
-│   ├── calibration.py       # EnergyCalibrator (linear/quadratic channel→energy)
-│   ├── fitters.py           # GaussianFitter + PeakFitter for cal sources
-│   └── corrections.py       # DeadTimeCorrection, PileUpCorrection, BackgroundSubtractor
-├── fitter/                  ← analysis framework
-│   ├── __init__.py          # SpectrumModel, SpectrumFitter, CWExtractor, etc.
-│   ├── model.py             # SpectrumModel (BetaSpectrum + DetectorResponse wrapper)
-│   ├── fit_engine.py        # SpectrumFitter + CurveFitter + FitConfig + FitResult
-│   ├── extractor.py         # CWExtractor, GVAExtractor, CWExtractionResult, GVAExtractionResult
-│   └── result.py            # AnalysisResult + FitSummary
-├── tests/
-│   ├── exp_data/            # exp_data module tests
-│   ├── fitter/              # fitter module tests
-│   ├── physics/             # Physics component tests
-│   ├── quality/             # Code quality tests (logging, API, nuclear data)
-│   ├── integration/         # Full pipeline integration tests
-│   └── common/              # Shared hypothesis + property tests
+├── docs/                          ← Obsidian knowledge base (physics theory)
+│   ├── README.md                  ← Vault overview + navigation index
+│   ├── 00-beta-spectrum-overview.md
+│   ├── 02-phase-space.md
+│   └── corrections/               ← One note per correction component
+├── beta_spectrum/                 ← theory-only core package
+│   ├── __init__.py                ← Public API
+│   ├── base.py                    ← Abstract SpectrumComponent base class
+│   ├── constants.py               ← Physical constants (natural units, m_e = 1)
+│   ├── utils.py                   ← T<->W conversion, momentum, nuclear radius
+│   ├── spectrum.py                ← BetaSpectrum + SpectrumConfig + BranchConfig
+│   ├── nuclear_data.py            ← paceENSDF integration + JSON input support
+│   ├── cli.py                     ← Command-line interface (bs_pnpi)
+│   ├── logging_utils.py           ← Structured logging helpers
+│   ├── components/
+│   │   ├── phase_space.py         ← Phase space shape (p·W·(W₀−W)²)
+│   │   ├── fermi.py               ← Coulomb correction (loggamma for stability)
+│   │   ├── 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
+│   ├── data/
+│   │   └── exchange_coeff.csv     ← Tabulated coefficients for X(Z,W), Z=2..120
+│   └── output/                    ← visualization, CSV export
+│       ├── core.py                ← BetaSpectrumAnalyzer
+│       ├── csv_export.py          ← CSV export
+│       └── visualization.py       ← Plotting built from exported CSV
 ├── data/
-│   ├── exchange_coeff.csv         # Tabulated coefficients for X(Z,W), Z=2..120
-│   └── custom_input_example.json  # Sample JSON input file
-├── output/                      # Test run outputs (gitignored except .gitkeep)
-├── AGENTS.md                    # Agent instructions: stack, git workflow, TDD, quality gates
-├── CONVENTIONS.md               # Physics conventions: natural units, energy representation, component design
-├── TODO.md                      # Development roadmap
+│   └── custom_input_example.json  ← Sample JSON input file
+├── notebooks/                     ← Calculator notebook + nbmake conftest
+├── output/                        ← Test run outputs (gitignored except .gitkeep)
+├── tests/
+│   ├── physics/                   ← Physics component tests
+│   ├── quality/                   ← Code quality tests (logging, API, nuclear data)
+│   ├── integration/               ← CLI and pipeline integration tests
+│   └── common/                    ← Shared hypothesis + property tests
+├── scripts/                       ← docs_check.py, notebook runners
+├── AGENTS.md                      ← Agent instructions: stack, git workflow, TDD, quality gates
+├── CONVENTIONS.md                 ← Physics conventions: natural units, energy representation, component design
+├── TODO.md                        ← Development roadmap
 ```
 
 ______________________________________________________________________
@@ -122,13 +104,20 @@ 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
-- `get_data()` — return all numerical data for custom analysis
+  - 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 downstream analysis
 
 ______________________________________________________________________
 
@@ -156,27 +145,21 @@ analyzer.plot_analysis("debug.png", show_components=True)  # Debug mode (4-panel
 analyzer.export_to_csv("spectrum.csv")   # Export data to CSV
 ```
 
-### Fitting and C(W) extraction
+### From the command line
 
-```python
-from beta_spectrum import SpectrumConfig, BetaSpectrum
-from fitter import SpectrumModel, SpectrumFitter, CWExtractor
-from exp_data.spectrum import ExpSpectrum
-
-# Build theoretical model
-spectrum = BetaSpectrum.from_config(config)
-model = SpectrumModel(spectrum, detector_response=detector)
-
-# Fit to experimental data
-exp = ExpSpectrum(energies=energies_keV, counts=counts, errors=errors)
-fitter = SpectrumFitter(model, exp.energies, exp.counts, exp.errors)
-result = fitter.fit(x0=[1.0, 0.0])
-
-# Extract C(W) shape factor
-extractor = CWExtractor(config)
-cw_result = extractor.extract(measured, errors, energies_keV, model_values, endpoint_keV=Q)
+```bash
+bs_pnpi --nuclide Tc99 -vv --output spectrum.csv   # ENSDF lookup + debug export
+bs_pnpi --input data/custom_input_example.json     # JSON-driven configuration
 ```
 
+CLI flags (`--e-step`, `--neutrino-mass`, `--intensity-cutoff`, `--csv-separator`, `--separate-branches`) override JSON input values only when explicitly given; every long option has a single-letter alias.
+
+______________________________________________________________________
+
+## Nuclear Data (paceENSDF)
+
+`beta_spectrum.nuclear_data` retrieves decay data automatically: `get_decay_info_from_paceENSDF()`, `decay_info_to_config()`, and `create_config_from_source("paceENSDF", nuclide="Tc99")`. Z, A, endpoint energies, half-life, spin/parity, branch intensities, log*ft and forbiddenness are populated from ENSDF; per-branch endpoints are the measured transition endpoint energies (max endpoint across branches defines the spectrum endpoint). paceENSDF is an optional extra: `pip install ".[ensdf]"`.
+
 ______________________________________________________________________
 
 ## Implementation Status
@@ -184,9 +167,6 @@ ______________________________________________________________________
 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
 - **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` with full paceENSDF integration, structured logging (`-v`/`-vv`/`-q`), `--dry-run`, `--version`, `--log-file`
+- **CSV metadata headers** with nuclide, decay and calculation parameters plus software version and git commit hash

TODO.md

@@ -1,298 +1,119 @@
 # Strategic Development Roadmap
 
-## Scientific Goal
+## Scientific Scope
 
-Extract the shape factor C(W) from the experimentally measured beta spectrum of ⁹⁹Tc, and parametrize C(W) via the vector coupling constant g_V from the theory of beta decay.
+High-precision calculation of beta-decay spectrum shapes. Given a nuclide, decay data and a set of correction toggles, the calculator produces C(W)-level spectral shapes accurate to a few parts in 10⁻⁴ down to 1 keV for low-to-medium Z nuclei (Hayen et al. approach).
 
-This requires repeating the analysis pipeline of [Paulsen et al., Phys. Rev. C **110**, 055503 (2024)](https://doi.org/10.1103/PhysRevC/110/055503), which used Metallic Magnetic Calorimeters (MMCs) to measure the ⁹⁹Tc beta spectrum with sub-keV threshold and ~100 eV resolution.
-
-See `docs/refs/2024_Paulsen_T99_beta_spectrum/summary.md` for full analysis details.
+For forbidden nonunique decays the shape factor C(W) enters the calculated spectrum as an **input**. Extracting C(W) from measured spectra and parametrizing it through coupling constants (SSM/SSNME, g_V/g_A fits, detector response modeling) is analysis work that lives in separate dedicated projects and is **out of scope** for this repository since v0.5.0.
 
 ______________________________________________________________________
 
-## Package Architecture (v0.4.0 — In Progress)
+## Package Architecture (v0.5.0)
 
-### Current State (v0.3.x)
-
-Everything lives in `beta_spectrum/` — a monolithic package where theory, fitting, detector response, and analysis are intermingled. This is becoming unwieldy.
-
-### Target Architecture (v0.4.0)
+Single distributable, `beta_spectrum`:
 
 ```
-beta_spectrum/              ← theory-only core
-├── __init__.py
-├── base.py                 ← SpectrumComponent base
-├── constants.py            ← physical constants
-├── utils.py                ← T_to_W, W_to_T, momentum
-├── spectrum.py             ← BetaSpectrum, SpectrumConfig, BranchConfig
-├── components/             ← individual correction factors
-│   ├── phase_space.py
-│   ├── fermi.py
-│   ├── finite_size.py
-│   ├── screening.py
-│   ├── exchange.py
-│   ├── radiative.py
-│   └── detector_response.py   ← stays here (part of theory)
-└── nuclear_data.py         ← paceENSDF integration
-
-exp_data/                   ← experimental data handling
-├── __init__.py
-├── raw_reader.py           ← raw data loading (ROOT, CSV, ASCII)
-├── calibration.py          ← energy calibration, line fitting
-├── fitters.py              ← simple Gaussian fitters for cal sources
-├── corrections.py          ← dead-time, pile-up, background
-└── spectrum.py             ← experimental spectrum container
-
-fitter/                     ← main analysis framework
-├── __init__.py
-├── model.py                ← theoretical model + detector convolution
-├── fit_engine.py           ← multi-parameter fitting orchestration
-├── extractor.py            ← C(W), g_V, g_A extraction
-├── result.py               ← FitResult, analysis results
-# └── report.py               ← PDF report generation (deferred)
-
-output/                     ← output, visualization, CSV export
-├── __init__.py             ← BetaSpectrumAnalyzer, CSVEmitter, Visualizer
-├── core.py                 ← BetaSpectrumAnalyzer base class
-├── csv_export.py           ← CSV export functionality
-└── visualization.py        ← Plotting and visualization functionality
+beta_spectrum/               ← theory-only core package
+│   ├── __init__.py          ← public API
+│   ├── base.py              ← SpectrumComponent base
+│   ├── constants.py         ← physical constants
+│   ├── utils.py             ← T_to_W, W_to_T, momentum
+│   ├── spectrum.py          ← BetaSpectrum, SpectrumConfig, BranchConfig
+│   ├── nuclear_data.py      ← paceENSDF integration + JSON input
+│   ├── cli.py               ← bs_pnpi command-line interface
+│   ├── logging_utils.py     ← structured logging helpers
+│   ├── components/          ← individual correction factors
+│   │   ├── phase_space.py
+│   │   ├── fermi.py
+│   │   ├── finite_size.py
+│   │   ├── screening.py
+│   │   ├── exchange.py
+│   │   └── radiative.py
+│   ├── data/
+│   │   └── exchange_coeff.csv
+│   └── output/              ← visualization, CSV export subpackage
+│       ├── core.py          ← BetaSpectrumAnalyzer
+│       ├── csv_export.py    ← CSVEmitter
+│       └── visualization.py ← Visualizer, parse_spectrum_csv
 ```
 
 ### Rationale
 
-- **`beta_spectrum`** — pure theory. Given a nuclide and corrections, produces a theoretical spectrum. No fitting, no data, no detector convolution.
-- **`exp_data`** — experimental data handling. Loading raw detector outputs, calibrating energy scale, fitting calibration peaks, applying instrumental corrections.
-- **`fitter`** — the analysis workhorse. Takes theory model + detector response + experimental data → performs convolution + fitting → extracts physics parameters.
+- **Pure theory.** Given a nuclide and corrections, the package produces a theoretical spectrum. No fitting, no experimental data, no detector convolution.
+- **`beta_spectrum.output`** is a calculation-output subpackage (plotting and CSV export of computed spectra), not data analysis.
 
-### Migration Plan
+### Scope Reduction (v0.5.0)
 
-#### Phase 1: Create `exp_data` module ✅ COMPLETE
-
-- [x] `exp_data/spectrum.py` — `ExpSpectrum` dataclass: energies, counts, errors, metadata
-- [x] `exp_data/calibration.py` — `EnergyCalibrator` + `CalibrationResult`: linear/quadratic channel→energy
-- [x] `exp_data/fitters.py` — `GaussianFitter` + `PeakFitter`: peak detection and fitting
-- [x] `exp_data/corrections.py` — `DeadTimeCorrection`, `PileUpCorrection`, `BackgroundSubtractor`
-- [x] 43 tests, ruff + mypy clean
-
-#### Phase 2: Create `fitter` module ✅ COMPLETE
-
-- [x] `fitter/model.py` — `SpectrumModel`: wraps `BetaSpectrum` + `DetectorResponse`, callable model with convolution
-- [x] `fitter/fit_engine.py` — `SpectrumFitter` + `FitConfig` + `FitResult`: χ² fitting, covariance, profile likelihood, confidence intervals
-- [x] `fitter/extractor.py` — `CWExtractor` + `GVAExtractor`: C(W) and g_V/g_A extraction
-- [x] `fitter/result.py` — `AnalysisResult` + `FitSummary`: physics-specific result container
-- [x] 34 tests, ruff + mypy clean
-- [x] Detector convolution handles W < 1.0 edge cases
-
-#### Phase 3: Clean up `beta_spectrum` ✅ COMPLETE
-
-- [x] Remove `cw_extractor.py` (moved to `fitter/extractor.py`)
-- [x] Remove `fitter.py` (moved to `fitter/fit_engine.py`)
-- [x] Clean up `__init__.py` — only exports theory classes
-- [x] Add `CurveFitter` general-purpose fitter to `fitter/fit_engine.py` for backward compat
-- [x] Add `CWExtractor.fit_parametrization()` and `fit_gV_gA()` compatibility methods
-- [x] Remove old `test_cw_extractor.py` (replaced by `tests/fitter/test_extractor.py` — file was removed)
-- [x] 322 tests pass, ruff + mypy clean on all 26 source files
-
-#### Phase 4: Update package structure ✅ COMPLETE
-
-- [x] `pyproject.toml` auto-discovers all 4 packages via `setuptools.packages.find`
-- [x] `beta_spectrum.cli` keeps CLI entry point (no re-exports needed)
-- [x] All imports updated across 26 source files
-- [x] 322 tests pass
-- [x] README updated with new architecture diagram and Quick Start examples
-- [x] Version bumped to 0.4.0
+- `exp_data/` and `fitter/` packages (experimental data handling, calibration, χ² fit engine, C(W)/g_V/g_A extraction), their tests, the `bs_exp` CLI and the ROOT extra: removed from this repository (restore from history up to `67de4e3` if an analysis project needs them).
+- `detector_response.py`: removed earlier (up to `5f99566`); the deletion decision is now final.
+- `data/raw/` local copies of detector runs and Am-241 calibration files: removed (canonical copies on experiment storage).
+- `docs/gA-extraction/`, `docs/nuclear-structure/`, `docs/refs/`, `docs/corrections/07-detector-response.md`: removed; literature and nuclear-structure research is to be redone in dedicated projects.
 
 ______________________________________________________________________
 
-## Track A: ⁹⁹Tc Spectrum-Shape Analysis
+## Roadmap: Calculator Capability
 
-### A1. Research & Data Collection
+### C1. Nuclear Data & Input Flexibility
 
-- [ ] Parse additional PDF sources on ⁹⁹Tc beta decay and spectrum-shape method
-  - Kostensalo & Suhonen, Phys. Rev. C **96**, 024317 (2017) — theoretical spectrum shape predictions
-  - Suhonen, Frontier. Phys. **5**, 55 (2017) — review of g_A quenching
-  - Haaranen et al., Phys. Rev. C **93**, 034308 (2016) — effective coupling constants
-  - Additional experimental ⁹⁹Tc spectra (historical data for comparison)
-- [ ] Research methods to extract g_V from fitted C(W) shape
-  - Study the relationship between C(W) parametrization and g_V, g_A effective values
-  - Understand how CVC (Conserved Vector Current) constrains g_V
-  - Review the spectrum-shape method formalism (Behrens-Bühring + shell model)
-- [ ] Document the full analysis pipeline: experimental data → theoretical model → C(W) extraction → g_V extraction
+- [x] paceENSDF integration: `get_decay_info_from_paceENSDF()`, `decay_info_to_config()`, `create_config_from_source()`, FORBIDDENNESS_MAP, `--nuclide` CLI flag
+- [x] Measured per-branch transition endpoints from ENSDF (max endpoint across branches defines the spectrum endpoint)
+- [x] JSON custom input (`load_json_input()`, `json_to_config()`), CLI flags override file values only when explicitly given
+- [ ] Input sanity checks:
+  - `|Z_parent - Z_daughter|` must equal 1 (beta decay)
+  - For `beta_minus`: Z_daughter = Z_parent + 1; for `beta_plus`/`ec`: Z_daughter = Z_parent - 1
+  - Warn if decay_type contradicts Z values
+  - Validate `endpoint_MeV > 0` and `endpoint_MeV > level_energy_keV/1000`
+  - Cross-check transition_type against ENSDF forbiddenness
 
-### A2. Detector Response Function
+### C2. Multiple Branch Decays
 
-- [x] ~~Implement analytical detector response function for 4π semiconductor detector~~
-  - ~~Gaussian core with low-energy tail (charge collection effects)~~
-  - ~~Parameters: energy resolution σ(E), tail fraction, tail shape~~
-  - ~~Energy-dependent resolution σ(E) = a + b·√E~~
-  - ~~Implemented: Gaussian, Gaussian+tail, Tikhonov models with energy-dependent resolution and Fano factor~~
-- [x] ~~Implement convolution routines~~
-  - ~~Convolve theoretical spectrum with detector response for comparison to data~~
-  - ~~Implemented: convolve() with tabulated response support~~
-- [ ] Monte Carlo simulation (GEANT4) — out of scope for this project but noted for future
-  - Full detector geometry simulation
-  - Energy loss in source, dead layers, absorber
-  - Response matrix for unfolding (as done in Paulsen et al. with EGSNRC)
+- [x] `BranchConfig` dataclass, multi-branch `BetaSpectrum` with per-branch calculators
+- [x] Intensity normalization to 100% across branches; per-branch spectra masked beyond their endpoint
+- [x] `--intensity-cutoff` CLI filter
+- [x] Core invariant: phase space × Fermi function included unconditionally (`phase_space` and `fermi` have no toggles)
+- [x] Physical branch spectra and universal factor in CSV export; `--separate-branches` normalization mode
 
-### A3. Data Processing & Analysis
+### C3. Recoil and Weak Magnetism
 
-- [ ] Data processing pipeline for 1-hour run spectra
-  - Background subtraction
-  - Energy calibration (using known X-ray/gamma lines)
-  - Dead-time correction
-  - Pulse pile-up correction
-- [ ] Endpoint fitting and energy calibration
-  - Kurie plot analysis for endpoint determination
-  - Polynomial calibration using known spectral lines
-  - Linearity correction across energy range
-- [ ] Data quality assessment procedures
-  - Residual analysis (data vs model)
-  - χ²/ndf evaluation
-  - Consistency checks across multiple runs
-  - Statistical combination of long-exposure spectra
+Theory documented in `docs/corrections/06-recoil.md`; implementation pending.
 
-### A4. Fitter for C(W) Extraction
+- [ ] Implement recoil correction (R_N, Q terms, weak magnetism) as a `SpectrumComponent`
 
-- [x] ~~Implement fitter routine to extract experimental C(W) from data~~
-  - ~~χ² minimization: theoretical spectrum (with all corrections) vs measured spectrum~~
-  - ~~Free parameters: C(W) shape parameters, endpoint energy, normalization, background~~
-  - ~~Covariance matrix and uncertainty propagation~~
-  - ~~Implemented: `CurveFitter` with `least_squares` optimization, `FitResult` with covariance, chi2, residuals~~
-- [x] ~~Parametrize C(W) in terms of g_V (and g_A)~~
-  - ~~Fit C(W) data to theoretical parametrization~~
-  - ~~Extract g_V^eff and g_A^eff from the fit~~
-  - ~~Compare with Paulsen et al. results: g_A^eff = 0.574(36), g_V^eff = 0.376(5)~~
-  - ~~Implemented: `CWExtractor` with Kurie plot analysis, parametrized fitting, g_V/g_A extraction~~
-- [ ] Systematic uncertainty analysis
-  - Vary correction implementations within uncertainties
-  - Test sensitivity to detector response model
-  - Background model dependence
+### C4. Beta+ / Electron Capture Support
 
-______________________________________________________________________
+- [ ] Fermi function sign change for positrons
+- [ ] Radiative corrections for β⁺ (Sirlin corrections differ for e⁻ vs e⁺)
+- [ ] Screening correction for positrons (different atomic interaction)
+- [ ] Exchange correction (negligible for β⁺ but check)
 
-## Track B: Package Usability Improvements
+### C5. Neutrino Spectrum Support
 
-### B1. Nuclear Data Integration
+- [ ] Convenience analyser methods for plotting neutrino spectra (convert electron spectrum, plot both scaled)
+- [ ] Neutrino radiative correction h(W̃, W0) per Sirlin 2011; conversion f_ν(W_ν) = f_e(W̃) · [1 + (α/2π)(h − g)]
 
-- [x] ~~paceENSDF integration for automated nuclear data retrieval~~
-  - ~~`get_decay_info_from_paceENSDF(nuclide, decay_type)` — retrieves Q-value, half-life, spin/parity, branches, forbiddenness from ENSDF~~
-  - ~~`decay_info_to_config()` — converts DecayInfo to SpectrumConfig with branch support~~
-  - ~~`create_config_from_source("paceENSDF", nuclide="Tc99")` — one-liner config creation~~
-  - ~~`_parse_nuclide_symbol()` — parses "Tc99" → (Tc, 43, 99) with full Z lookup for Z=1..98~~
-  - ~~`_resolve_decay_index()` — finds correct decay index (ground state by default)~~
-  - ~~`_get_decay_mode_symbol()` — maps beta_minus→BM, beta_plus/ec→ECBP~~
-  - ~~FORBIDDENNESS_MAP — translates paceENSDF codes (0A, 1F, 2F, …) to transition_type~~
-  - ~~CLI `--nuclide Tc99` uses paceENSDF automatically~~
-  - ~~Branch intensities, log_ft, level energies, transition types all populated from ENSDF~~
-  - ~~Tested with Tc99 (2 branches, F2) and Co60 (multi-branch)~~
+### C6. Neutrino Mass Effects
 
-### B2. Input Flexibility
+Implemented: `neutrino_mass_MeV` config field, modified phase space factor √(1 − m_ν²/(W₀−W)²), `--neutrino-mass` CLI option.
 
-- [x] ~~Generalize parameter inputs via custom input file (YAML/JSON)~~
-  - ~~Declare isotope, transition type, detector parameters~~
-  - ~~Toggle corrections on/off per-component~~
-  - ~~Reproducible calculation configurations~~
-  - ~~Implemented: `load_json_input()`, `json_to_config()`, `DEFAULT_JSON_SCHEMA` with full detector param support, CLI `--input` flag~~
-- [x] ~~Retrieve parameters directly from ENSDF database~~
-  - ~~Auto-populate W0, Z, A, transition type from isotope name~~
-  - ~~Reduce manual configuration errors~~
-  - ~~Implemented: `get_decay_info_from_paceENSDF()` returns DecayInfo with Z_parent, Z_daughter, A_number, endpoint_MeV, transition_type, forbiddenness_code, half_life, branches; CLI `--nuclide` uses this automatically~~
+- [ ] Document which corrections need m_ν as an additional input parameter
+- [ ] Study impact on radiative corrections (mass singularities)
 
-### B3. Multiple Branch Decays
+### C7. Automated Report Generation
 
-- [x] ~~Add support for multiple branch decays~~
-  - ~~Calculate total spectral shape weighted by branch intensities~~
-  - ~~Handle decays with significant branching to excited states~~
-  - ~~Sum contributions from all branches with proper weighting~~
-  - ~~Implemented: `BranchConfig` dataclass, multi-branch `BetaSpectrum` with per-branch calculators~~
-  - ~~Branch data from paceENSDF: intensity, endpoint (Q-value - excitation energy), transition type~~
-  - ~~Intensity normalization to 100% across all branches~~
-  - ~~CSV export: total spectrum, per-branch spectra, per-branch components (PhaseSpace, Radiative), universal components (Fermi, Screening, etc.)~~
-  - ~~Plot: total spectrum + individual branches with legend (transition type, endpoint energy)~~
-  - ~~Debug plot: vertical layout with per-branch panels, universal components panel~~
-  - ~~Energy grid extends to max branch endpoint, per-branch spectra masked beyond their endpoint~~
-  - ~~Single-branch = trivial case of multi-branch (no separate mode)~~
-  - ~~Per-branch transition types extracted from ENSDF (not defaulting to 'A')~~
-  - ~~`--intensity-cutoff` CLI option to filter out negligible branches~~
-  - ~~Branch spectra plotted as raw (un-normalized) contributions so total = sum of branches~~
-
-### B4. Beta+ Decay Support
-
-- [ ] Modify corrections to support β⁺ decays
-  - Fermi function sign change for positrons
-  - Radiative corrections for β⁺ (Sirlin corrections differ for e⁻ vs e⁺)
-  - Screening correction for positrons (different atomic interaction)
-  - Exchange correction (negligible for β⁺ but check)
-
-### B5. Neutrino Spectrum Support
-
-- [ ] Add convenience analyser methods for plotting neutrino spectra
-  - Convert electron spectrum to neutrino spectrum
-  - Plot both on same figure with proper scaling
-- [ ] Adopt corrections that differ for neutrino spectrum
-  - Implement neutrino radiative correction h(W̃, W0) per Sirlin 2011
-  - Conversion formula: f_ν(W_ν) = f_e(W̃) · [1 + (α/2π)(h − g)]
-  - Neutrino-specific corrections (no Fermi function, different radiative terms)
-
-### B6. Neutrino Mass Effects
-
-- [ ] Implement option to produce spectra with non-zero neutrino mass
-  - Modify phase space: (W0 − W)² → (W0 − W)² · √(1 − m_ν²/(W0−W)²)
-  - Add kink feature near endpoint for spectral distortion
-  - Useful for neutrino mass search applications
-- [ ] Research and document implications of significant neutrino mass
-  - Which corrections need m_ν as input parameter?
-  - Impact on radiative corrections (mass singularities)
-  - Impact on endpoint fitting procedures
-
-### B7. Automated Report Generation
-
-- [ ] Implement extensive procedurally generated PDF report
-  - Comprehensive details of the calculation
+- [ ] Procedurally generated PDF report for a calculation:
+  - Parameter table (nuclide, decay data, enabled corrections, version, commit)
   - Component plots (phase space, Fermi, each correction factor)
   - Total spectrum with all corrections applied
-  - Residual analysis, fit quality metrics
-  - Parameter table with uncertainties
   - References to theoretical sources
-- [ ] Extend report generation to fitters
-  - Fit results with covariance matrix
-  - Parameter correlations
-  - Confidence intervals
-  - Goodness-of-fit statistics
 
-### B8. CLI & Output Improvements
+### C8. CLI & Output Improvements
 
-- [x] Implement structured logging system
-  - Log to stdout with configurable verbosity levels (-v=INFO, -vv=DEBUG)
-  - Visualize workflow steps: data loading, parameter resolution, calculation progress
-  - Log all used parameters: Z_parent, Z_daughter, A_number, endpoint, transition type, enabled corrections
-  - Optional log file output with timestamped filename
-  - Integrate Python `logging` module with custom formatter
-- [x] Add metadata header to CSV output files
-  - Date and time of calculation (ISO 8601 format)
-  - Nuclide information: parent/daughter symbol, Z, A
-  - Decay/transition info: endpoint energy, transition type, forbiddenness
-  - Calculation parameters: energy step, enabled corrections, detector settings
-  - Software version and git commit hash (7 chars)
-  - Format: YAML-style comments at top of CSV file
-- [ ] Optimize CLI argument design and add sanity checks
-  - Auto-deduce `transition_type` from `decay_type` and nuclear data — remove as explicit CLI parameter
-  - Auto-deduce `decay_type` from Z_parent vs Z_daughter difference — remove as explicit CLI parameter
-  - Add sanity checks for custom input:
-    - `|Z_parent - Z_daughter|` must equal 1 (beta decay)
-    - For `beta_minus`: Z_daughter must be Z_parent + 1
-    - For `beta_plus`/`ec`: Z_daughter must be Z_parent - 1
-    - Warn if decay_type contradicts Z values (e.g., beta+ with Z_daughter > Z_parent)
-    - Validate endpoint_MeV > 0 and endpoint_MeV > level_energy_keV/1000
-    - Cross-check transition_type against ENSDF forbiddenness (if paceENSDF source)
-- [x] Add `--dry-run` option to validate input and display resolved parameters without calculation
-- [x] Add `--version` flag
-- [x] Add `-q/--quiet` flag to suppress terminal output
-- [x] Add `./output/` directory for debug verification artifacts (gitignored except `.gitkeep`)
-- [x] Add debug verification step to development workflow: `bs_pnpi -vv` with parameter consistency checks across all components, documented in DEVELOPMENT.md §4.5
-- [x] Add plot output with ID tracking (commit hash + UTC timestamp) and two display modes:
-  - Normal mode (default): single spectrum plot with nuclear data header
-  - Debug mode (-vv): 4-panel view with all correction factors and deviations
+- [ ] Auto-deduce `transition_type` from nuclear data (remove as explicit parameter where possible)
+- [ ] Auto-deduce `decay_type` from Z_parent vs Z_daughter difference
+- [x] Structured logging (`-v`, `-vv`, `-q`, `--log-file`)
+- [x] CSV metadata headers (ISO timestamp, nuclide, decay, calculation params, version, commit hash)
+- [x] `--dry-run`, `--version`, single-letter aliases for every long option
+- [x] Plots built from the exported CSV; commit hash + UTC timestamp stamps; normal/debug display modes
 
 ______________________________________________________________________
 
@@ -300,58 +121,46 @@ ______________________________________________________________________
 
 ### Additional Corrections
 
-Not implementing at this time. Additional corrections (weak magnetism, recoil kinematics, etc.) require a complex implementation with project-wide consequences. This is a far-future plan, irrelevant to the current ⁹⁹Tc analysis goal where the shape factor C(W) absorbs all unmodelled energy-dependent effects.
+Additional corrections beyond the implemented set require complex implementation with project-wide consequences; deferred as far-future work.
 
-### Nuclear Shape Factors C(Z,W)
+### Theoretical Shape Factors C(Z,W)
 
-Not implementing theoretical shape factor calculations. The entire point of this project is to **extract** C(W) from experimental data, not to compute it from nuclear theory (shell model, QRPA, etc.). Theoretical shape factor calculations require:
-
-- Nuclear structure models (NUSHELLX, shell model codes)
-- One-body transition densities (OBTDs)
-- Coulomb displacement energy calculations
-- Multiple effective interactions
-
-These are handled by the nuclear physics community. Our role is to provide the analysis framework that extracts C(W) from data, then compare the extracted shape to theoretical predictions.
+Never computed here. Shell-model / QRPA transition densities and effective interactions belong to nuclear-structure codes; this repository treats C(W) as an input parameter of the spectrum model. The dedicated research into nuclear shape factors will be redone in separate projects.
 
 ______________________________________________________________________
 
 ## Current Status
 
-**Version:** 0.3.1 → **In progress: 0.4.0 — Package restructuring**
+**Version:** v0.5.0 — theory-only scope reduction. Calculation logic unchanged from v0.4.1.
 
-### Completed (v0.3.x)
+### Completed (v0.3.x — theory core)
 
-- All theoretical corrections: phase space, Fermi function, finite size, screening, exchange, radiative (with delta_cut resummation)
-- Detector response module with analytical models (Gaussian, Gaussian+tail, Tikhonov), convolution API
-- χ² curve fitting framework (CurveFitter) with confidence intervals, profile likelihood, correlation analysis
-- C(W) shape factor extraction pipeline (CWExtractor) with Kurie plot analysis, parametrized fitting, g_V/g_A extraction
+- All implemented theoretical corrections: phase space, Fermi function, finite size, screening, exchange, radiative (with delta_cut resummation)
 - Multi-branch decay support (unified single/multi-branch)
-- CLI interface (`bs_pnpi`) with paceENSDF integration (`--nuclide`), structured logging (-v/-vv/-q), --dry-run, --version, --log-file
+- CLI (`bs_pnpi`) with paceENSDF integration, structured logging, `--dry-run`, `--version`, `--log-file`
 - CSV metadata headers with branch info
-- JSON input with full detector param support
-- paceENSDF integration: `get_decay_info_from_paceENSDF()` returns DecayInfo with Z, A, endpoint, transition type, half-life, spin/parity, branches
-- Comprehensive test suite (262 tests, including 6 paceENSDF integration tests)
+- paceENSDF `DecayInfo` retrieval (Z, A, endpoint, transition type, half-life, spin/parity, branches)
 - Notebook quality control with nbmake and auto-save plot hooks
-- Plot output: unified multi-branch approach, externalized headers, compact legends, consistent titles
+- Plot output: unified multi-branch approach, externalized headers, consistent titles
 
-### Completed (v0.4.0)
+### Completed (v0.4.0 — package split, later obsolete)
 
-- **Track B: Package restructuring** — separating into `beta_spectrum` (theory), `exp_data` (experimental data), `fitter` (analysis framework)
-  - [x] Phase 1: Create `exp_data` module (raw_data.py, calibration.py, fitters.py, corrections.py, spectrum.py) — 43 tests
-  - [x] Phase 2: Create `fitter` module (model.py, fit_engine.py, extractor.py, result.py) — 34 tests
-  - [x] Phase 3: Clean up `beta_spectrum` (remove cw_extractor.py, fitter.py, **init**.py cleanup) — 322 tests
-  - [x] Phase 4: Package auto-discovered by setuptools, all imports updated, black formatting applied
-- **v0.4.0 is COMPLETE**
-  - 322 tests pass (262 original + 60 new)
-  - ruff + mypy clean on all 26 source files
-  - black formatting applied to all new code
-  - `pyproject.toml` auto-discovers all 4 packages: `beta_spectrum`, `beta_spectrum.output`, `exp_data`, `fitter`
-  - README updated with new architecture diagram and Quick Start examples
-  - Version bumped to 0.4.0
+- Interim split into `beta_spectrum` / `exp_data` / `fitter` packages
+- Superseded by the v0.5.0 scope reduction: the analysis packages were removed entirely
+
+### Completed (v0.4.1 — normalization & CLI)
+
+- 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: `--separate-branches` rename, single-letter aliases, dry-run cleanup, stale detector_* JSON schema keys removed
+
+### Completed (v0.5.0 — scope reduction)
+
+- Removal of `exp_data/`, `fitter/`, raw detector data, analysis notebooks/scripts
+- Documentation rewrite to theory-only scope
+- 280 tests pass (theory-only suite), ruff + mypy clean
 
 ### Next Immediate Step
 
-- **Track A: ⁹⁹Tc spectrum-shape analysis** — start real data processing pipeline
-  - Data processing pipeline for experimental runs
-  - Endpoint fitting and energy calibration
-  - Systematic uncertainty analysis
+- C1 input sanity checks (decay/Z consistency, endpoint validation)
+- C3 recoil & weak magnetism correction implementation

beta_spectrum/__init__.py

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

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
 }