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