unzhakov
876e1e3ab8
- Add neutrino mass parameter to SpectrumConfig and CLI interface - Update PhaseSpace component to handle neutrino mass effects properly - Update CSV headers to show basic endpoint and all transition types - Add --csv-verbosity option to CLI with 'full' and 'concise' options - Fix comparison plots to properly show linear scale variations - Add comprehensive tests for neutrino mass functionality
97 lines
3.3 KiB
Python
97 lines
3.3 KiB
Python
# components/phase_space.py
|
|
import numpy as np
|
|
from beta_spectrum.base import SpectrumComponent
|
|
from beta_spectrum.utils import momentum
|
|
|
|
|
|
class PhaseSpace(SpectrumComponent):
|
|
"""
|
|
3-body phase space calculation (assuming heavy nucleus and no recoil)
|
|
"""
|
|
|
|
def __init__(
|
|
self,
|
|
W0: float,
|
|
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.
|
|
|
|
Args:
|
|
|
|
W0: Endpoint energy in natural units m_e=c=h_bar=1
|
|
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._logger:
|
|
self._logger.debug(
|
|
"PhaseSpace: W0=%.5f, transition_type=%s, m_e=%.6f, m_nu=%.6f",
|
|
W0,
|
|
transition_type,
|
|
m_e,
|
|
m_nu,
|
|
)
|
|
if self.tr_type not in ["A", "F1", "F1U", "F2", "F2U", "F3", "F3U", "F4"]:
|
|
raise ValueError(
|
|
f"Cannot calculate phase space for forbiddance '{self.tr_type}'.\nShould be [A, F1, F1U, F2, F2U, F3, F3U, F4].\n"
|
|
)
|
|
self.W0 = W0
|
|
self.m_e = m_e
|
|
self.m_nu = m_nu
|
|
|
|
def __call__(self, W: np.ndarray) -> np.ndarray:
|
|
w_e = W
|
|
p_e = momentum(W, self.m_e)
|
|
w_nu = self.W0 - W
|
|
|
|
# Calculate p_nu, handling the case where w_nu < self.m_nu
|
|
# In this case, the neutrino energy is less than its mass, which is unphysical
|
|
# The phase space should be 0 in this region
|
|
if self.m_nu > 0:
|
|
# Check where w_nu < self.m_nu
|
|
below_threshold = w_nu < self.m_nu
|
|
# Calculate p_nu only where it's physically valid
|
|
p_nu = np.zeros_like(w_nu)
|
|
valid_mask = ~below_threshold
|
|
if np.any(valid_mask):
|
|
p_nu[valid_mask] = momentum(w_nu[valid_mask], self.m_nu)
|
|
else:
|
|
below_threshold = np.zeros_like(w_nu, dtype=bool)
|
|
p_nu = momentum(w_nu, self.m_nu)
|
|
|
|
# Calculate phase space factor
|
|
phase_space = p_e * w_e * p_nu * w_nu
|
|
|
|
# Calculate forbid factor
|
|
if self.tr_type in ["A", "F1"]:
|
|
forbid_factor: float | np.ndarray = 1
|
|
elif self.tr_type in ["F1U", "F2"]:
|
|
forbid_factor = p_nu**2 + p_e**2
|
|
elif self.tr_type in ["F2U", "F3"]:
|
|
forbid_factor = p_nu**4 + 3 / 10 * p_nu**2 * p_e**2 + p_e**4
|
|
else: # F3U, F4
|
|
forbid_factor = (
|
|
p_nu**6 + 7 * p_nu**4 * p_e**2 + 7 * p_nu**2 * p_e**4 + p_e**6
|
|
)
|
|
|
|
result = np.asarray(phase_space * forbid_factor, dtype=np.float64)
|
|
|
|
# For regions below the neutrino mass threshold, the phase space should be 0
|
|
if self.m_nu > 0:
|
|
result[below_threshold] = 0.0
|
|
|
|
if self._logger:
|
|
self._logger.debug(
|
|
"PhaseSpace: output range=[%.6e, %.6e]",
|
|
result.min(),
|
|
result.max(),
|
|
)
|
|
return result
|