unzhakov
5c2b12c0e6
- Fix chi2 calculation bug in calibration.py (remove incorrect coeffs[::-1] reversal in np.polyval) - Fix broken beta_spectrum/cli.py (remove orphaned detector response add_argument calls and unused imports) - Add ExpSpectrum.apply_calibration() and EnergyCalibrator.from_am241_spectrum() - Add LivetimeDetermination class for pulser-based livetime - Add SpectrumPipeline orchestration (read→subtract→calibrate→livetime→validate) - Add ROOT I/O with PyROOT primary + uproot fallback - Add bs_exp CLI with process, livetime, export subcommands - Add exp_calibration notebook with pre-calibration visualization - Add comprehensive tests (73 passing) - Add bs_exp console script and [root] extra to pyproject.toml - Remove obsolete detector_response module and unused notebooks
333 lines
12 KiB
Python
333 lines
12 KiB
Python
"""
|
||
Physics tests for DetectorResponse — unique physical behavior.
|
||
|
||
Tests detector response shapes, resolution, and convolution physics.
|
||
"""
|
||
|
||
import numpy as np
|
||
import pytest
|
||
|
||
try:
|
||
from beta_spectrum.components.detector_response import DetectorResponse
|
||
except ImportError:
|
||
DetectorResponse = None # type: ignore[misc, assignment]
|
||
|
||
|
||
@pytest.mark.skipif(DetectorResponse is None, reason="DetectorResponse not yet implemented")
|
||
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)
|
||
|
||
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 = channel_energies < energy
|
||
assert (
|
||
np.sum(response[below_peak]) > 0
|
||
), "Response below peak should be non-zero with tail enabled"
|
||
|
||
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
|
||
|
||
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]
|
||
|
||
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"
|
||
|
||
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])
|
||
|
||
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 = resp.get_response(2.5)
|
||
assert np.any(response > 0), "Interpolated response should be non-zero"
|
||
|
||
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,
|
||
)
|
||
|
||
response = resp.get_response(0.5)
|
||
assert np.all(response == 0), "Response below calibration range should be zero"
|
||
|
||
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)
|
||
|
||
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])
|
||
|
||
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)
|
||
|
||
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)
|
||
|
||
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"
|