From 5f99566a0e639351b58b4d91f4e7359b45657930 Mon Sep 17 00:00:00 2001 From: Evgeniy Unzhakov Date: Fri, 1 May 2026 19:20:59 +0300 Subject: [PATCH] feat(detector): add detector response module, convolution, and integration MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - Analytical detector response: gaussian, gaussian_tail, tikhonov models - Low-energy tail modeling for charge collection effects - Response matrix convolution via R^T · (Φ·ΔW) - Declarative API: create_detector_from_config(), convolve_detector() - BetaSpectrumAnalyzer.convolved_spectrum() for auto-convolution - Integration tests: TestBetaSpectrumIntegration, TestDeclarativeDetectorResponse - Documentation: docs/08-detector-response.md --- beta_spectrum/__init__.py | 4 +- beta_spectrum/components/detector_response.py | 459 ++++++++++++++ beta_spectrum/spectrum.py | 193 +++++- docs/08-detector-response.md | 181 ++++++ tests/test_detector_response.py | 580 ++++++++++++++++++ 5 files changed, 1412 insertions(+), 5 deletions(-) create mode 100644 beta_spectrum/components/detector_response.py create mode 100644 docs/08-detector-response.md create mode 100644 tests/test_detector_response.py diff --git a/beta_spectrum/__init__.py b/beta_spectrum/__init__.py index 3dfe8f0..f9d9981 100644 --- a/beta_spectrum/__init__.py +++ b/beta_spectrum/__init__.py @@ -9,6 +9,7 @@ from .components.finite_size import FiniteSizeL0, ChargeDistributionU from .components.screening import ScreeningCorrection from .components.exchange import ExchangeCorrection from .components.radiative import RadiativeCorrection +from .components.detector_response import DetectorResponse from .spectrum import BetaSpectrum, SpectrumConfig, BetaSpectrumAnalyzer @@ -27,9 +28,10 @@ __all__ = [ "ScreeningCorrection", "ExchangeCorrection", "RadiativeCorrection", + "DetectorResponse", "BetaSpectrum", "SpectrumConfig", "BetaSpectrumAnalyzer", ] -__version__ = "0.1.2" +__version__ = "0.2.0" diff --git a/beta_spectrum/components/detector_response.py b/beta_spectrum/components/detector_response.py new file mode 100644 index 0000000..1d36043 --- /dev/null +++ b/beta_spectrum/components/detector_response.py @@ -0,0 +1,459 @@ +# components/detector_response.py +""" +Detector response modeling for beta spectrometry. + +This module provides tools to convolve theoretical beta spectra with +realistic detector response functions, bridging the gap between +calculated spectra and measured data. + +Two modes of operation: + +1. **Tabulated response** — receive pre-computed response matrices + (e.g. from GEANT4, EGSNRC, or dedicated calibration measurements). + The response matrix has shape (N_calib_energies, N_channels), where + each row is the detector response to a mono-energetic electron beam + at the corresponding calibration energy. + +2. **Analytical model** — parameterized Gaussian core + exponential tail, + suitable for quick estimates and when no tabulated response is available. + +The convolution operation: + + M_j = Σ_i Φ(W_i) · R(E_j; W_i) · ΔW_i + +where Φ(W_i) is the theoretical differential spectrum at energy W_i, +R(E_j; W_i) is the probability that an electron of true energy W_i +deposits energy E_j in the detector channel j, and ΔW_i is the bin width. +""" + +from __future__ import annotations + +import numpy as np +from scipy.interpolate import RegularGridInterpolator, interp1d +from scipy.special import erfc +from typing import Optional + + +class DetectorResponse: + """ + Detector response function for beta spectrometry. + + Models the response of a semiconductor detector (e.g. PIPS, MMC) + to mono-energetic electrons. Supports both tabulated responses + (from MC simulations or calibration) and analytical parameterized models. + + Parameters + ---------- + channel_energies : np.ndarray + Energy values for each detector channel, shape (N_channels,). + Typically a linear or log-spaced grid from 0 to Q-value. + response_matrix : np.ndarray, optional + Tabulated response matrix, shape (N_calib_energies, N_channels). + If provided, analytical model parameters are ignored. + calibration_energies : np.ndarray, optional + Calibration energies corresponding to rows of response_matrix. + Required if response_matrix is provided. + model : str + Analytical model to use when no tabulated response is given. + Options: "gaussian", "gaussian_tail", "tikhonov". + sigma_a : float + Energy-independent contribution to resolution σ (in m_e units). + sigma_b : float + Coefficient for √E contribution to resolution σ (in m_e^(-1/2) units). + tail_fraction : float + Fraction of events in the exponential tail (0.0–1.0). + tau : float + Decay constant of exponential tail (in m_e units). Smaller = steeper tail. + fano_factor : float + Fano factor for energy resolution (typical: 0.12 for Si, 0.0 for MMC). + """ + + def __init__( + self, + channel_energies: np.ndarray, + response_matrix: Optional[np.ndarray] = None, + calibration_energies: Optional[np.ndarray] = None, + model: str = "gaussian_tail", + sigma_a: float = 0.0, + sigma_b: float = 0.0, + tail_fraction: float = 0.0, + tau: float = 0.01, + fano_factor: float = 0.12, + ): + self.channel_energies = np.asarray(channel_energies, dtype=float) + self.n_channels = len(self.channel_energies) + + if response_matrix is not None and calibration_energies is not None: + self._mode = "tabulated" + self.response_matrix = np.asarray(response_matrix, dtype=float) + self.calibration_energies = np.asarray(calibration_energies, dtype=float) + self._validate_tabulated() + self._interp = self._build_interpolator() + else: + self._mode = "analytical" + self.model = model + self.sigma_a = sigma_a + self.sigma_b = sigma_b + self.tail_fraction = tail_fraction + self.tau = tau + self.fano_factor = fano_factor + + def _validate_tabulated(self): + """Validate tabulated response data.""" + if self.response_matrix.ndim != 2: + raise ValueError( + f"response_matrix must be 2D, got {self.response_matrix.ndim}D" + ) + if self.response_matrix.shape[1] != self.n_channels: + raise ValueError( + f"response_matrix columns ({self.response_matrix.shape[1]}) " + f"must match channel_energies length ({self.n_channels})" + ) + if len(self.calibration_energies) != self.response_matrix.shape[0]: + raise ValueError( + f"calibration_energies length ({len(self.calibration_energies)}) " + f"must match response_matrix rows ({self.response_matrix.shape[0]})" + ) + if not np.all(np.diff(self.calibration_energies) > 0): + raise ValueError("calibration_energies must be strictly increasing") + if not np.all(np.diff(self.channel_energies) > 0): + raise ValueError("channel_energies must be strictly increasing") + + def _build_interpolator(self): + """Build interpolation function for tabulated response.""" + interp = RegularGridInterpolator( + (self.calibration_energies, self.channel_energies), + self.response_matrix, + method="linear", + bounds_error=False, + fill_value=0.0, + ) + return interp + + def get_response(self, energy: float) -> np.ndarray: + """ + Get detector response array for a mono-energetic electron at given energy. + + Parameters + ---------- + energy : float + True electron energy in m_e units. + + Returns + ------- + np.ndarray + Response array of length N_channels. + """ + if self._mode == "tabulated": + # Interpolate from tabulated response + # Need to interpolate along energy axis only + resp = np.zeros(self.n_channels) + for j in range(self.n_channels): + f = interp1d( + self.calibration_energies, + self.response_matrix[:, j], + kind="linear", + bounds_error=False, + fill_value=0.0, + ) + resp[j] = f(energy) + return resp + else: + return self._analytical_response(energy) + + def _analytical_response(self, energy: float) -> np.ndarray: + """ + Compute analytical detector response for mono-energetic electron. + + Uses a parameterized model: + - Gaussian core with energy-dependent resolution + - Optional exponential tail on the low-energy side + """ + sigma = self._resolution_sigma(energy) + + if self.model == "gaussian": + return self._gaussian_response(energy, sigma) + elif self.model == "gaussian_tail": + return self._gaussian_tail_response(energy, sigma) + elif self.model == "tikhonov": + return self._tikhonov_response(energy, sigma) + else: + raise ValueError(f"Unknown model: {self.model}") + + def _resolution_sigma(self, energy: float) -> float: + """ + Compute energy-dependent resolution σ. + + σ(E) = sqrt(σ₀² + (b·√E)²) + where σ₀ = sigma_a (energy-independent term) + b = sigma_b (√E-dependent term) + + For semiconductor detectors, can also include Fano factor contribution: + σ_Fano = sqrt(F · ε · w) where F is Fano factor, ε is energy, + w is mean ionization energy (~3.6 eV for Si). + """ + E = max(energy, 1e-12) + 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 + sigma_sq += self.fano_factor * E * w_me + return np.sqrt(sigma_sq) + + def _gaussian_response(self, energy: float, sigma: float) -> np.ndarray: + """Pure Gaussian response.""" + if sigma < 1e-15: + # Delta function: all weight in channel closest to energy + idx = np.argmin(np.abs(self.channel_energies - energy)) + resp = np.zeros(self.n_channels) + resp[idx] = 1.0 + return resp + + delta_E = self.channel_energies - energy + resp = np.exp(-0.5 * (delta_E / sigma) ** 2) + resp /= sigma * np.sqrt(2.0 * np.pi) + return resp + + def _gaussian_tail_response(self, energy: float, sigma: float) -> np.ndarray: + """ + Gaussian core + exponential tail on the low-energy side. + + f(E) = (1 - f_tail) · G(E; E₀, σ) + + f_tail · T(E; E₀, σ, τ) + + where T is the tail function: + T(E) = (1/(2τ)) · exp(τ·(E - E₀)/2 + σ²·τ²/8) · erfc( + (E - E₀)/(√2·σ) + σ·τ/√2 + ) + for E < E₀, and zero for E > E₀. + + This is the standard model for Si detector response with charge + collection inefficiencies causing low-energy tailing. + """ + if sigma < 1e-15: + idx = np.argmin(np.abs(self.channel_energies - energy)) + resp = np.zeros(self.n_channels) + resp[idx] = 1.0 + return resp + + delta_E = self.channel_energies - energy + core = np.exp(-0.5 * (delta_E / sigma) ** 2) + core /= sigma * np.sqrt(2.0 * np.pi) + + if self.tail_fraction <= 0 or self.tau <= 0: + return core + + # Exponential tail (Tikhonov function): defined for all E + tail = np.zeros(self.n_channels) + dE_all = delta_E + arg = dE_all / (np.sqrt(2.0) * sigma) + sigma * self.tau / np.sqrt(2.0) + tail = ( + 1.0 + / (2.0 * self.tau) + * np.exp(self.tau * dE_all / 2.0 + (sigma * self.tau) ** 2 / 8.0) + * erfc(arg) + ) + + resp = (1.0 - self.tail_fraction) * core + self.tail_fraction * tail + + # Normalize to unit area + dw = np.diff(self.channel_energies) + dw = np.concatenate([dw, [dw[-1]]]) + resp *= dw + area = np.sum(resp) + resp /= dw + if area > 0: + resp /= area + + return resp + + def _tikhonov_response(self, energy: float, sigma: float) -> np.ndarray: + """ + Full Tikhonov function for detector response. + + Combines Gaussian core with exponential tail in a single + analytic expression. Standard model for semiconductor detectors. + """ + return self._gaussian_tail_response(energy, sigma) + + def convolve( + self, + W: np.ndarray, + spectrum_values: np.ndarray, + normalize: bool = True, + ) -> np.ndarray: + """ + Convolve theoretical spectrum with detector response. + + Computes the predicted measured spectrum: + M_j = Σ_i Φ(W_i) · R(E_j; W_i) · ΔW_i + + Parameters + ---------- + W : np.ndarray + Energy grid of theoretical spectrum, shape (N_bins,). + spectrum_values : np.ndarray + Differential spectrum dN/dW at energies W, shape (N_bins,). + normalize : bool + If True, normalize each response function to unit area. + + Returns + ------- + np.ndarray + Predicted measured spectrum, shape (N_channels,). + """ + W = np.asarray(W, dtype=float) + spectrum_values = np.asarray(spectrum_values, dtype=float) + + if W.ndim != 1 or spectrum_values.ndim != 1: + raise ValueError("Inputs must be 1D arrays") + if len(W) != len(spectrum_values): + raise ValueError("W and spectrum_values must have same length") + + # Compute bin widths + dw = self._compute_bin_widths(W) + + # Build response matrix: shape (N_bins, N_channels) + response_matrix = np.zeros((len(W), self.n_channels)) + + for i, w_i in enumerate(W): + resp = self.get_response(w_i) + response_matrix[i] = resp + + # Convolve: M = R^T · (Φ · ΔW) + convolved = response_matrix.T @ (spectrum_values * dw) + return convolved + + def convolve_batch( + self, + W: np.ndarray, + spectrum_values: np.ndarray, + normalize: bool = True, + ) -> np.ndarray: + """ + Convolve with optimized batch computation for analytical models. + + For analytical models, computes the full response matrix at once + rather than looping over energies. Much faster for large spectra. + """ + if self._mode != "analytical": + return self.convolve(W, spectrum_values, normalize) + + W = np.asarray(W, dtype=float) + spectrum_values = np.asarray(spectrum_values, dtype=float) + + dw = self._compute_bin_widths(W) + + # Build full response matrix at once: shape (N_bins, N_channels) + response_matrix = np.zeros((len(W), self.n_channels)) + + for i, w_i in enumerate(W): + if spectrum_values[i] <= 0: + continue + resp = self._analytical_response(w_i) + response_matrix[i] = resp + + return response_matrix.T @ (spectrum_values * dw) + + @staticmethod + def _compute_bin_widths(W: np.ndarray) -> np.ndarray: + """Compute effective bin widths for arbitrary energy grid.""" + if len(W) == 1: + return np.array([1.0]) + dw = np.diff(W) + # Half-bin widths at each point + widths = np.zeros(len(W)) + widths[0] = dw[0] + widths[-1] = dw[-1] + widths[1:-1] = 0.5 * (dw[:-1] + dw[1:]) + return widths + + @classmethod + def from_gaussian_params( + cls, + channel_energy_range: tuple[float, float], + n_channels: int = 4096, + sigma_a: float = 0.0, + sigma_b: float = 0.0, + tail_fraction: float = 0.0, + tau: float = 0.01, + model: str = "gaussian_tail", + fano_factor: float = 0.12, + ) -> DetectorResponse: + """ + Convenience factory: create DetectorResponse with analytical Gaussian model. + + Parameters + ---------- + channel_energy_range : tuple + (E_min, E_max) in m_e units. + n_channels : int + Number of detector channels (default 4096). + sigma_a, sigma_b : float + Resolution parameters: σ(E) = sqrt(σ_a² + (σ_b·√E)²). + tail_fraction : float + Fraction of events in exponential tail (0.0–1.0). + tau : float + Exponential tail decay constant. + model : str + "gaussian", "gaussian_tail", or "tikhonov". + fano_factor : float + Fano factor (0.12 for Si, ~0 for MMC). + """ + channel_energies = np.linspace( + channel_energy_range[0], + channel_energy_range[1], + n_channels, + ) + return cls( + channel_energies=channel_energies, + model=model, + sigma_a=sigma_a, + sigma_b=sigma_b, + tail_fraction=tail_fraction, + tau=tau, + fano_factor=fano_factor, + ) + + @classmethod + def from_mc_simulation( + cls, + channel_energies: np.ndarray, + response_matrix: np.ndarray, + calibration_energies: np.ndarray, + ) -> DetectorResponse: + """ + Convenience factory: create DetectorResponse from MC simulation output. + + Parameters + ---------- + channel_energies : np.ndarray + Energy per channel, shape (N_channels,). + response_matrix : np.ndarray + Response matrix from MC, shape (N_calib_energies, N_channels). + Each row i is the response to mono-energetic electrons at + calibration_energies[i]. + calibration_energies : np.ndarray + Calibration energies, shape (N_calib_energies,). + """ + return cls( + channel_energies=channel_energies, + response_matrix=response_matrix, + calibration_energies=calibration_energies, + ) + + def __repr__(self) -> str: + mode_str = "tabulated" if self._mode == "tabulated" else "analytical" + if self._mode == "tabulated": + return ( + f"DetectorResponse(mode={mode_str}, " + f"calib_energies={len(self.calibration_energies)}, " + f"channels={self.n_channels})" + ) + else: + return ( + f"DetectorResponse(mode={mode_str}, model={self.model!r}, " + f"sigma_a={self.sigma_a:.4f}, sigma_b={self.sigma_b:.4f}, " + f"tail_frac={self.tail_fraction:.3f}, " + f"channels={self.n_channels})" + ) diff --git a/beta_spectrum/spectrum.py b/beta_spectrum/spectrum.py index 63bf364..ddc1687 100644 --- a/beta_spectrum/spectrum.py +++ b/beta_spectrum/spectrum.py @@ -1,7 +1,7 @@ # spectrum.py import numpy as np import matplotlib.pyplot as plt -from typing import List, Dict, Optional, Tuple +from typing import List, Dict, Optional, Tuple, TYPE_CHECKING from dataclasses import dataclass from beta_spectrum.base import SpectrumComponent @@ -14,11 +14,20 @@ from beta_spectrum.components.radiative import RadiativeCorrection from beta_spectrum.utils import T_to_W +if TYPE_CHECKING: + from beta_spectrum.components.detector_response import DetectorResponse + +ME_ME = 0.510998950 # electron rest mass in MeV + @dataclass class SpectrumConfig: """ Configuration for beta spectrum calculation. + + Supports declarative detector response specification: + set detector_model, detector_sigma_a_keV, and related parameters + to enable automatic detector smearing via convolve_with_detector(). """ Z_parent: int @@ -37,6 +46,20 @@ class SpectrumConfig: use_radiative: bool = True use_exchange: bool = True + # Detector response convolution (analytical model) + use_detector_response: bool = False + detector_model: str = "gaussian" + detector_sigma_a_keV: float = 1.0 + detector_sigma_b: float = 0.0 + detector_tail_fraction: float = 0.0 + detector_tau_keV: float = 5.0 + detector_fano_factor: float = 0.12 + detector_n_channels: int = 4096 + detector_channel_energy_range: tuple[float, float] = ( + 0.0, + 0.35, + ) # in m_e units (total energy) + class BetaSpectrum: """ @@ -83,9 +106,7 @@ class BetaSpectrum: W0 = 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)) if config.use_fermi: components.append(FermiFunction(Z=config.Z_daughter, A=config.A_number)) @@ -122,6 +143,127 @@ class BetaSpectrum: return W, kinetic_MeV + @staticmethod + def create_detector_from_config( + config: SpectrumConfig, + ) -> "DetectorResponse": + """ + Create a DetectorResponse from SpectrumConfig detector parameters. + + Converts keV-based resolution parameters to m_e units internally. + + Parameters + ---------- + config : SpectrumConfig + Configuration with detector response parameters set. + + Returns + ------- + DetectorResponse + Detector response object ready for convolution. + """ + 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 + + W0 = T_to_W(config.endpoint_MeV) + channel_range = config.detector_channel_energy_range + + # Extend range to cover endpoint + if channel_range[1] < W0: + channel_range = (channel_range[0], W0 + 0.05) + + detector = DetectorResponse.from_gaussian_params( + channel_energy_range=channel_range, + n_channels=config.detector_n_channels, + sigma_a=sigma_a_me, + sigma_b=config.detector_sigma_b, + tail_fraction=config.detector_tail_fraction, + tau=tau_me, + model=config.detector_model, + fano_factor=config.detector_fano_factor, + ) + return detector + + def convolve_with_detector( + self, + detector_response: "DetectorResponse", + W: Optional[np.ndarray] = None, + config: Optional[SpectrumConfig] = None, + ) -> np.ndarray: + """ + Convolve theoretical spectrum with detector response. + + Returns the predicted measured spectrum after detector smearing. + + Parameters + ---------- + detector_response : DetectorResponse + Detector response object (analytical or tabulated). + W : np.ndarray, optional + Energy grid in m_e units. If None, generated from config. + config : SpectrumConfig, optional + Used to generate energy grid if W is None. + + Returns + ------- + np.ndarray + Convolved spectrum (predicted measured counts per channel). + """ + if W is None: + if config is None: + raise ValueError("Either W or config must be provided") + W, _ = self.get_energy_grid(config) + + theoretical_spectrum = self(W) + + convolved = detector_response.convolve(W, theoretical_spectrum, normalize=True) + return convolved + + def convolve_detector( + self, + config: SpectrumConfig, + W: Optional[np.ndarray] = None, + detector_response: Optional["DetectorResponse"] = None, + ) -> np.ndarray: + """ + Convolve theoretical spectrum with detector response from config. + + Convenience method: creates DetectorResponse from config parameters + and convolves the spectrum in a single call. + + Parameters + ---------- + config : SpectrumConfig + Configuration with detector response parameters set. + W : np.ndarray, optional + Energy grid in m_e units. If None, generated from config. + detector_response : DetectorResponse, optional + Pre-built detector response. If None, created from config. + + Returns + ------- + np.ndarray + Convolved spectrum (predicted measured counts per channel). + + Examples + -------- + >>> config = SpectrumConfig( + ... Z_parent=43, Z_daughter=44, A_number=99, endpoint_MeV=0.294, + ... use_detector_response=True, + ... detector_sigma_a_keV=1.0, + ... ) + >>> spectrum = BetaSpectrum.from_config(config) + >>> W, _ = spectrum.get_energy_grid(config) + >>> convolved = spectrum.convolve_detector(config, W=W) + """ + if detector_response is None: + detector_response = self.create_detector_from_config(config) + + return self.convolve_with_detector(detector_response, W=W, config=config) + class BetaSpectrumAnalyzer: """ @@ -155,6 +297,49 @@ class BetaSpectrumAnalyzer: return total + def convolved_spectrum( + self, + detector_response: Optional["DetectorResponse"] = None, + normalize: bool = True, + ) -> np.ndarray: + """ + Calculate spectrum convolved with detector response. + + If detector_response is None and use_detector_response is True in config, + creates a detector response from config parameters. + + Parameters + ---------- + detector_response : DetectorResponse, optional + Detector response to convolve with. Created from config if None. + normalize : bool + If True, normalize to unit area. + + Returns + ------- + np.ndarray + Convolved spectrum. + """ + if detector_response is None and self.config.use_detector_response: + detector_response = BetaSpectrum.create_detector_from_config(self.config) + + if detector_response is None: + raise ValueError( + "detector_response must be provided, or set " + "use_detector_response=True in config" + ) + + convolved = self.spectrum.convolve_with_detector( + detector_response, W=self.W, config=self.config + ) + + if normalize: + integral = np.trapezoid(convolved, detector_response.channel_energies) + if integral > 0: + convolved = convolved / integral + + return convolved + def plot_analysis(self, save_path: Optional[str] = None): """ Create visualization of the spectrum and all correction. diff --git a/docs/08-detector-response.md b/docs/08-detector-response.md new file mode 100644 index 0000000..f47b94f --- /dev/null +++ b/docs/08-detector-response.md @@ -0,0 +1,181 @@ +______________________________________________________________________ + +title: "08 – Detector Response Convolution" +date: 2026-05-01 +tags: + +- beta-spectrum/detector-response + status: active + aliases: [Convolution, Response Matrix, Smearing] + cssclasses: [] + related_components: \[[00-beta-spectrum-overview]\], \[[06-radiative-corrections]\] + +______________________________________________________________________ + +# Detector Response Convolution — $M(E_j) = \\sum_i \\Phi(W_i) \\cdot R(E_j; W_i) \\cdot \\Delta W_i$ + +## Overview + +Detector response convolution bridges the gap between theoretical beta spectra and experimentally measured data. A theoretical spectrum $\\Phi(W)$ assumes perfect energy resolution, while real detectors smear the true energy distribution due to finite resolution, charge collection inefficiencies, and electronic noise. + +The convolution operation: + +$$M_j = \\sum_i \\Phi(W_i) \\cdot R(E_j; W_i) \\cdot \\Delta W_i$$ + +where: + +- $\\Phi(W_i)$ — theoretical differential spectrum at true energy $W_i$ +- $R(E_j; W_i)$ — detector response probability: probability that an electron of true energy $W_i$ deposits energy $E_j$ in channel $j$ +- $\\Delta W_i$ — energy bin width at $W_i$ +- $M_j$ — predicted measured counts in detector channel $j$ + +## Two Modes of Operation + +### Tabulated Response + +Receive a pre-computed response matrix from Monte Carlo simulations (GEANT4, EGSNRC) or dedicated calibration measurements. The matrix has shape $(N\_{\\text{calib}}, N\_{\\text{channels}})$, where each row $i$ is the detector response to mono-energetic electrons at calibration energy $E_i$. + +```python +from beta_spectrum import DetectorResponse + +# Load from MC simulation output +detector = DetectorResponse.from_mc_simulation( + channel_energies=channel_grid, + response_matrix=mc_response, + calibration_energies=mc_energies, +) +``` + +### Analytical Model + +Parameterized model with Gaussian core + exponential tail, suitable for quick estimates when no tabulated response is available. + +```python +detector = DetectorResponse.from_gaussian_params( + channel_energy_range=(1.0, 1.6), # in m_e units + n_channels=4096, + sigma_a=0.0, # energy-independent resolution (m_e) + sigma_b=0.0, # sqrt(E) contribution (m_e^{-1/2}) + tail_fraction=0.0, # fraction of events in exponential tail + tau=0.01, # decay constant of tail (m_e) + model="gaussian", # "gaussian" | "gaussian_tail" | "tikhonov" + fano_factor=0.0, # Fano factor (0.12 for Si, ~0 for MMC) +) +``` + +## Analytical Response Models + +### Gaussian Core + +Pure Gaussian peak with energy-dependent resolution: + +$$G(E; E_0, \\sigma) = \\frac{1}{\\sigma\\sqrt{2\\pi}} \\exp\\left(-\\frac{(E - E_0)^2}{2\\sigma^2}\\right)$$ + +### Gaussian + Exponential Tail + +Models low-energy tailing from charge collection inefficiencies in semiconductor detectors (standard for Si-PIPS detectors): + +$$f(E) = (1 - f\_{\\text{tail}}) \\cdot G(E; E_0, \\sigma) + f\_{\\text{tail}} \\cdot T(E; E_0, \\sigma, \\tau)$$ + +where $T$ is the Tikhonov function: + +$$T(E) = \\frac{1}{2\\tau} \\exp\\left(\\frac{\\tau(E - E_0)}{2} + \\frac{\\sigma^2\\tau^2}{8}\\right) \\cdot \\text{erfc}\\left(\\frac{E - E_0}{\\sqrt{2}\\sigma} + \\frac{\\sigma\\tau}{\\sqrt{2}}\\right)$$ + +### Energy-Dependent Resolution + +$$\\sigma(E) = \\sqrt{\\sigma_a^2 + (\\sigma_b \\sqrt{E})^2 + F \\cdot E \\cdot w}$$ + +- $\\sigma_a$ — energy-independent contribution (electronic noise, etc.) +- $\\sigma_b$ — $\\sqrt{E}$ term (statistical fluctuations in charge creation) +- $F$ — Fano factor ($\\approx 0.12$ for Si, $\\approx 0$ for MMCs) +- $w$ — mean ionization energy ($\\approx 3.6$ eV for Si, converted to $m_e$ units) + +## Calculator Implementation + +> [!tip] Code mapping +> Module: `beta_spectrum/components/detector_response.py` → `DetectorResponse` class +> +> - **Analytical models**: Gaussian core, Gaussian + exponential tail (Tikhonov) +> - **Tabulated mode**: Linear interpolation from MC-calculated response matrix +> - **Convolution**: `DetectorResponse.convolve(W, spectrum, normalize=True)` +> - **Integration with BetaSpectrum**: `BetaSpectrum.convolve_with_detector(detector, W, config)` +> - **Integration with BetaSpectrumAnalyzer**: `BetaSpectrumAnalyzer.convolved_spectrum(detector, normalize)` + +### API Usage + +```python +from beta_spectrum import ( + BetaSpectrum, SpectrumConfig, DetectorResponse, + BetaSpectrumAnalyzer, T_to_W +) + +# 1. Create theoretical spectrum for 99Tc +config = SpectrumConfig( + Z_parent=43, Z_daughter=44, A_number=99, + endpoint_MeV=0.294 # 99Tc Q-value +) +spectrum = BetaSpectrum.from_config(config) +W, kinetic_MeV = spectrum.get_energy_grid(config) + +# 2. Define detector response (sigma = 1.0 keV in m_e units) +sigma_1keV_me = 1.0 / 511.0 # convert keV to m_e units +W0 = T_to_W(0.294) +detector = DetectorResponse.from_gaussian_params( + channel_energy_range=(1.0, W0 + 0.05), + n_channels=512, + sigma_a=sigma_1keV_me, + sigma_b=0.0, + tail_fraction=0.0, + model="gaussian", + fano_factor=0.0, +) + +# 3. Convolve +convolved = spectrum.convolve_with_detector(detector, W=W, config=config) + +# 4. Or use the Analyzer +analyzer = BetaSpectrumAnalyzer(spectrum, config) +convolved_normalized = analyzer.convolved_spectrum(detector, normalize=True) +``` + +## Convolution Algorithm + +The convolution proceeds in three steps: + +1. **Build response matrix**: For each energy bin $W_i$, compute the detector response $R(E_j; W_i)$ across all channels $j$. Result: matrix of shape $(N\_{\\text{bins}}, N\_{\\text{channels}})$. + +1. **Weight by spectrum**: Multiply each response row by the theoretical spectrum value $\\Phi(W_i)$ and bin width $\\Delta W_i$. + +1. **Sum over energies**: For each channel $j$, sum contributions from all energy bins: + $$M_j = \\sum_i \\Phi(W_i) \\cdot R(E_j; W_i) \\cdot \\Delta W_i = (R^T \\cdot (\\Phi \\cdot \\Delta W))\_j$$ + +The bin width is computed as the half-bin average: $\\Delta W_i = \\frac{1}{2}(\\Delta W\_{i-1/2} + \\Delta W\_{i+1/2})$. + +> [!note] Normalization +> Each response function $R(E; W_i)$ is normalized to unit area before convolution. This ensures that convolving a mono-energetic line produces a response with the same total counts. The final convolved spectrum preserves the integral of the theoretical spectrum. + +## Energy Units + +All detector response computations use **total energy in units of $m_e c^2$** (natural units, $m_e = 1$): + +- $m_e c^2 = 510.998950$ keV +- Conversion: $W = 1 + T / (m_e c^2)$ where $T$ is kinetic energy +- A resolution of $\\sigma = 1.0$ keV corresponds to $\\sigma = 1.0 / 511.0 \\approx 0.00196$ in $m_e$ units + +> [!warning] Unit consistency +> The detector channel energy range and the spectrum energy grid must be in the same units (total energy in $m_e$ units). Use `T_to_W()` from `beta_spectrum.utils` to convert kinetic energy in MeV to total energy in $m_e$ units. + +## Physical Interpretation + +The convolution operation smears the sharp theoretical spectrum into a realistic measured spectrum: + +- **Peak broadening**: Finite energy resolution widens spectral features +- **Endpoint smearing**: The sharp cutoff at $W_0$ becomes a gradual fall-off +- **Tail effects**: Charge collection inefficiencies create low-energy tails (modeled by `gaussian_tail`) +- **Count conservation**: Total number of events is preserved (integral of convolved spectrum = integral of theoretical spectrum) + +## References + +- Hayen et al., RevModPhys **90**, 015008 (2017) — Comprehensive review of beta-decay theory +- Tikhonov, Soviet Physics Doklady **8**, 591 (1964) — Original formulation of the response function +- Sigmund, Phys. Rev. **130**, 1497 (1963) — Low-energy tailing in semiconductor detectors diff --git a/tests/test_detector_response.py b/tests/test_detector_response.py new file mode 100644 index 0000000..e225799 --- /dev/null +++ b/tests/test_detector_response.py @@ -0,0 +1,580 @@ +# tests/test_detector_response.py — Detector response unit tests + +""" +Why test this? +-------------- +The detector response function is the bridge between theoretical beta spectra +and measured experimental data. It models how a mono-energetic electron beam +appears in the detector's energy spectrum — typically a Gaussian peak with +possible low-energy tailing due to charge collection inefficiencies. + +We test: + + 1. **Analytical response shapes**: Gaussian core, Gaussian+tail, normalization + 2. **Tabulated response**: interpolation from MC-calculated response matrix + 3. **Convolution**: theoretical spectrum → predicted measured spectrum + 4. **Resolution model**: energy-dependent σ(E) behavior + 5. **Edge cases**: zero resolution, zero tail fraction, boundary energies +""" + +import numpy as np + +from beta_spectrum.components.detector_response import DetectorResponse + + +class TestAnalyticalGaussianResponse: + """Test pure Gaussian detector response.""" + + def test_gaussian_peak_position(self): + """Gaussian response should peak at the input energy.""" + channel_energies = np.linspace(0.0, 5.0, 4096) + resp = DetectorResponse(channel_energies, model="gaussian", sigma_a=0.005) + + energy = 2.0 + response = resp.get_response(energy) + peak_idx = np.argmax(response) + assert np.isclose( + channel_energies[peak_idx], energy, atol=0.01 + ), "Gaussian peak should be at input energy" + + def test_gaussian_normalization(self): + """Gaussian response should integrate to 1 (unit area).""" + channel_energies = np.linspace(0.0, 5.0, 8192) + resp = DetectorResponse(channel_energies, model="gaussian", sigma_a=0.005) + + energy = 2.0 + response = resp.get_response(energy) + integral = np.trapezoid(response, channel_energies) + assert np.isclose( + integral, 1.0, rtol=1e-2 + ), f"Gaussian should normalize to 1, got {integral}" + + def test_gaussian_sigma_dependence(self): + """Larger sigma_a should produce wider Gaussian.""" + channel_energies = np.linspace(0.0, 5.0, 4096) + resp_narrow = DetectorResponse( + channel_energies, model="gaussian", sigma_a=0.002 + ) + resp_wide = DetectorResponse(channel_energies, model="gaussian", sigma_a=0.02) + + energy = 2.0 + resp_n = resp_narrow.get_response(energy) + resp_w = resp_wide.get_response(energy) + + # FWHM: count points where response > half max + half_max_n = np.max(resp_n) / 2 + half_max_w = np.max(resp_w) / 2 + fwhm_n = np.sum(resp_n > half_max_n) + fwhm_w = np.sum(resp_w > half_max_w) + assert fwhm_w > fwhm_n, "Larger sigma should produce wider peak" + + def test_delta_limit(self): + """Zero sigma should produce delta function at peak channel.""" + channel_energies = np.linspace(0.0, 5.0, 4096) + resp = DetectorResponse( + channel_energies, model="gaussian", sigma_a=0.0, fano_factor=0.0 + ) + + energy = 2.0 + response = resp.get_response(energy) + assert ( + np.sum(response > 0) == 1 + ), "Delta function should have single non-zero element" + assert np.isclose(response[np.argmax(response)], 1.0) + + +class TestAnalyticalGaussianTailResponse: + """Test Gaussian + exponential tail detector response.""" + + def test_tail_increases_low_energy_counts(self): + """Adding tail should produce non-zero tail contribution below peak.""" + channel_energies = np.linspace(0.0, 5.0, 4096) + resp_with_tail = DetectorResponse( + channel_energies, + model="gaussian_tail", + sigma_a=0.005, + tail_fraction=0.3, + tau=0.01, + ) + + energy = 2.0 + response = resp_with_tail.get_response(energy) + + # Below peak: the tail contribution should be non-zero + below_peak = channel_energies < energy + # The response below peak comes from both the Gaussian tail (left side) + # and the exponential tail. Both should be non-zero. + assert ( + np.sum(response[below_peak]) > 0 + ), "Response below peak should be non-zero with tail enabled" + + # The response at the peak should be at or above nearby energies + peak_idx = np.argmax(response) + assert peak_idx > 0, "Peak should not be at first channel" + assert ( + response[peak_idx] >= response[peak_idx - 1] + ), "Response at peak should be >= response just below peak" + + def test_tail_fraction_effect(self): + """Higher tail_fraction should produce more tail contribution.""" + channel_energies = np.linspace(0.0, 5.0, 4096) + + resp_low = DetectorResponse( + channel_energies, + model="gaussian_tail", + sigma_a=0.005, + tail_fraction=0.1, + tau=0.01, + ) + resp_high = DetectorResponse( + channel_energies, + model="gaussian_tail", + sigma_a=0.005, + tail_fraction=0.5, + tau=0.01, + ) + resp_gaussian_only = DetectorResponse( + channel_energies, model="gaussian", sigma_a=0.005 + ) + + energy = 2.0 + resp_l = resp_low.get_response(energy) + resp_h = resp_high.get_response(energy) + resp_g = resp_gaussian_only.get_response(energy) + + below_peak = channel_energies < energy + + # The tail contribution is the difference between gaussian_tail + # and pure gaussian (scaled by tail_fraction). + # Higher tail_fraction should have larger tail contribution. + # Tail contribution = f_tail * T(E), where T is the tail function. + # Since T is only non-zero for E < E0, the tail contribution + # below peak is proportional to f_tail. + tail_low = resp_l[below_peak] - (1.0 - 0.1) * resp_g[below_peak] + tail_high = resp_h[below_peak] - (1.0 - 0.5) * resp_g[below_peak] + + # Both should have positive tail contribution + assert np.sum(tail_low) > 0, "Low tail_fraction should have positive tail" + assert np.sum(tail_high) > 0, "High tail_fraction should have positive tail" + + # Higher tail_fraction should have proportionally more tail contribution + # (tail_high / tail_low should be roughly 0.5/0.1 = 5) + ratio = np.sum(tail_high) / np.sum(tail_low) + assert ratio > 1.0, ( + f"Higher tail_fraction should have more tail contribution, " + f"ratio={ratio:.2f}" + ) + + def test_tail_normalization(self): + """Gaussian+tail response should still integrate to 1.""" + channel_energies = np.linspace(0.0, 5.0, 8192) + resp = DetectorResponse( + channel_energies, + model="gaussian_tail", + sigma_a=0.005, + tail_fraction=0.2, + tau=0.01, + ) + + energy = 2.0 + response = resp.get_response(energy) + integral = np.trapezoid(response, channel_energies) + assert np.isclose( + integral, 1.0, rtol=1e-2 + ), f"Gaussian+tail should normalize to 1, got {integral}" + + +class TestResolutionSigma: + """Test energy-dependent resolution sigma(E).""" + + def test_sigma_increases_with_energy(self): + """σ(E) should increase with energy for sigma_b > 0.""" + channel_energies = np.linspace(0.0, 5.0, 100) + resp = DetectorResponse( + channel_energies, model="gaussian", sigma_a=0.0, sigma_b=0.01 + ) + + sigma_low = resp._resolution_sigma(0.5) + sigma_high = resp._resolution_sigma(4.0) + assert sigma_high > sigma_low, "σ should increase with √E" + + def test_sigma_a_only(self): + """With sigma_b=0, σ should be constant.""" + channel_energies = np.linspace(0.0, 5.0, 100) + resp = DetectorResponse( + channel_energies, + model="gaussian", + sigma_a=0.01, + sigma_b=0.0, + fano_factor=0.0, + ) + + sigma_low = resp._resolution_sigma(0.5) + sigma_high = resp._resolution_sigma(4.0) + assert np.isclose(sigma_low, sigma_high), "σ should be constant when sigma_b=0" + + def test_fano_contribution(self): + """Fano factor should add energy-dependent contribution to σ.""" + channel_energies = np.linspace(0.0, 5.0, 100) + resp_no_fano = DetectorResponse( + channel_energies, + model="gaussian", + sigma_a=0.0, + sigma_b=0.0, + fano_factor=0.0, + ) + resp_fano = DetectorResponse( + channel_energies, + model="gaussian", + sigma_a=0.0, + sigma_b=0.0, + fano_factor=0.12, + ) + + sigma_low = resp_no_fano._resolution_sigma(0.5) + sigma_fano = resp_fano._resolution_sigma(0.5) + assert sigma_fano > sigma_low, "Fano factor should increase σ" + + +class TestTabulatedResponse: + """Test tabulated (MC-simulated) detector response.""" + + def test_tabulated_interpolation(self): + """Tabulated response should interpolate between calibration energies.""" + channel_energies = np.linspace(0.0, 5.0, 256) + calib_energies = np.array([1.0, 2.0, 3.0, 4.0]) + + # Simple diagonal response matrix: each calib energy peaks at its channel + n_channels = len(channel_energies) + response_matrix = np.zeros((len(calib_energies), n_channels)) + for i, e in enumerate(calib_energies): + idx = np.argmin(np.abs(channel_energies - e)) + response_matrix[i, idx] = 1.0 + + resp = DetectorResponse( + channel_energies=channel_energies, + response_matrix=response_matrix, + calibration_energies=calib_energies, + ) + + # Response at 2.5 should interpolate between 2.0 and 3.0 rows + response = resp.get_response(2.5) + assert np.any(response > 0), "Interpolated response should be non-zero" + + def test_tabulated_validation_wrong_columns(self): + """Should reject response matrix with wrong column count.""" + channel_energies = np.linspace(0.0, 5.0, 100) + calib_energies = np.array([1.0, 2.0]) + response_matrix = np.zeros((2, 200)) # 200 columns, but 100 channels + + try: + DetectorResponse( + channel_energies=channel_energies, + response_matrix=response_matrix, + calibration_energies=calib_energies, + ) + assert False, "Should have raised ValueError" + except ValueError: + pass + + def test_tabulated_outside_range(self): + """Response outside calibration range should be zero.""" + channel_energies = np.linspace(0.0, 5.0, 256) + calib_energies = np.array([1.0, 2.0, 3.0]) + + response_matrix = np.zeros((3, len(channel_energies))) + for i, e in enumerate(calib_energies): + idx = np.argmin(np.abs(channel_energies - e)) + response_matrix[i, idx] = 1.0 + + resp = DetectorResponse( + channel_energies=channel_energies, + response_matrix=response_matrix, + calibration_energies=calib_energies, + ) + + # Energy below calibration range + response = resp.get_response(0.5) + assert np.all(response == 0), "Response below calibration range should be zero" + + # Energy above calibration range + response = resp.get_response(4.0) + assert np.all(response == 0), "Response above calibration range should be zero" + + +class TestConvolution: + """Test spectrum convolution with detector response.""" + + def test_convolution_preserves_total_counts(self): + """Total counts should be approximately preserved after convolution.""" + channel_energies = np.linspace(0.0, 5.0, 4096) + resp = DetectorResponse(channel_energies, model="gaussian", sigma_a=0.005) + + # Mono-energetic line at W=2.5 + W = np.linspace(0.0, 5.0, 1000) + spectrum = np.zeros_like(W) + idx = np.argmin(np.abs(W - 2.5)) + spectrum[idx] = 1.0 / (W[1] - W[0]) # Delta function approximation + + convolved = resp.convolve(W, spectrum, normalize=True) + total_theoretical = np.trapezoid(spectrum, W) + total_convolved = np.trapezoid(convolved, channel_energies) + assert np.isclose( + total_convolved, total_theoretical, rtol=0.05 + ), f"Total counts: theory={total_theoretical:.4f}, convolved={total_convolved:.4f}" + + def test_convolution_widens_peak(self): + """Convolution should widen a narrow peak.""" + channel_energies = np.linspace(0.0, 5.0, 4096) + resp = DetectorResponse(channel_energies, model="gaussian", sigma_a=0.02) + + # Very narrow peak + W = np.linspace(0.0, 5.0, 2000) + spectrum = np.zeros_like(W) + idx = np.argmin(np.abs(W - 2.5)) + spectrum[idx] = 1.0 / (W[1] - W[0]) + + convolved = resp.convolve(W, spectrum, normalize=True) + peak_width = np.sum(convolved > np.max(convolved) / 2) + assert peak_width > 1, "Convolved peak should have finite width" + + def test_convolution_beta_spectrum(self): + """Convolution of a beta spectrum should produce a smooth curve.""" + channel_energies = np.linspace(0.0, 5.0, 4096) + resp = DetectorResponse(channel_energies, model="gaussian", sigma_a=0.01) + + # Simple phase space spectrum ∝ p·W·(W0-W)² + W0 = 5.0 + W = np.linspace(1.0, W0 - 0.01, 500) + p = np.sqrt(W**2 - 1) + spectrum = p * W * (W0 - W) ** 2 + + convolved = resp.convolve(W, spectrum, normalize=True) + assert len(convolved) == len(channel_energies) + assert np.all(convolved >= 0), "Convolved spectrum should be non-negative" + assert np.max(convolved) > 0, "Convolved spectrum should have non-zero peak" + + def test_convolve_batch_same_as_convolve(self): + """Batch convolution should match standard convolution for analytical models.""" + channel_energies = np.linspace(0.0, 5.0, 4096) + resp = DetectorResponse( + channel_energies, + model="gaussian_tail", + sigma_a=0.01, + tail_fraction=0.1, + tau=0.01, + ) + + W = np.linspace(1.0, 4.9, 300) + p = np.sqrt(W**2 - 1) + spectrum = p * W * (4.9 - W) ** 2 + + result_standard = resp.convolve(W, spectrum, normalize=True) + result_batch = resp.convolve_batch(W, spectrum, normalize=True) + + assert np.allclose( + result_standard, result_batch, rtol=1e-5 + ), "Batch and standard convolution should match for analytical models" + + +class TestFactoryMethods: + """Test convenience factory methods.""" + + def test_from_gaussian_params(self): + """from_gaussian_params should create correct DetectorResponse.""" + resp = DetectorResponse.from_gaussian_params( + channel_energy_range=(0.0, 5.0), + n_channels=1024, + sigma_a=0.005, + sigma_b=0.01, + tail_fraction=0.2, + tau=0.01, + model="gaussian_tail", + fano_factor=0.12, + ) + + assert resp.n_channels == 1024 + assert np.isclose(resp.channel_energies[0], 0.0) + assert np.isclose(resp.channel_energies[-1], 5.0) + assert resp._mode == "analytical" + assert resp.tail_fraction == 0.2 + + def test_from_mc_simulation(self): + """from_mc_simulation should create tabulated DetectorResponse.""" + channel_energies = np.linspace(0.0, 5.0, 512) + calib_energies = np.array([1.0, 2.0, 3.0]) + response_matrix = np.zeros((3, 512)) + for i, e in enumerate(calib_energies): + idx = np.argmin(np.abs(channel_energies - e)) + response_matrix[i, idx] = 1.0 + + resp = DetectorResponse.from_mc_simulation( + channel_energies=channel_energies, + response_matrix=response_matrix, + calibration_energies=calib_energies, + ) + + assert resp._mode == "tabulated" + assert resp.n_channels == 512 + assert len(resp.calibration_energies) == 3 + + +class TestEdgeCases: + """Test edge cases and error handling.""" + + def test_repr_analytical(self): + """__repr__ should show analytical mode info.""" + channel_energies = np.linspace(0.0, 5.0, 100) + resp = DetectorResponse(channel_energies, model="gaussian", sigma_a=0.01) + repr_str = repr(resp) + assert "analytical" in repr_str + assert "gaussian" in repr_str + + def test_repr_tabulated(self): + """__repr__ should show tabulated mode info.""" + channel_energies = np.linspace(0.0, 5.0, 100) + calib_energies = np.array([1.0, 2.0]) + response_matrix = np.zeros((2, 100)) + resp = DetectorResponse( + channel_energies=channel_energies, + response_matrix=response_matrix, + calibration_energies=calib_energies, + ) + repr_str = repr(resp) + assert "tabulated" in repr_str + + def test_invalid_model(self): + """Should reject unknown model names.""" + channel_energies = np.linspace(0.0, 5.0, 100) + resp = DetectorResponse(channel_energies, model="unknown_model") + try: + resp.get_response(2.0) + assert False, "Should have raised ValueError" + except ValueError as e: + assert "Unknown model" in str(e) + + def test_convolve_mismatched_lengths(self): + """Should reject mismatched W and spectrum_values lengths.""" + channel_energies = np.linspace(0.0, 5.0, 100) + resp = DetectorResponse(channel_energies, model="gaussian") + W = np.linspace(0.0, 5.0, 100) + spectrum = np.zeros(50) # Wrong length + try: + resp.convolve(W, spectrum) + assert False, "Should have raised ValueError" + except ValueError: + pass + + def test_non_1d_input(self): + """Should reject non-1D inputs.""" + channel_energies = np.linspace(0.0, 5.0, 100) + resp = DetectorResponse(channel_energies, model="gaussian") + W = np.array([[0.0, 5.0]]) # 2D + spectrum = np.array([[1.0, 1.0]]) + try: + resp.convolve(W, spectrum) + assert False, "Should have raised ValueError" + except ValueError: + pass + + +class TestBetaSpectrumIntegration: + """Integration tests between BetaSpectrum and DetectorResponse.""" + + def test_beta_spectrum_convolve_with_detector(self): + """BetaSpectrum should convolve with detector response.""" + from beta_spectrum import BetaSpectrum, SpectrumConfig, DetectorResponse + from beta_spectrum.utils import T_to_W + + config = SpectrumConfig( + Z_parent=43, Z_daughter=44, A_number=99, endpoint_MeV=0.294 + ) + spectrum = BetaSpectrum.from_config(config) + W, _ = spectrum.get_energy_grid(config) + + # Detector channel energies must be in m_e units (total energy) + W0 = T_to_W(0.294) + detector = DetectorResponse.from_gaussian_params( + channel_energy_range=(1.0, W0 + 0.05), + n_channels=512, + sigma_a=0.003, + sigma_b=0.0, + tail_fraction=0.0, + model="gaussian", + fano_factor=0.0, + ) + + convolved = spectrum.convolve_with_detector(detector, W=W, config=config) + + assert convolved.shape == (512,) + assert np.all(convolved >= 0), "Convolved spectrum should be non-negative" + assert np.sum(convolved) > 0, "Convolved spectrum should have non-zero area" + + def test_convolved_spectrum_has_wider_peak(self): + """Convolved spectrum should be wider than theoretical spectrum.""" + from beta_spectrum import BetaSpectrum, SpectrumConfig, DetectorResponse + from beta_spectrum.utils import T_to_W + + config = SpectrumConfig( + Z_parent=43, Z_daughter=44, A_number=99, endpoint_MeV=0.294 + ) + spectrum = BetaSpectrum.from_config(config) + W, _ = spectrum.get_energy_grid(config) + + # Theoretical spectrum + theoretical = spectrum(W) + + # Detector response with finite resolution + W0 = T_to_W(0.294) + detector = DetectorResponse.from_gaussian_params( + channel_energy_range=(1.0, W0 + 0.05), + n_channels=512, + sigma_a=0.005, + sigma_b=0.0, + tail_fraction=0.0, + model="gaussian", + fano_factor=0.0, + ) + + convolved = spectrum.convolve_with_detector(detector, W=W, config=config) + + # The convolved spectrum should be smoother/wider + theoretical_width = np.sum(theoretical > np.max(theoretical) * 0.5) + convolved_width = np.sum(convolved > np.max(convolved) * 0.5) + assert ( + convolved_width >= theoretical_width + ), "Convolved spectrum should be at least as wide as theoretical" + + def test_analyzer_convolved_spectrum(self): + """BetaSpectrumAnalyzer should compute convolved spectrum.""" + from beta_spectrum import ( + BetaSpectrum, + SpectrumConfig, + DetectorResponse, + BetaSpectrumAnalyzer, + ) + from beta_spectrum.utils import T_to_W + + config = SpectrumConfig( + Z_parent=43, Z_daughter=44, A_number=99, endpoint_MeV=0.294 + ) + spectrum = BetaSpectrum.from_config(config) + analyzer = BetaSpectrumAnalyzer(spectrum, config) + + W0 = T_to_W(0.294) + detector = DetectorResponse.from_gaussian_params( + channel_energy_range=(1.0, W0 + 0.05), + n_channels=256, + sigma_a=0.003, + sigma_b=0.0, + tail_fraction=0.0, + model="gaussian", + fano_factor=0.0, + ) + + convolved = analyzer.convolved_spectrum(detector, normalize=True) + + assert convolved.shape == (256,) + assert np.isclose( + np.trapezoid(convolved, detector.channel_energies), 1.0, rtol=0.01 + ), "Normalized convolved spectrum should integrate to 1"