From 330ebddaef831ace65721180b3012b539f581ce6 Mon Sep 17 00:00:00 2001 From: Evgeniy Unzhakov Date: Tue, 5 May 2026 09:08:08 +0300 Subject: [PATCH] refactor: refine CLI options and enhance debug logging - Remove --mode, --transition-type, and --decay-index CLI flags (parameters now come from ENSDF data or custom JSON input) - Clarify energy units in all CLI help text (MeV throughout, keV for detector params) - Enhance logging verbosity: - -v (INFO): shows source, component initialization, workflow steps - -vv (DEBUG): shows detailed internals of all spectrum components (parameters, evaluation counts, input/output ranges) - Improve CSV headers with element notation (e.g. Tc99 -> Ru99) - Add CONVENTIONS.md documenting energy unit and logging conventions - Update TODO.md with completion status --- CONVENTIONS.md | 67 ++++++++++- TODO.md | 2 +- beta_spectrum/base.py | 8 ++ beta_spectrum/cli.py | 87 ++++++--------- beta_spectrum/components/exchange.py | 39 ++++++- beta_spectrum/components/fermi.py | 32 +++++- beta_spectrum/components/finite_size.py | 65 ++++++++++- beta_spectrum/components/phase_space.py | 3 + beta_spectrum/components/radiative.py | 35 +++++- beta_spectrum/components/screening.py | 32 +++++- beta_spectrum/spectrum.py | 142 ++++++++++++++++++++++-- notebooks/conftest.py | 1 + scripts/run_notebooks.py | 4 +- tests/test_logging_utils.py | 31 +++++- 14 files changed, 470 insertions(+), 78 deletions(-) diff --git a/CONVENTIONS.md b/CONVENTIONS.md index 2d01cbe..e0a295b 100644 --- a/CONVENTIONS.md +++ b/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: diff --git a/TODO.md b/TODO.md index 6433a35..db6ffad 100644 --- a/TODO.md +++ b/TODO.md @@ -205,6 +205,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). **Next immediate step:** A3 — data processing pipeline for experimental spectra (background subtraction, energy calibration, dead-time correction, pulse pile-up correction). diff --git a/beta_spectrum/base.py b/beta_spectrum/base.py index ed784f4..b48d623 100644 --- a/beta_spectrum/base.py +++ b/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: diff --git a/beta_spectrum/cli.py b/beta_spectrum/cli.py index 7987b95..0564cb9 100644 --- a/beta_spectrum/cli.py +++ b/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 @@ -42,8 +42,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 +71,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 +102,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 +134,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 +146,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", @@ -228,12 +211,10 @@ def _run(args: argparse.Namespace) -> None: source_type: str 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 +226,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 +240,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,13 +256,13 @@ 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 = BetaSpectrumAnalyzer(spectrum, config, logger=logger) analyzer.plot_analysis(save_path=args.plot) logger.info("Analysis plot saved to %s", args.plot) @@ -322,7 +293,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 +418,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() diff --git a/beta_spectrum/components/exchange.py b/beta_spectrum/components/exchange.py index 73d5c6c..1e75d7c 100644 --- a/beta_spectrum/components/exchange.py +++ b/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 diff --git a/beta_spectrum/components/fermi.py b/beta_spectrum/components/fermi.py index d63f6b5..bb48972 100644 --- a/beta_spectrum/components/fermi.py +++ b/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 diff --git a/beta_spectrum/components/finite_size.py b/beta_spectrum/components/finite_size.py index 8648e34..47d8b25 100644 --- a/beta_spectrum/components/finite_size.py +++ b/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 diff --git a/beta_spectrum/components/phase_space.py b/beta_spectrum/components/phase_space.py index 5eaee28..669cd42 100644 --- a/beta_spectrum/components/phase_space.py +++ b/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,7 +26,9 @@ 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.tr_type not in ["A", "F1", "F1U", "F2", "F2U", "F3", "F3U", "F4"]: raise ValueError( diff --git a/beta_spectrum/components/radiative.py b/beta_spectrum/components/radiative.py index 176a626..2c979a2 100644 --- a/beta_spectrum/components/radiative.py +++ b/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: """ diff --git a/beta_spectrum/components/screening.py b/beta_spectrum/components/screening.py index 9571ac5..12d9d40 100644 --- a/beta_spectrum/components/screening.py +++ b/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 diff --git a/beta_spectrum/spectrum.py b/beta_spectrum/spectrum.py index e59107a..22adb19 100644 --- a/beta_spectrum/spectrum.py +++ b/beta_spectrum/spectrum.py @@ -94,6 +94,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 +135,42 @@ 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, use_endpoint_resummation=True, logger=logger) + ) return cls(components, logger=logger) @@ -663,11 +680,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)}", diff --git a/notebooks/conftest.py b/notebooks/conftest.py index 0d84d19..aaaf6dd 100644 --- a/notebooks/conftest.py +++ b/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 diff --git a/scripts/run_notebooks.py b/scripts/run_notebooks.py index 7ffd4c2..08e9cea 100644 --- a/scripts/run_notebooks.py +++ b/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/") diff --git a/tests/test_logging_utils.py b/tests/test_logging_utils.py index 6b3c429..78399a4 100644 --- a/tests/test_logging_utils.py +++ b/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"])