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

@@ -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`

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)
 

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 = R0_FM * A ** (1.0 / 3.0)  # rms charge radius in fm
-        a_fm = r_0 * 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

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
+ME_MEV = 0.510998950  # electron rest mass in MeV / c^2
-MP_MEV = 938.272088816  # 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

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
+        sigma_a_me = config.detector_sigma_a_keV / ME_MEV
-        tau_me = config.detector_tau_keV / ME_ME
+
+        tau_me = config.detector_tau_keV / ME_MEV
 
         W0 = T_to_W(config.endpoint_MeV)
         channel_range = config.detector_channel_energy_range

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:
+def nuclear_radius(A: float, r0: float = R0_FM) -> float:
-    r0_fm = 1.2  # fm
+    R_fm = r0 * A ** (1 / 3)
-    R_fm = r0_fm * 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:
+def T_to_W(T_MeV: float | np.ndarray) -> np.ndarray:
-    return T_keV / ME_MEV + 1.0
+    return T_MeV / ME_MEV + 1.0