Merge dev/cli-logging-refinement to main

- Fix paceENSDF Q-value extraction (was reading half-life as Q-value)
- Fix RadiativeCorrection Z/A parameters
- Add debug verification workflow (bs_pnpi -vv)
- Add plot output with ID tracking and two display modes
- Fix pyright/mypy type errors across codebase
- Add ./output/ directory for test artifacts
- Update docs and development workflow

.gitignore

@@ -81,6 +81,8 @@ ipython_config.py
 # data/
 outputs/
 results/
+output/*
+!output/.gitkeep
 *.npy
 *.npz
 *.h5

CONVENTIONS.md

@@ -66,7 +66,35 @@ from beta_spectrum.constants import ME_MEV, ME_KEV, ALPHA, MP_MEV
 - Never use manual arithmetic for unit conversion — use `T_to_W()`, `W_to_T()`
 - Never mix MeV and m_e units in the same calculation
 
-## 3. SpectrumComponent Design
+## 3. Energy Unit Convention
+
+### MeV Throughout
+
+All energy parameters in the public API use **MeV** as the standard unit:
+
+- `endpoint_MeV` — always in MeV
+- `e_step_MeV` — always in MeV
+- `detector_sigma_a_keV` — in keV (detector-specific, exception noted)
+- `detector_tau_keV` — in keV (detector-specific, exception noted)
+
+Detector parameters use keV because detector resolutions are typically specified in keV in the literature. All other energies use MeV.
+
+### CLI Unit Convention
+
+All CLI energy parameters use MeV unless explicitly marked otherwise:
+
+- `--e-step` — MeV (default: 0.001)
+- `--sigma` — keV (detector resolution, default: 1.0)
+- `--tau` — keV (tail decay constant, default: 5.0)
+
+### Consistency Rule
+
+When reading or writing values, always check the unit suffix in the parameter name:
+- `_MeV` suffix → value is in MeV
+- `_keV` suffix → value is in keV
+- No suffix on energy parameter → value is in m_e units (internal)
+
+## 4. SpectrumComponent Design
 
 All spectral components (Fermi function, screening, radiative correction, etc.) follow a common interface:
 
@@ -87,7 +115,42 @@ class MyComponent(SpectrumComponent):
         return ...
 ```
 
-## 4. Physics-Specific Testing Guidelines
+## 5. Logging Convention
+
+### Verbosity Levels
+
+- **No verbosity (default)**: Only errors and warnings
+- **`-v` (INFO)**: Shows source selection, component initialization, workflow steps, and output locations
+- **`-vv` (DEBUG)**: Shows all INFO messages plus detailed internals of each spectrum component (parameters, input/output ranges, evaluation counts)
+
+### Component Logging
+
+Each `SpectrumComponent` logs at two levels:
+
+- **INFO** at initialization: Component name and key parameters
+- **DEBUG** at evaluation: Number of energy points, input/output value ranges
+
+This allows debugging without cluttering normal output.
+
+## 6. Nuclear Data Notation
+
+### Element and Nuclide Notation
+
+- Element symbol: `Tc`, `Ru`, `U`, etc. (from atomic number)
+- Nuclide notation: `Tc99`, `Ru99`, `U238` (element symbol + mass number)
+- Decay notation: `Tc99 -> Ru99` (parent -> daughter)
+
+### CSV Header Convention
+
+CSV metadata headers use element notation for readability:
+
+```
+# nuclide: Tc99 -> Ru99 (Z=43->44)
+```
+
+This provides both the human-readable notation and the machine-parseable Z values.
+
+## 7. Physics-Specific Testing Guidelines
 
 When writing tests for physics components, verify the following:
 

DEVELOPMENT.md

@@ -111,6 +111,103 @@ Notebooks in `notebooks/` are executed as pytest tests using `nbmake`. This ensu
    pytest notebooks/ --nbmake  # or just `pytest` (includes notebooks/)
    ```
 
+______________________________________________________________________
+
+### 4.5. Debug Verification (Final Step of Each Dev Iteration)
+
+After all quality gates pass, run a **debug-level verification** as the final step of every development iteration. This catches configuration mismatches, missing parameter propagation, and silent logic errors before merging.
+
+**Purpose:** Verify that all parameters (endpoint energy, energy grid, Z, A, physical constants) are consistent throughout the pipeline, and that every spectrum component emits sufficient debug logging for explicit comparison.
+
+**Procedure:**
+
+1. **Run `bs_pnpi` with `-vv` (DEBUG logging) and output to `./output/`:**
+   ```bash
+   bs_pnpi --nuclide Tc99 -vv --output ./output/verification.csv
+   ```
+
+2. **Capture full debug output to a log file for analysis:**
+   ```bash
+   bs_pnpi --nuclide Tc99 -vv --output ./output/verification.csv --log-file ./output/debug.log 2>&1 | tee ./output/debug_stdout.log
+   ```
+
+3. **Verify consistency of key parameters in the debug output:**
+
+   | Parameter | What to check | Where to find it |
+   |---|---|---|
+   | `endpoint_MeV` | Matches expected Q-value | CLI dry-run, config creation, energy grid |
+   | `Z_parent` / `Z_daughter` | Consistent across all components | Component init logs |
+   | `A_number` | Consistent across all components | Component init logs |
+   | `e_step_MeV` | Matches grid size × (endpoint - e_step) | Energy grid creation |
+   | `W0` (endpoint in natural units) | `W0 = 1.0 + endpoint_MeV / m_e` | PhaseSpace init |
+   | Physical constants | `m_e = 0.51099895` MeV, `α = 1/137.036` | Constants module, Fermi function |
+   | Energy grid bounds | `e_step` to `endpoint - e_step` | `get_energy_grid()` output |
+
+4. **Verify each component's debug logging:**
+
+   Every enabled component must log its initialization parameters and a representative evaluation. Check the debug output for entries like:
+
+   ```
+   DEBUG ... Initializing PhaseSpace: W0=1.57664, transition_type=A, m_e=1.0
+   DEBUG ... Evaluating PhaseSpace at 294 energy points
+   DEBUG ... PhaseSpace range: [min=..., max=...]
+   ```
+
+   Required log entries per component:
+
+   | Component | Must log |
+   |---|---|
+   | `PhaseSpace` | `W0`, `transition_type`, `m_e`, `m_nu`, energy range |
+   | `FermiFunction` | `Z`, `A`, `α`, Coulomb approximation used |
+   | `FiniteSizeL0` | `Z`, `A`, nuclear radius `R`, L0 approximation |
+   | `ChargeDistributionU` | `Z`, `A`, charge radius, U correction range |
+   | `ScreeningCorrection` | `Z_parent`, screening model, S(Z,W) range |
+   | `ExchangeCorrection` | `Z_parent`, Hayen coefficients used, X(Z,W) range |
+   | `RadiativeCorrection` | `W0`, `δ_r` formula, resummation flag, `delta_cut` value |
+
+5. **Report any missing logging:**
+
+   If a component's debug output does not allow explicit verification of its parameters or computed range, **add logging to that component** before the iteration is considered complete. The goal is that a developer (or an LLM) reading the debug log can reconstruct every parameter and verify correctness without reading source code.
+
+6. **Save debug artifacts to `./output/`:**
+
+   After verification, ensure all artifacts are in `./output/`:
+   - `debug_stdout.log` — full stdout+stderr from `bs_pnpi -vv`
+   - `debug.log` — log file output (if `--log-file` was used)
+   - `verification.csv` — the spectrum CSV with metadata header
+   - `verification.png` — analysis plot (if `--plot` was used)
+
+   These files serve as a **verifiable record** of the last working state and can be used for regression comparison.
+
+**Example verification command:**
+
+```bash
+# Full debug run with all outputs
+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
+
+# Quick parameter extraction from debug log
+grep -E '(Initializing|Evaluating|W0=|Z=|endpoint|e_step|PhaseSpace|Fermi|FiniteSize|Screening|Exchange|Radiative)' ./output/debug_stdout.log
+```
+
+**Success criteria:**
+
+- [ ] All parameters (Z, A, endpoint, e_step, W0) appear consistently in debug output
+- [ ] Each enabled component logs its initialization parameters
+- [ ] Each enabled component logs its evaluation energy range
+- [ ] Physical constants match expected values
+- [ ] No silent parameter coercion or unit conversion errors detected
+- [ ] All artifacts saved to `./output/`
+
+**If any check fails:**
+
+1. Add missing debug logging to the component
+2. Fix any parameter mismatch
+3. Re-run the verification until all checks pass
+4. Only then consider the iteration complete
+
 ## 5. Agent Instructions
 
 When working on this project, agents must:
@@ -123,6 +220,7 @@ When working on this project, agents must:
 6. **Work on a dev branch** — create `dev/<description>` and switch to it
 7. **Run quality gates** after making changes and report results
 8. **Update docs** when changing public API or adding features
+9. **Run debug verification** (`bs_pnpi -vv`) as the final step of each iteration, saving all artifacts to `./output/` and verifying parameter consistency across all components
 
 ## 6. Release Checklist
 
@@ -130,6 +228,8 @@ When working on this project, agents must:
 - [ ] Code is formatted (`black --check`)
 - [ ] No lint errors (`ruff check`)
 - [ ] Type checks pass (`mypy`)
+- [ ] Debug verification passes (`bs_pnpi -vv`) — all parameters consistent, all components logged
+- [ ] Debug artifacts saved to `./output/` (debug_stdout.log, verification.csv, etc.)
 - [ ] Update version in `pyproject.toml`
 - [ ] Update `Development Status` section in `README.md`
 - [ ] Update `TODO.md` marking corresponding items as complete.

README.md

@@ -8,7 +8,11 @@ A modular Python package for calculating high-precision **beta decay energy spec
 - Selective inclusion of higher-order corrections via `SpectrumConfig` toggles
 - Clean, composable component architecture
 - Built-in visualization and CSV export via `BetaSpectrumAnalyzer`
+  - 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
 - Comprehensive scientific reference in [`docs/`](docs/)
+- Debug verification workflow: `bs_pnpi -vv` validates all parameters end-to-end
 
 ______________________________________________________________________
 
@@ -74,6 +78,7 @@ beta_shape_pnpi/
 └── 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)
 ```
 
 ______________________________________________________________________
@@ -104,8 +109,11 @@ Defined in `beta_spectrum/spectrum.py`. Orchestrates the calculation:
 Also in `spectrum.py`. Provides debugging and visualization tools:
 
 - `total_spectrum(normalize=True/False)` — computed spectrum array
-- `plot_analysis(save_path=None)` — 4-panel matplotlib figure (spectrum, components, cumulative effect, deviations from unity)
-- `export_to_csv(filename)` — export via pandas DataFrame
+- `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
 
 ______________________________________________________________________
@@ -128,8 +136,9 @@ W, E_MeV = spectrum.get_energy_grid(config)
 values = spectrum(W)
 
 analyzer = BetaSpectrumAnalyzer(spectrum, config)
-analyzer.plot_analysis()          # Show/ save 4-panel figure
-analyzer.export_to_csv("output.csv")   # Export data to CSV
+analyzer.plot_analysis("analysis.png")          # Normal mode (single spectrum plot)
+analyzer.plot_analysis("debug.png", show_components=True)  # Debug mode (4-panel)
+analyzer.export_to_csv("spectrum.csv")   # Export data to CSV
 ```
 
 ______________________________________________________________________

TODO.md

@@ -178,6 +178,11 @@ ______________________________________________________________________
 - [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
 
 ______________________________________________________________________
 
@@ -205,6 +210,6 @@ ______________________________________________________________________
 **Version:** 0.3.0\
 **Implemented:** Phase space, Fermi function, finite size, screening, exchange, radiative corrections (with delta_cut resummation). Detector response module with analytical models (Gaussian, Gaussian+tail, Tikhonov), convolution API, declarative config integration. χ² curve fitting framework (CurveFitter) with confidence intervals, profile likelihood, and correlation analysis. C(W) shape factor extraction pipeline (CWExtractor) with Kurie plot analysis, parametrized fitting, and g_V/g_A extraction. CLI interface (`bs_pnpi`) with paceENSDF integration, structured logging (-v/-vv/-q), --dry-run, --version, --log-file. CSV metadata headers. JSON input with full detector param support. Comprehensive test suite (195 tests). Notebook quality control with nbmake and auto-save plot hooks.
 
-**Completed:** A2 — detector response function and convolution routines. A4 — fitter routine and C(W) extraction pipeline. B2 — input flexibility (JSON input, ENSDF auto-population). B8 — CLI & output improvements (logging, CSV headers, dry-run, version).
+**Completed:** A2 — detector response function and convolution routines. A4 — fitter routine and C(W) extraction pipeline. B2 — input flexibility (JSON input, ENSDF auto-population). B8 — CLI & output improvements (logging, CSV headers, dry-run, version). B8.2 — CLI refinement: removed --mode, --transition-type, --decay-index flags; clarified units (MeV everywhere); enhanced logging (INFO shows components, DEBUG shows internals); CSV headers use element notation (Tc99 -> Ru99). B8.3 — debug verification framework: `./output/` directory for artifacts, development workflow step 4.5 with `-vv` parameter consistency checks across all components. B8.4 — plot output with ID tracking (commit hash + UTC timestamp) and two display modes (normal spectrum plot with nuclear data header; debug 4-panel view with all corrections).
 
 **Next immediate step:** A3 — data processing pipeline for experimental spectra (background subtraction, energy calibration, dead-time correction, pulse pile-up correction).

beta_spectrum/base.py

@@ -23,6 +23,14 @@ class SpectrumComponent(ABC):
 
     def __init__(self, logger: Optional["logging.Logger"] = None) -> None:
         self._logger = logger
+        self._log_component_info()
+
+    def _log_component_info(self) -> None:
+        """Log component identification at INFO level."""
+        if self._logger is None:
+            return
+        class_name = self.__class__.__name__
+        self._logger.info("Component initialized: %s", class_name)
 
     @abstractmethod
     def __call__(self, W: ArrayLike) -> ArrayLike:

beta_spectrum/cli.py

@@ -9,13 +9,13 @@ Usage examples:
     bs_pnpi --input custom_input.json --output spectrum.csv
 
     # Full analysis with plotting
-    bs_pnpi --nuclide Co60 --mode BM --plot analysis.png --output spectrum.csv
+    bs_pnpi --nuclide Co60 --plot analysis.png --output spectrum.csv
 
     # With detector response
-    bs_pnpi --nuclide Tc99 --detector sigma=1.0 --output convolved.csv
+    bs_pnpi --nuclide Tc99 --detector gaussian --output convolved.csv
 
     # Verbose with log file
-    bs_pnpi --nuclide Tc99 --output spectrum.csv -vv --log-file calc.log
+    bs_pnpi --nuclide Tc99 --output spectrum.csv -vv --log-file /tmp/calc.log
 
     # Dry run to inspect resolved configuration
     bs_pnpi --nuclide Tc99 --dry-run
@@ -28,6 +28,7 @@ import sys
 from typing import Any, Dict
 
 from beta_spectrum.logging_utils import LoggingConfig, setup_logging
+from beta_spectrum.spectrum import SpectrumConfig
 
 
 def _build_parser() -> argparse.ArgumentParser:
@@ -42,8 +43,8 @@ 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 sigma=1.0 --output convolved.csv\n"
-            "  bs_pnpi --nuclide Tc99 --output spectrum.csv -vv --log-file calc.log\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"
         ),
     )
@@ -71,31 +72,11 @@ def _build_parser() -> argparse.ArgumentParser:
     )
 
     # Common parameters
-    parser.add_argument(
-        "--mode",
-        type=str,
-        default="beta_minus",
-        choices=["beta_minus", "beta_plus", "ec"],
-        help="Decay mode for paceENSDF source (default: beta_minus).",
-    )
-    parser.add_argument(
-        "--decay-index",
-        type=int,
-        default=None,
-        help="Decay index for paceENSDF (0 = ground state, default: auto-select).",
-    )
-    parser.add_argument(
-        "--transition-type",
-        type=str,
-        default=None,
-        choices=["A", "F1", "F1U", "F2", "F2U", "F3", "F3U", "F4"],
-        help="Override transition type (forbiddenness classification).",
-    )
     parser.add_argument(
         "--e-step",
         type=float,
         default=0.001,
-        help="Energy step size in MeV (default: 0.001).",
+        help="Energy step size in MeV (default: 0.001). All energies are in MeV throughout the toolkit.",
     )
 
     # Output
@@ -122,28 +103,29 @@ def _build_parser() -> argparse.ArgumentParser:
         default=None,
         metavar="MODEL",
         help="Enable detector response convolution. "
-        "Specify model: 'gaussian', 'gaussian_tail', 'tikhonov', 'tabulated'.",
+        "Specify model: 'gaussian', 'gaussian_tail', 'tikhonov', 'tabulated'. "
+        "All detector parameters are in keV.",
     )
     detector_group.add_argument(
         "--sigma",
         type=float,
         default=1.0,
         metavar="KEV",
-        help="Resolution parameter sigma in keV (default: 1.0).",
+        help="Resolution parameter sigma in keV (default: 1.0). Energy-dependent part: sigma(E) = sigma_a + sigma_b * sqrt(E).",
     )
     detector_group.add_argument(
         "--tau",
         type=float,
         default=5.0,
         metavar="KEV",
-        help="Tail decay constant tau in keV (default: 5.0).",
+        help="Low-energy tail decay constant tau in keV (default: 5.0). Only used with 'gaussian_tail' model.",
     )
     detector_group.add_argument(
         "--tail-fraction",
         type=float,
         default=0.0,
         metavar="FRAC",
-        help="Low-energy tail fraction (0.0–1.0, default: 0.0).",
+        help="Low-energy tail fraction (0.0–1.0, default: 0.0). Only used with 'gaussian_tail' model.",
     )
 
     # Logging and output control
@@ -153,7 +135,8 @@ def _build_parser() -> argparse.ArgumentParser:
         "--verbose",
         action="count",
         default=0,
-        help="Increase verbosity: -v=INFO, -vv=DEBUG (repeatable).",
+        help="Increase verbosity: -v=INFO (shows source, components, and workflow), "
+        "-vv=DEBUG (detailed internals of all spectrum components). Repeatable.",
     )
     log_group.add_argument(
         "-q",
@@ -164,8 +147,9 @@ def _build_parser() -> argparse.ArgumentParser:
     log_group.add_argument(
         "--log-file",
         type=str,
-        metavar="FILE",
-        help="Append log messages to the specified file.",
+        metavar="PATH",
+        help="Write log messages to the specified file path (absolute or relative). "
+        "Log file is always opened in append mode at DEBUG level.",
     )
     log_group.add_argument(
         "--dry-run",
@@ -227,13 +211,12 @@ def _run(args: argparse.Namespace) -> None:
     detector_kwargs = _parse_detector_args(args)
 
     source_type: str
+    config: SpectrumConfig
     if args.nuclide:
-        logger.info("Source: paceENSDF nuclide=%s mode=%s", args.nuclide, args.mode)
+        logger.info("Source: paceENSDF nuclide=%s", args.nuclide)
         config = create_config_from_source(
             "paceENSDF",
             nuclide=args.nuclide,
-            decay_type=args.mode,
-            decay_index=args.decay_index,
             e_step_MeV=args.e_step,
             use_detector_response=detector_kwargs.get("use_detector_response", False),
             detector_model=detector_kwargs.get("detector_model", "gaussian"),
@@ -245,19 +228,12 @@ def _run(args: argparse.Namespace) -> None:
     elif args.input:
         logger.info("Source: JSON file=%s", args.input)
         data = load_json_input(args.input)
-        if args.transition_type:
-            data["transition_type"] = args.transition_type
         config = json_to_config(data)
         source_type = "json"
     else:
         logger.error("No source specified (use --nuclide or --input)")
         sys.exit(1)
 
-    # Override transition type if explicitly provided
-    if args.transition_type:
-        config.transition_type = args.transition_type
-        logger.info("Override transition_type=%s", args.transition_type)
-
     # Dry run: display resolved config and exit
     if args.dry_run:
         _print_dry_run_output(config, source_type)
@@ -266,15 +242,12 @@ def _run(args: argparse.Namespace) -> None:
     # Create and calculate spectrum
     start_time = time.time()
     logger.info(
-        "Creating BetaSpectrum for %s -> Z=%d,%d, A=%d, endpoint=%.3f MeV",
-        _element_symbol(config.Z_parent),
-        config.Z_parent,
-        config.Z_daughter,
-        config.A_number,
+        "Creating BetaSpectrum for %s, endpoint=%.3f MeV",
+        _decay_notation(config.Z_parent, config.Z_daughter, config.A_number),
         config.endpoint_MeV,
     )
 
-    spectrum = BetaSpectrum.from_config(config)
+    spectrum = BetaSpectrum.from_config(config, logger=logger)
     _W, _kinetic_MeV = spectrum.get_energy_grid(config)
     _ = spectrum(_W)
 
@@ -285,14 +258,17 @@ def _run(args: argparse.Namespace) -> None:
 
     # Export to CSV
     if args.output:
-        analyzer = BetaSpectrumAnalyzer(spectrum, config)
+        analyzer = BetaSpectrumAnalyzer(spectrum, config, logger=logger)
         analyzer.export_to_csv(args.output, source_type=source_type)
         logger.info("Spectrum exported to %s", args.output)
 
     # Generate plot
     if args.plot:
-        analyzer = BetaSpectrumAnalyzer(spectrum, config)
-        analyzer.plot_analysis(save_path=args.plot)
+        analyzer = BetaSpectrumAnalyzer(spectrum, config, logger=logger)
+        analyzer.plot_analysis(
+            save_path=args.plot,
+            show_components=(args.verbose >= 2),
+        )
         logger.info("Analysis plot saved to %s", args.plot)
 
     logger.info("Done.")
@@ -322,7 +298,9 @@ def _print_dry_run_output(config, source_type: str) -> None:
     print(f"Package:    beta-spectrum v{__import__('beta_spectrum').__version__}")
     print(f"Git commit: {get_git_short_hash()}")
     print(f"Source:     {source_type}")
-    print(f"Nuclide:    {config.Z_parent} -> {config.Z_daughter}, A={config.A_number}")
+    print(
+        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")
@@ -445,6 +423,16 @@ def _element_symbol(Z: int) -> str:
     return symbols[Z] if Z < len(symbols) else f"Z{Z}"
 
 
+def _nuclide_symbol(Z: int, A: int) -> str:
+    """Get nuclide notation like 'Tc99' from atomic and mass numbers."""
+    return f"{_element_symbol(Z)}{A}"
+
+
+def _decay_notation(Z_parent: int, Z_daughter: int, A: int) -> str:
+    """Get decay notation like 'Tc99 -> Ru99' from atomic and mass numbers."""
+    return f"{_nuclide_symbol(Z_parent, A)} -> {_nuclide_symbol(Z_daughter, A)}"
+
+
 def main() -> None:
     """Main entry point for the CLI."""
     parser = _build_parser()

beta_spectrum/components/exchange.py

@@ -1,5 +1,7 @@
 from pathlib import Path
-from typing import Dict
+from typing import Dict, Optional
+
+import logging
 
 import numpy as np
 
@@ -11,7 +13,14 @@ class ExchangeCorrection(SpectrumComponent):
     Atomic xchange correction using Hayen 2018 Table X coefficients (App. G)
     """
 
-    def __init__(self, Z: int, filename: str | None = None, eps: float = 1e-6):
+    def __init__(
+        self,
+        Z: int,
+        filename: Optional[str] = None,
+        eps: float = 1e-6,
+        logger: Optional[logging.Logger] = None,
+    ):
+        super().__init__(logger=logger)
         self.Z = Z
         self.eps = eps
 
@@ -20,6 +29,13 @@ class ExchangeCorrection(SpectrumComponent):
 
         self.coeffs = self._load_coeffs(filename)
 
+        if self._logger:
+            self._logger.debug(
+                "ExchangeCorrection: Z=%d, coefficients loaded from %s",
+                Z,
+                filename,
+            )
+
     def _load_coeffs(self, filename: str) -> Dict[str, float]:
         data = np.genfromtxt(filename, delimiter=",", names=True)
 
@@ -42,6 +58,14 @@ class ExchangeCorrection(SpectrumComponent):
         return Path(__file__).resolve().parent.parent / "data" / "exchange_coeff.csv"
 
     def __call__(self, W: np.ndarray) -> np.ndarray:
+        if self._logger:
+            self._logger.debug(
+                "ExchangeCorrection: evaluating at %d points, W range=[%.4f, %.4f]",
+                len(W),
+                W.min(),
+                W.max(),
+            )
+
         W = np.asarray(W)
 
         # Define physical cutoff
@@ -71,4 +95,13 @@ class ExchangeCorrection(SpectrumComponent):
 
         X_safe = 1.0 + term1 + term2 + term3 + term4
 
-        return np.asarray(X_safe, dtype=np.float64)
+        result = np.asarray(X_safe, dtype=np.float64)
+
+        if self._logger:
+            self._logger.debug(
+                "ExchangeCorrection: output range=[%.6e, %.6e]",
+                result.min(),
+                result.max(),
+            )
+
+        return result

beta_spectrum/components/fermi.py

@@ -1,4 +1,5 @@
 # components/fermi.py
+import logging
 import numpy as np
 from scipy.special import loggamma
 
@@ -8,14 +9,32 @@ from ..utils import nuclear_radius
 
 
 class FermiFunction(SpectrumComponent):
-    def __init__(self, Z: int, A: int):
+    def __init__(self, Z: int, A: int, logger: logging.Logger | None = None):
+        super().__init__(logger=logger)
         self.Z = Z
         self.A = A
         self.alphaZ = ALPHA * Z
         self.gamma = np.sqrt(1.0 - self.alphaZ**2)
         self.R = nuclear_radius(A)
+        if self._logger:
+            self._logger.debug(
+                "FermiFunction: Z=%d, A=%d, alpha*Z=%.4f, gamma=%.6f, R=%.4f",
+                Z,
+                A,
+                self.alphaZ,
+                self.gamma,
+                self.R,
+            )
 
     def __call__(self, W: np.ndarray) -> np.ndarray:
+        if self._logger:
+            self._logger.debug(
+                "FermiFunction: evaluating at %d points, W range=[%.4f, %.4f]",
+                len(W),
+                W.min(),
+                W.max(),
+            )
+
         p = np.sqrt(W**2 - 1.0)
 
         # Avoid division by zero
@@ -43,4 +62,13 @@ class FermiFunction(SpectrumComponent):
 
         logF = log_prefactor + log_power + log_gamma_abs_sq - log_gamma_denom + log_exp
 
-        return np.asarray(np.exp(logF), dtype=np.float64)
+        result = np.asarray(np.exp(logF), dtype=np.float64)
+
+        if self._logger:
+            self._logger.debug(
+                "FermiFunction: output range=[%.6e, %.6e]",
+                result.min(),
+                result.max(),
+            )
+
+        return result

beta_spectrum/components/finite_size.py

@@ -1,3 +1,5 @@
+import logging
+
 import numpy as np
 
 from beta_spectrum.constants import ALPHA
@@ -13,7 +15,8 @@ class FiniteSizeL0(SpectrumComponent):
     Includes (1 + gamma)/2 prefactor
     """
 
-    def __init__(self, Z: int, A: int):
+    def __init__(self, Z: int, A: int, logger: logging.Logger | None = None):
+        super().__init__(logger=logger)
         self.Z = Z
         self.A = A
 
@@ -25,20 +28,74 @@ class FiniteSizeL0(SpectrumComponent):
         # beta-minus decay
         self.sign = -1
 
+        if self._logger:
+            self._logger.debug(
+                "FiniteSizeL0: Z=%d, A=%d, R=%.4f, alphaZ=%.4f, gamma=%.6f",
+                Z,
+                A,
+                self.R,
+                self.alphaZ,
+                self.gamma,
+            )
+
     def __call__(self, W: np.ndarray) -> np.ndarray:
+        if self._logger:
+            self._logger.debug(
+                "FiniteSizeL0: evaluating at %d points, W range=[%.4f, %.4f]",
+                len(W),
+                W.min(),
+                W.max(),
+            )
+
         # Terms from expansion
         term1 = self.sign * self.alphaZ * W * self.R
         term2 = (13 / 60) * self.alphaZ**2
         term3 = -0.5 * (self.alphaZ * self.R) / W
 
-        return 1 + term1 + term2 + term3
+        result = 1 + term1 + term2 + term3
+
+        if self._logger:
+            self._logger.debug(
+                "FiniteSizeL0: output range=[%.6e, %.6e]",
+                result.min(),
+                result.max(),
+            )
+
+        return result
 
 
 class ChargeDistributionU(SpectrumComponent):
-    def __init__(self, Z: int, A: int):
+    def __init__(self, Z: int, A: int, logger: logging.Logger | None = None):
+        super().__init__(logger=logger)
         self.Z = Z
         self.R = nuclear_radius(A)
         self.alphaZ = ALPHA * Z
 
+        if self._logger:
+            self._logger.debug(
+                "ChargeDistributionU: Z=%d, A=%d, R=%.4f, alphaZ=%.4f",
+                Z,
+                A,
+                self.R,
+                self.alphaZ,
+            )
+
     def __call__(self, W: np.ndarray) -> np.ndarray:
-        return 1.0 + (1 / 5) * (self.alphaZ * W * self.R) ** 2
+        if self._logger:
+            self._logger.debug(
+                "ChargeDistributionU: evaluating at %d points, W range=[%.4f, %.4f]",
+                len(W),
+                W.min(),
+                W.max(),
+            )
+
+        result = 1.0 + (1 / 5) * (self.alphaZ * W * self.R) ** 2
+
+        if self._logger:
+            self._logger.debug(
+                "ChargeDistributionU: output range=[%.6e, %.6e]",
+                result.min(),
+                result.max(),
+            )
+
+        return result

beta_spectrum/components/phase_space.py

@@ -15,6 +15,7 @@ class PhaseSpace(SpectrumComponent):
         transition_type: str = "A",
         m_e: float = 1.0,
         m_nu: float = 0.0,
+        logger=None,
     ):
         """
         Create instance of phase space with given decay parameters.
@@ -25,8 +26,18 @@ class PhaseSpace(SpectrumComponent):
             transition_type: allowed or N-order forbidden, should be [A, F1, F1U, F2, F2U, F3, F3U, F4]
             m_e: 1st body mass (electron by default: m_e = 1)
             m_nu:2nd body mass (massless neutrino by default: m_nu = 0)
+            logger: Optional logger for debug/info output
         """
+        super().__init__(logger=logger)
         self.tr_type = transition_type
+        if self._logger:
+            self._logger.debug(
+                "PhaseSpace: W0=%.5f, transition_type=%s, m_e=%.6f, m_nu=%.6f",
+                W0,
+                transition_type,
+                m_e,
+                m_nu,
+            )
         if self.tr_type not in ["A", "F1", "F1U", "F2", "F2U", "F3", "F3U", "F4"]:
             raise ValueError(
                 f"Cannot calculate phase space for forbiddance '{self.tr_type}'.\nShould be [A, F1, F1U, F2, F2U, F3, F3U, F4].\n"
@@ -44,17 +55,21 @@ class PhaseSpace(SpectrumComponent):
         phase_space = p_e * w_e * p_nu * w_nu
 
         if self.tr_type in ["A", "F1"]:
-            forbid_factor = 1
-
+            forbid_factor: float | np.ndarray = 1
         elif self.tr_type in ["F1U", "F2"]:
             forbid_factor = p_nu**2 + p_e**2
-
         elif self.tr_type in ["F2U", "F3"]:
             forbid_factor = p_nu**4 + 3 / 10 * p_nu**2 * p_e**2 + p_e**4
-
-        elif self.tr_type in ["F3U", "F4"]:
+        else:  # F3U, F4
             forbid_factor = (
                 p_nu**6 + 7 * p_nu**4 * p_e**2 + 7 * p_nu**2 * p_e**4 + p_e**6
             )
 
-        return np.asarray(phase_space * forbid_factor, dtype=np.float64)
+        result = np.asarray(phase_space * forbid_factor, dtype=np.float64)
+        if self._logger:
+            self._logger.debug(
+                "PhaseSpace: output range=[%.6e, %.6e]",
+                result.min(),
+                result.max(),
+            )
+        return result

beta_spectrum/components/radiative.py

@@ -1,4 +1,6 @@
 # components/radiative.py
+import logging
+
 import numpy as np
 from scipy.special import spence
 
@@ -41,6 +43,7 @@ class RadiativeCorrection(SpectrumComponent):
         A: int | None = None,
         use_endpoint_resummation: bool = True,
         delta_cut: float = 1e-3,
+        logger: logging.Logger | None = None,
     ):
         """
         Args:
@@ -49,7 +52,9 @@ class RadiativeCorrection(SpectrumComponent):
             A: Mass number (optional, affects nuclear model correction)
             use_endpoint_resummation: Apply Eq.(53) to handle ln(W0 - W) divergence
             delta_cut: Only apply resummation when (W0 - W) < delta_cut (spec Section 3.4)
+            logger: Optional logger for debug/info output
         """
+        super().__init__(logger=logger)
         self.W0 = W0
         self.Z = Z
         self.A = A
@@ -65,6 +70,16 @@ class RadiativeCorrection(SpectrumComponent):
         # Precompute O(Z*alpha^2) energy-independent nuclear part if A is given
         self._delta_F = self._compute_nuclear_model_correction(Z, A)
 
+        if self._logger:
+            self._logger.debug(
+                "RadiativeCorrection: W0=%.4f, Z=%d, A=%s, delta_cut=%.2e, resummation=%s",
+                W0,
+                Z,
+                A,
+                delta_cut,
+                use_endpoint_resummation,
+            )
+
     def __call__(self, W: np.ndarray) -> np.ndarray:
         """
         Calculate R(W, W0) = 1 + delta_1(W, W0) + delta_2(Z, W).
@@ -75,6 +90,15 @@ class RadiativeCorrection(SpectrumComponent):
         Returns:
             Radiative correction factor array
         """
+        if self._logger:
+            self._logger.debug(
+                "RadiativeCorrection: evaluating at %d points, W range=[%.4f, %.4f], W0=%.4f",
+                len(W),
+                W.min(),
+                W.max(),
+                self.W0,
+            )
+
         W = np.asarray(W, dtype=float)
 
         # Avoid numerical issues exactly at endpoint W0
@@ -95,7 +119,16 @@ class RadiativeCorrection(SpectrumComponent):
             delta_R_finite = delta_1 + delta_2
             delta_R[mask_finite] = delta_R_finite
 
-        return np.asarray(1.0 + delta_R, dtype=np.float64)
+        result = np.asarray(1.0 + delta_R, dtype=np.float64)
+
+        if self._logger:
+            self._logger.debug(
+                "RadiativeCorrection: output range=[%.6e, %.6e]",
+                result.min(),
+                result.max(),
+            )
+
+        return result
 
     def _delta_1(self, W: np.ndarray) -> np.ndarray:
         """

beta_spectrum/components/screening.py

@@ -1,5 +1,7 @@
 from typing import Any
 
+import logging
+
 import numpy as np
 
 from beta_spectrum.base import SpectrumComponent
@@ -26,6 +28,7 @@ class ScreeningCorrection(SpectrumComponent):
         Delta: float = 0.01,
         eps: float = 1e-6,
         C: float = 0.015,
+        logger: logging.Logger | None = None,
     ):
         """
         Parameters
@@ -39,6 +42,7 @@ class ScreeningCorrection(SpectrumComponent):
         Delta : switching width
         eps : numerical floor for W_tilde
         """
+        super().__init__(logger=logger)
         self.fermi = fermi_function
         self.Z = fermi_function.Z
         self.Ws = Ws
@@ -47,10 +51,27 @@ class ScreeningCorrection(SpectrumComponent):
         self._C = C
         self.V0 = V0 if V0 is not None else self._estimate_V0()
 
+        if self._logger:
+            self._logger.debug(
+                "ScreeningCorrection: Z=%d, V0=%.6f, Ws=%.4f, Delta=%.4f",
+                self.Z,
+                self.V0,
+                self.Ws,
+                self.Delta,
+            )
+
     # --------------------
     # Core evaluation
     # --------------------
     def __call__(self, W: np.ndarray) -> np.ndarray:
+        if self._logger:
+            self._logger.debug(
+                "ScreeningCorrection: evaluating at %d points, W range=[%.4f, %.4f]",
+                len(W),
+                W.min(),
+                W.max(),
+            )
+
         W = np.asarray(W)
 
         # Energy shift by screening potential V0
@@ -70,7 +91,16 @@ class ScreeningCorrection(SpectrumComponent):
         # Smooth switching
         f = self._switching_function(W)
 
-        return np.asarray(1.0 + f * (S_raw - 1.0), dtype=np.float64)
+        result = np.asarray(1.0 + f * (S_raw - 1.0), dtype=np.float64)
+
+        if self._logger:
+            self._logger.debug(
+                "ScreeningCorrection: output range=[%.6e, %.6e]",
+                result.min(),
+                result.max(),
+            )
+
+        return result
 
     # --------------------
     # Helpers

beta_spectrum/cw_extractor.py

@@ -256,9 +256,9 @@ class CWExtractor:
         if self.config.use_radiative:
             from .components.radiative import RadiativeCorrection
 
-            W0 = T_to_W(self.config.endpoint_MeV)
+            W0_float = float(T_to_W(self.config.endpoint_MeV))
             rad = RadiativeCorrection(
-                W0=W0,
+                W0=W0_float,
                 Z=self.config.Z_parent,
                 A=self.config.A_number,
                 use_endpoint_resummation=True,
@@ -291,9 +291,9 @@ class CWExtractor:
             Extracted C(W) values with uncertainties.
         """
         if endpoint_keV is not None:
-            W0 = T_to_W(endpoint_keV / 1000.0)
+            W0_float = float(T_to_W(endpoint_keV / 1000.0))
         else:
-            W0 = self._determine_endpoint()
+            W0_float = float(self._determine_endpoint())
 
         # Compute theoretical factors (excluding C(W))
         factors = self._theoretical_factors(self.energies_W)
@@ -310,12 +310,12 @@ class CWExtractor:
             product *= factors[key]
 
         # Endpoint factor: (W0 - W)^2
-        endpoint_factor = (W0 - self.energies_W) ** 2
+        endpoint_factor = (W0_float - self.energies_W) ** 2
 
         # Overall normalization
         if normalize_method == "endpoint":
             # Use flux at lowest energy where C(W) ≈ 1 for allowed decay
-            low_E_mask = self.energies_W < W0 - 0.02
+            low_E_mask = self.energies_W < W0_float - 0.02
             if np.any(low_E_mask):
                 C0 = np.median(
                     self.flux[low_E_mask]
@@ -325,7 +325,7 @@ class CWExtractor:
                 C0 = 1.0
         else:
             # Fit normalization
-            C0 = self._fit_normalization(W0, product, endpoint_factor)
+            C0 = self._fit_normalization(W0_float, product, endpoint_factor)
 
         # Extract C(W)
         cw = self.flux / (C0 * product * endpoint_factor)
@@ -347,8 +347,8 @@ class CWExtractor:
             cw_errors=cw_errors,
             flux=self.flux.copy(),
             flux_errors=self.flux_errors.copy(),
-            endpoint_W=W0,
-            endpoint_keV=W_to_T(W0) * 1000.0,
+            endpoint_W=W0_float,
+            endpoint_keV=float(W_to_T(W0_float)) * 1000.0,
             fit_result=None,
         )
 
@@ -585,7 +585,7 @@ class CWExtractor:
             p = momentum(W)
             from .components.fermi import FermiFunction
 
-            F = FermiFunction()(W, self.config.Z_parent)
+            F = FermiFunction(Z=self.config.Z_parent, A=self.config.A_number)(W)
 
             kurie = np.sqrt(flux / (p * W * F))
 
@@ -605,7 +605,7 @@ class CWExtractor:
             except np.linalg.LinAlgError:
                 continue
 
-        return T_to_W(best_endpoint / 1000.0)
+        return float(T_to_W(best_endpoint / 1000.0))
 
     def _fit_kurie_endpoint(self, cw_result: CWExtractionResult) -> FitResult:
         """Fit Kurie plot to determine endpoint energy."""
@@ -620,7 +620,7 @@ class CWExtractor:
         p = momentum(W)
         from .components.fermi import FermiFunction
 
-        F = FermiFunction()(W, self.config.Z_parent)
+        F = FermiFunction(Z=self.config.Z_parent, A=self.config.A_number)(W)
 
         kurie = np.sqrt(flux / (p * W * F))
         E_keV = self.energies_keV[mask]
@@ -639,7 +639,7 @@ class CWExtractor:
         )
 
         return fitter.fit(
-            x0=[np.mean(kurie), -np.mean(kurie) / E_keV[-1]],
+            x0=[float(np.mean(kurie)), float(-np.mean(kurie) / E_keV[-1])],
             param_names=["a", "b"],
         )
 

beta_spectrum/fitter.py

@@ -390,18 +390,18 @@ class CurveFitter:
                 np.asarray(bounds[1], dtype=np.float64),
             )
 
-        ls_kwargs = dict(
-            fun=self._residuals,
-            x0=x0_array,
-            method=self.config.method,
-            max_nfev=self.config.max_nfev,
-            x_scale=self.config.x_scale,
-            loss=self.config.loss,
-            f_scale=self.config.f_scale,
-            ftol=self.config.ftol,
-            xtol=self.config.xtol,
-            gtol=self.config.gtol,
-        )
+        ls_kwargs: dict[str, object] = {
+            "fun": self._residuals,
+            "x0": x0_array,
+            "method": self.config.method,
+            "max_nfev": self.config.max_nfev,
+            "x_scale": self.config.x_scale,
+            "loss": self.config.loss,
+            "f_scale": self.config.f_scale,
+            "ftol": self.config.ftol,
+            "xtol": self.config.xtol,
+            "gtol": self.config.gtol,
+        }
         if bounds is not None:
             ls_kwargs["bounds"] = bounds_array
 

beta_spectrum/logging_utils.py

@@ -137,9 +137,9 @@ def get_git_short_hash(length: int = 7) -> str:
     str
         Short git hash, or "unknown" if not in a git repo.
     """
-    try:
-        import subprocess
+    import subprocess
 
+    try:
         result = subprocess.run(
             ["git", "rev-parse", f"--short={length}", "HEAD"],
             capture_output=True,

beta_spectrum/nuclear_data.py

@@ -309,6 +309,7 @@ def _resolve_decay_index(
 
     e = pe.ENSDF()
     pairs = e.ensdf_pairs(edata, mode)
+    assert pairs is not None, "paceENSDF returned None for ensdf_pairs"
 
     matching = [(k, v) for k, v in pairs.items() if str(k[0]) == parent_symbol]
 
@@ -378,13 +379,31 @@ def get_decay_info_from_paceENSDF(
     jpi_data = e.get_parent_jpi(edata, nuclide, decay_idx, mode=mode)
     beta_data = e.get_beta_minus(edata, nuclide, decay_idx, units="best")
 
-    # Extract Q-value and half-life
-    q_value_keV = 0.0
+    # 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"
-    for k, v in halflife_data.items():
-        q_value_keV = v[0]
-        half_life_str = str(v[2]) if len(v) > 2 else "unknown"
-        break
+
+    if beta_data:
+        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 halflife_data:
+        for k, v in halflife_data.items():
+            half_life_str = str(v[2]) if len(v) > 2 else "unknown"
+            break
 
     # Extract parent spin/parity
     parent_J = jpi_data[0][0] if jpi_data else 0.0
@@ -398,10 +417,19 @@ def get_decay_info_from_paceENSDF(
     forbidden_code = "0A"  # default
     transition_type = "A"  # default
 
+    if beta_data is None:
+        beta_data = {}
+
     for k, v in beta_data.items():
         for state in v:
-            # state[15] contains forbiddenness code (e.g. '2F', '1UF')
-            if len(state) > 15 and state[15] is not None:
+            # Use ground state branch (level_index == 0) for forbiddenness
+            if state[0] == 0 and len(state) > 15 and state[15] is not None:
+                code = str(state[15])
+                if code in FORBIDDENNESS_MAP:
+                    forbidden_code = code
+                    transition_type = FORBIDDENNESS_MAP[code]
+            # Fallback: use first branch's forbiddenness if no ground state
+            elif transition_type == "A" and len(state) > 15 and state[15] is not None:
                 code = str(state[15])
                 if code in FORBIDDENNESS_MAP:
                     forbidden_code = code
@@ -420,8 +448,8 @@ def get_decay_info_from_paceENSDF(
         Z_parent=Z_parent,
         Z_daughter=Z_daughter,
         A_number=A_number,
-        endpoint_keV=q_value_keV,
-        endpoint_MeV=q_value_keV / 1000.0,
+        endpoint_keV=endpoint_keV,
+        endpoint_MeV=endpoint_keV / 1000.0,
         transition_type=transition_type,
         forbiddenness_code=forbidden_code,
         parent_J=parent_J,

beta_spectrum/spectrum.py

@@ -2,9 +2,14 @@
 from __future__ import annotations
 
 import logging
+from datetime import datetime, timezone
+from typing import TYPE_CHECKING
 
 import numpy as np
 import matplotlib.pyplot as plt
+
+if TYPE_CHECKING:
+    from matplotlib.axes import Axes
 from typing import TYPE_CHECKING, Any, Dict, List, Optional, Tuple
 from dataclasses import dataclass
 
@@ -94,6 +99,10 @@ class BetaSpectrum:
         """
         result: np.ndarray = np.ones_like(W, dtype=float)
         for comp in self.components:
+            if self._logger:
+                self._logger.debug(
+                    "Evaluating %s at %d energy points", comp.__class__.__name__, len(W)
+                )
             result = result * comp(W)
         return result.astype(np.float64)
 
@@ -131,29 +140,48 @@ class BetaSpectrum:
         W0 = float(T_to_W(config.endpoint_MeV))
 
         if config.use_phase_space:
-            components.append(PhaseSpace(W0=W0, transition_type=config.transition_type))
+            components.append(
+                PhaseSpace(W0=W0, transition_type=config.transition_type, logger=logger)
+            )
 
         if config.use_fermi:
-            components.append(FermiFunction(Z=config.Z_daughter, A=config.A_number))
+            components.append(
+                FermiFunction(Z=config.Z_daughter, A=config.A_number, logger=logger)
+            )
 
         if config.use_finite_size:
-            components.append(FiniteSizeL0(Z=config.Z_daughter, A=config.A_number))
+            components.append(
+                FiniteSizeL0(Z=config.Z_daughter, A=config.A_number, logger=logger)
+            )
 
         if config.use_charge_dist:
             components.append(
-                ChargeDistributionU(Z=config.Z_daughter, A=config.A_number)
+                ChargeDistributionU(
+                    Z=config.Z_daughter, A=config.A_number, logger=logger
+                )
             )
 
         if config.use_screening:
             components.append(
-                ScreeningCorrection(FermiFunction(Z=config.Z_parent, A=config.A_number))
+                ScreeningCorrection(
+                    FermiFunction(Z=config.Z_parent, A=config.A_number, logger=logger),
+                    logger=logger,
+                )
             )
 
         if config.use_exchange:
-            components.append(ExchangeCorrection(Z=config.Z_parent))
+            components.append(ExchangeCorrection(Z=config.Z_parent, logger=logger))
 
         if config.use_radiative:
-            components.append(RadiativeCorrection(W0=W0, use_endpoint_resummation=True))
+            components.append(
+                RadiativeCorrection(
+                    W0=W0,
+                    Z=config.Z_daughter,
+                    A=config.A_number,
+                    use_endpoint_resummation=True,
+                    logger=logger,
+                )
+            )
 
         return cls(components, logger=logger)
 
@@ -389,13 +417,136 @@ class BetaSpectrumAnalyzer:
 
         return convolved
 
-    def plot_analysis(self, save_path: Optional[str] = None) -> None:
+    def plot_analysis(
+        self,
+        save_path: Optional[str] = None,
+        show_components: bool = True,
+    ) -> None:
         """
-        Create visualization of the spectrum and all correction.
+        Create visualization of the spectrum and all correction factors.
+
+        Parameters
+        ----------
+        save_path : str, optional
+            Path to save the figure.
+        show_components : bool
+            If True, show full debug view (4-panel with all corrections).
+            If False, show only the spectrum plot with nuclear data header.
         """
         total = self.total_spectrum(normalize=True)
         components = self.components
+        commit = get_git_short_hash(6)
+        timestamp = datetime.now(timezone.utc).strftime("%Y-%m-%d %H:%M:%S UTC")
 
+        if show_components:
+            self._plot_debug_view(total, components, save_path, commit, timestamp)
+        else:
+            self._plot_spectrum_only(total, save_path, commit, timestamp)
+
+    def _add_id_textbox(self, ax: "Axes", commit: str, timestamp: str) -> None:
+        """Add ID text box with commit hash and timestamp to plot."""
+        id_text = f"commit: {commit}  |  {timestamp}"
+        ax.text(
+            0.01,
+            0.01,
+            id_text,
+            transform=ax.transAxes,
+            fontsize=8,
+            color="gray",
+            verticalalignment="bottom",
+            horizontalalignment="left",
+            fontfamily="monospace",
+        )
+
+    def _add_nuclear_data_header(self, ax: "Axes") -> None:
+        """Add nuclear data information header to plot."""
+        parent = self._element_symbol(self.config.Z_parent)
+        daughter = self._element_symbol(self.config.Z_daughter)
+
+        enabled_corrections = []
+        if self.config.use_phase_space:
+            enabled_corrections.append("phase_space")
+        if self.config.use_fermi:
+            enabled_corrections.append("fermi")
+        if self.config.use_screening:
+            enabled_corrections.append("screening")
+        if self.config.use_finite_size:
+            enabled_corrections.append("finite_size")
+        if self.config.use_charge_dist:
+            enabled_corrections.append("charge_dist")
+        if self.config.use_radiative:
+            enabled_corrections.append("radiative")
+        if self.config.use_exchange:
+            enabled_corrections.append("exchange")
+
+        header_lines = [
+            f"Nuclide:    {parent}{self.config.A_number} -> {daughter}{self.config.A_number}",
+            f"Endpoint:   {self.config.endpoint_MeV * 1000:.1f} keV",
+            f"Transition: {self.config.transition_type}",
+            f"Corrections: {', '.join(enabled_corrections)}",
+        ]
+
+        if self.config.use_detector_response:
+            header_lines.append(
+                f"Detector:   {self.config.detector_model} "
+                f"(σ={self.config.detector_sigma_a_keV} keV)"
+            )
+
+        header_text = "\n".join(header_lines)
+        ax.text(
+            0.99,
+            0.99,
+            header_text,
+            transform=ax.transAxes,
+            fontsize=8,
+            verticalalignment="top",
+            horizontalalignment="right",
+            bbox=dict(
+                boxstyle="round,pad=0.3", facecolor="white", alpha=0.8, edgecolor="none"
+            ),
+        )
+
+    def _plot_spectrum_only(
+        self,
+        total: np.ndarray,
+        save_path: Optional[str],
+        commit: str,
+        timestamp: str,
+    ) -> None:
+        """Plot only the total spectrum with nuclear data header."""
+        fig, ax = plt.subplots(figsize=(10, 6))
+
+        ax.plot(self.energies_MeV, total, "b-", lw=2, label="Normalized spectrum")
+        ax.set_xlabel(r"Electron kinetic energy $E$ [MeV]", fontsize=11)
+        ax.set_ylabel("Normalized Counts", fontsize=11)
+        ax.set_title(
+            f"Beta-decay: {self.config.Z_parent} -> {self.config.Z_daughter}, A={self.config.A_number}",
+            fontsize=13,
+        )
+        ax.set_yscale("log")
+        ax.grid(True, alpha=0.3)
+        ax.legend(loc="best", fontsize=10)
+
+        self._add_nuclear_data_header(ax)
+        self._add_id_textbox(ax, commit, timestamp)
+
+        plt.tight_layout()
+
+        if save_path:
+            fig.savefig(save_path, dpi=150, bbox_inches="tight")
+            print(f"Figure saved to {save_path}")
+
+        plt.show()
+
+    def _plot_debug_view(
+        self,
+        total: np.ndarray,
+        components: Dict[str, np.ndarray],
+        save_path: Optional[str],
+        commit: str,
+        timestamp: str,
+    ) -> None:
+        """Full debug view with 4-panel spectrum analysis."""
         fig = plt.figure(figsize=(14, 10))
 
         # 1. Main spectrum plot (top-left)
@@ -663,11 +814,120 @@ class BetaSpectrumAnalyzer:
         # Build metadata header
         from datetime import datetime, timezone
 
+        def _element_symbol(Z: int) -> str:
+            symbols = [
+                "",
+                "H",
+                "He",
+                "Li",
+                "Be",
+                "B",
+                "C",
+                "N",
+                "O",
+                "F",
+                "Ne",
+                "Na",
+                "Mg",
+                "Al",
+                "Si",
+                "P",
+                "S",
+                "Cl",
+                "Ar",
+                "K",
+                "Ca",
+                "Sc",
+                "Ti",
+                "V",
+                "Cr",
+                "Mn",
+                "Fe",
+                "Co",
+                "Ni",
+                "Cu",
+                "Zn",
+                "Ga",
+                "Ge",
+                "As",
+                "Se",
+                "Br",
+                "Kr",
+                "Rb",
+                "Sr",
+                "Y",
+                "Zr",
+                "Nb",
+                "Mo",
+                "Tc",
+                "Ru",
+                "Rh",
+                "Pd",
+                "Ag",
+                "Cd",
+                "In",
+                "Sn",
+                "Sb",
+                "Te",
+                "I",
+                "Xe",
+                "Cs",
+                "Ba",
+                "La",
+                "Ce",
+                "Pr",
+                "Nd",
+                "Pm",
+                "Sm",
+                "Eu",
+                "Gd",
+                "Tb",
+                "Dy",
+                "Ho",
+                "Er",
+                "Tm",
+                "Yb",
+                "Lu",
+                "Hf",
+                "Ta",
+                "W",
+                "Re",
+                "Os",
+                "Ir",
+                "Pt",
+                "Au",
+                "Hg",
+                "Tl",
+                "Pb",
+                "Bi",
+                "Po",
+                "At",
+                "Rn",
+                "Fr",
+                "Ra",
+                "Ac",
+                "Th",
+                "Pa",
+                "U",
+                "Np",
+                "Pu",
+                "Am",
+                "Cm",
+                "Bk",
+                "Cf",
+                "Es",
+                "Fm",
+            ]
+            return symbols[Z] if Z < len(symbols) else f"Z{Z}"
+
+        parent_symbol = _element_symbol(self.config.Z_parent)
+        daughter_symbol = _element_symbol(self.config.Z_daughter)
+
         header_lines = [
             f"# beta-spectrum v{__import__('beta_spectrum').__version__}",
             f"# timestamp: {datetime.now(timezone.utc).isoformat()}",
             f"# source: {source_type}",
-            f"# nuclide: {self.config.Z_parent}->{self.config.Z_daughter}, A={self.config.A_number}",
+            f"# nuclide: {parent_symbol}{self.config.A_number} -> {daughter_symbol}{self.config.A_number} (Z={self.config.Z_parent}->{self.config.Z_daughter})",
             f"# endpoint: {self.config.endpoint_MeV * 1000:.1f} keV",
             f"# transition: {self.config.transition_type}",
             f"# corrections: {', '.join(enabled)}",

notebooks/conftest.py

@@ -1,5 +1,6 @@
 # conftest.py — shared fixtures for notebooks
 import matplotlib
+
 matplotlib.use("Agg")  # non-interactive backend for CI
 
 import matplotlib.pyplot as plt

scripts/run_notebooks.py

@@ -98,7 +98,9 @@ def main():
                     failed.append(nb.name)
 
     print()
-    print(f"Done. {len(all_plots)} plot(s) extracted, {len(failed)} notebook(s) failed.")
+    print(
+        f"Done. {len(all_plots)} plot(s) extracted, {len(failed)} notebook(s) failed."
+    )
     if all_plots:
         print()
         print("Plots saved to: notebooks/_plots/")

tests/test_logging_utils.py

@@ -217,7 +217,8 @@ class TestCSVMetadataHeader:
         with open(filename, "r") as f:
             content = f.read()
 
-        assert "# nuclide: 43->44, A=99" in content
+        assert "# nuclide: Tc99 -> Ru99" in content
+        assert "Z=43->44" in content
         os.unlink(filename)
 
     def test_header_contains_git_commit(self):
@@ -316,9 +317,9 @@ class TestCLILogging:
 
         parser = _build_parser()
         args = parser.parse_args(
-            ["--nuclide", "Tc99", "--dry-run", "--log-file", "test.log"]
+            ["--nuclide", "Tc99", "--dry-run", "--log-file", "/tmp/test.log"]
         )
-        assert args.log_file == "test.log"
+        assert args.log_file == "/tmp/test.log"
 
     def test_dry_run_flag(self, capsys):
         """Test that --dry-run is recognized."""
@@ -336,3 +337,27 @@ class TestCLILogging:
         with pytest.raises(SystemExit) as exc_info:
             parser.parse_args(["--version"])
         assert exc_info.value.code == 0
+
+    def test_removed_mode_flag(self):
+        """Test that --mode option has been removed."""
+        from beta_spectrum.cli import _build_parser
+
+        parser = _build_parser()
+        with pytest.raises(SystemExit):
+            parser.parse_args(["--nuclide", "Tc99", "--mode", "beta_minus"])
+
+    def test_removed_transition_type_flag(self):
+        """Test that --transition-type option has been removed."""
+        from beta_spectrum.cli import _build_parser
+
+        parser = _build_parser()
+        with pytest.raises(SystemExit):
+            parser.parse_args(["--nuclide", "Tc99", "--transition-type", "A"])
+
+    def test_removed_decay_index_flag(self):
+        """Test that --decay-index option has been removed."""
+        from beta_spectrum.cli import _build_parser
+
+        parser = _build_parser()
+        with pytest.raises(SystemExit):
+            parser.parse_args(["--nuclide", "Tc99", "--decay-index", "0"])