From 8b4545fe3ed4f9627ba556c37667975b2e9d680e Mon Sep 17 00:00:00 2001 From: Evgeniy Unzhakov Date: Fri, 1 May 2026 22:34:15 +0300 Subject: [PATCH] refactor: centralize all physical constants in constants.py, fix unit naming - Add EULER_GAMMA, ME_KEV, R0_FM to constants.py - Replace hardcoded 0.511, 0.5772..., 1.2 with imports from constants - Fix T_to_W parameter name: T_kev -> T_MeV (it takes MeV, not keV) - Remove duplicate ME_ME from spectrum.py (import ME_MEV from constants) - Remove hardcoded r_0 from detector_response.py and radiative.py - Update CONVENTIONS.md to reference correct constant names --- CONVENTIONS.md | 2 +- beta_spectrum/components/detector_response.py | 4 +++- beta_spectrum/components/radiative.py | 14 ++++++++++---- beta_spectrum/constants.py | 7 +++++-- beta_spectrum/spectrum.py | 8 ++++---- beta_spectrum/utils.py | 11 +++++------ 6 files changed, 28 insertions(+), 18 deletions(-) diff --git a/CONVENTIONS.md b/CONVENTIONS.md index b36461a..2d01cbe 100644 --- a/CONVENTIONS.md +++ b/CONVENTIONS.md @@ -59,7 +59,7 @@ All conversions and constants must go through the central utilities: ```python from beta_spectrum.utils import T_to_W, W_to_T -from beta_spectrum.constants import ME_ME, ME_KEV, ALPHA, M_P +from beta_spectrum.constants import ME_MEV, ME_KEV, ALPHA, MP_MEV ``` - Never hardcode physical constants — import from `constants.py` diff --git a/beta_spectrum/components/detector_response.py b/beta_spectrum/components/detector_response.py index 1d36043..e9e1b85 100644 --- a/beta_spectrum/components/detector_response.py +++ b/beta_spectrum/components/detector_response.py @@ -33,6 +33,8 @@ from scipy.interpolate import RegularGridInterpolator, interp1d from scipy.special import erfc from typing import Optional +from beta_spectrum.constants import ME_MEV + class DetectorResponse: """ @@ -196,7 +198,7 @@ class DetectorResponse: sigma_sq = self.sigma_a**2 + (self.sigma_b * np.sqrt(E)) ** 2 if self.fano_factor > 0: # Fano contribution: σ = sqrt(F * E * w) in natural units - w_me = 3.6e-3 / 0.511 # 3.6 eV in m_e units + w_me = 3.6e-3 / ME_MEV # 3.6 eV in m_e units sigma_sq += self.fano_factor * E * w_me return np.sqrt(sigma_sq) diff --git a/beta_spectrum/components/radiative.py b/beta_spectrum/components/radiative.py index a1e1dd6..7564562 100644 --- a/beta_spectrum/components/radiative.py +++ b/beta_spectrum/components/radiative.py @@ -3,7 +3,14 @@ import numpy as np from scipy.special import spence from beta_spectrum.base import SpectrumComponent -from beta_spectrum.constants import MP_MEV, ME_MEV, ALPHA, HBAR_C_MEV_FM +from beta_spectrum.constants import ( + MP_MEV, + ME_MEV, + ALPHA, + HBAR_C_MEV_FM, + EULER_GAMMA, + R0_FM, +) class RadiativeCorrection(SpectrumComponent): @@ -49,7 +56,7 @@ class RadiativeCorrection(SpectrumComponent): self.use_endpoint_resummation = use_endpoint_resummation self.delta_cut = delta_cut self.m_p = MP_MEV / ME_MEV # ~1836.15 (proton mass in m_e units) - self.gamma_E = 0.57721566490153286 # Euler-Mascheroni constant + self.gamma_E = EULER_GAMMA # Euler-Mascheroni constant # Precompute constant terms self._outer_normalization = ALPHA / (2.0 * np.pi) # ~0.00116 @@ -206,8 +213,7 @@ class RadiativeCorrection(SpectrumComponent): return 0.0 # Nuclear radius from Wilkinson's r_0 - r_0 = 1.2 # fm - a_fm = r_0 * A ** (1.0 / 3.0) # rms charge radius in fm + a_fm = R0_FM * A ** (1.0 / 3.0) # rms charge radius in fm a_me = a_fm / (HBAR_C_MEV_FM / ME_MEV) # convert to m_e^{-1} units # Lambda = sqrt(6)/a diff --git a/beta_spectrum/constants.py b/beta_spectrum/constants.py index f5be66d..486a38e 100644 --- a/beta_spectrum/constants.py +++ b/beta_spectrum/constants.py @@ -1,5 +1,8 @@ # constants.py ALPHA = 1 / 137.035999084 # fine-structure constant +EULER_GAMMA = 0.57721566490153286 # Euler-Mascheroni constant HBAR_C_MEV_FM = 197.3269804 # MeV * fm -ME_MEV = 0.510998950 # MeV / c^2 -MP_MEV = 938.272088816 # MeV / c^2 +ME_MEV = 0.510998950 # electron rest mass in MeV / c^2 +ME_KEV = ME_MEV * 1000.0 # electron rest mass in keV +MP_MEV = 938.272088816 # proton rest mass in MeV / c^2 +R0_FM = 1.2 # nuclear radius parameter in fm diff --git a/beta_spectrum/spectrum.py b/beta_spectrum/spectrum.py index ddc1687..8bfd3cc 100644 --- a/beta_spectrum/spectrum.py +++ b/beta_spectrum/spectrum.py @@ -13,12 +13,11 @@ from beta_spectrum.components.exchange import ExchangeCorrection from beta_spectrum.components.radiative import RadiativeCorrection from beta_spectrum.utils import T_to_W +from beta_spectrum.constants import ME_MEV if TYPE_CHECKING: from beta_spectrum.components.detector_response import DetectorResponse -ME_ME = 0.510998950 # electron rest mass in MeV - @dataclass class SpectrumConfig: @@ -165,8 +164,9 @@ class BetaSpectrum: from beta_spectrum.components.detector_response import DetectorResponse from beta_spectrum.utils import T_to_W - sigma_a_me = config.detector_sigma_a_keV / ME_ME - tau_me = config.detector_tau_keV / ME_ME + sigma_a_me = config.detector_sigma_a_keV / ME_MEV + + tau_me = config.detector_tau_keV / ME_MEV W0 = T_to_W(config.endpoint_MeV) channel_range = config.detector_channel_energy_range diff --git a/beta_spectrum/utils.py b/beta_spectrum/utils.py index a7f326e..4af6047 100644 --- a/beta_spectrum/utils.py +++ b/beta_spectrum/utils.py @@ -1,12 +1,11 @@ # utils.py import numpy as np -from .constants import HBAR_C_MEV_FM, ME_MEV +from .constants import HBAR_C_MEV_FM, ME_MEV, R0_FM -def nuclear_radius(A: float, r0: float = 1.2) -> float: - r0_fm = 1.2 # fm - R_fm = r0_fm * A ** (1 / 3) +def nuclear_radius(A: float, r0: float = R0_FM) -> float: + R_fm = r0 * A ** (1 / 3) return R_fm / (HBAR_C_MEV_FM / ME_MEV) @@ -33,5 +32,5 @@ def W_to_T(W: float | np.ndarray) -> np.ndarray: return (W - 1.0) * ME_MEV -def T_to_W(T_keV: float | np.ndarray) -> np.ndarray: - return T_keV / ME_MEV + 1.0 +def T_to_W(T_MeV: float | np.ndarray) -> np.ndarray: + return T_MeV / ME_MEV + 1.0