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
This commit is contained in:
@@ -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:
|
||||
|
||||
+35
-52
@@ -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()
|
||||
|
||||
@@ -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
|
||||
|
||||
@@ -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
|
||||
|
||||
@@ -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
|
||||
|
||||
@@ -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(
|
||||
|
||||
@@ -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:
|
||||
"""
|
||||
|
||||
@@ -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
|
||||
|
||||
+134
-8
@@ -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)}",
|
||||
|
||||
Reference in New Issue
Block a user