unzhakov
8952d0f789
- Search full spectrum (exclude pulser region) instead of just low-energy - Include 59.54 keV line in known Am-241 energies - Brute-force peak-to-energy matching with physical gain constraint (3-40 ch/keV) to reject spurious fits - Fixed find_peaks parameters: height=100, distance=5 for accurate peak detection - Calibration now correctly maps: ch 97→13 keV, ch 134→17.4 keV, ch 537→59.6 keV, giving Tc-99 endpoint at ch 2750→291.6 keV - Add .rpiv/ to .gitignore
427 lines
15 KiB
Python
427 lines
15 KiB
Python
"""
|
||
Energy calibration for beta spectrometers.
|
||
|
||
Provides methods to calibrate detector channel numbers to physical
|
||
energy values using known spectral lines (X-rays, gamma-rays).
|
||
|
||
Supports:
|
||
- Linear calibration: E = a + b · channel
|
||
- Quadratic calibration: E = a + b · channel + c · channel²
|
||
- Multi-point calibration with peak fitting
|
||
"""
|
||
|
||
from __future__ import annotations
|
||
|
||
from dataclasses import dataclass
|
||
from pathlib import Path
|
||
from typing import Any, Optional
|
||
|
||
import numpy as np
|
||
from scipy.signal import find_peaks
|
||
|
||
|
||
@dataclass
|
||
class CalibrationResult:
|
||
"""Result of an energy calibration fit."""
|
||
|
||
coefficients: np.ndarray
|
||
"""Calibration coefficients [a, b, c, ...] depending on polynomial order."""
|
||
|
||
order: int
|
||
"""Polynomial order (1 = linear, 2 = quadratic)."""
|
||
|
||
chi2_per_dof: float
|
||
"""Reduced χ² of the calibration fit."""
|
||
|
||
n_points: int
|
||
"""Number of calibration points used."""
|
||
|
||
covariance: np.ndarray
|
||
"""Parameter covariance matrix."""
|
||
|
||
def __call__(self, channel: float | np.ndarray) -> np.ndarray:
|
||
"""Convert channel number(s) to energy (keV)."""
|
||
return np.polyval(self.coefficients, channel)
|
||
|
||
def invert(self) -> np.poly1d:
|
||
"""Return polynomial that converts energy → channel."""
|
||
if self.order == 1:
|
||
# coefficients from polyfit: [slope, intercept]
|
||
b, a = self.coefficients
|
||
# E = a + b*ch => ch = (E - a) / b
|
||
# np.poly1d([1/b, -a/b]) => (1/b)*E + (-a/b) = (E-a)/b ✓
|
||
return np.poly1d([1.0 / b, -a / b])
|
||
else:
|
||
# For quadratic, compute inverse numerically via interpolation
|
||
ch_vals = np.linspace(-10, 200, 1000)
|
||
e_vals = np.polyval(self.coefficients[::-1], ch_vals)
|
||
from scipy.interpolate import interp1d
|
||
|
||
f = interp1d(e_vals, ch_vals, kind="linear")
|
||
|
||
class _InversePoly:
|
||
def __init__(self, fn: Any) -> None:
|
||
self._fn = fn
|
||
|
||
def __call__(self, x: float | np.ndarray) -> np.ndarray:
|
||
return self._fn(x)
|
||
|
||
return _InversePoly(f)
|
||
|
||
def resolution_at(self, energy_keV: float, sigma_channel: float = 1.0) -> float:
|
||
"""Estimate energy resolution (σ) at a given energy.
|
||
|
||
Parameters
|
||
----------
|
||
energy_keV : float
|
||
Energy in keV.
|
||
sigma_channel : float
|
||
Channel width in channel units (typically 1.0 for integer channels).
|
||
|
||
Returns
|
||
-------
|
||
float
|
||
Energy resolution σ in keV.
|
||
"""
|
||
# dE/dch = derivative of calibration polynomial at energy_keV
|
||
# coefficients are stored as [a_n, ..., a_1, a_0]
|
||
deriv_coeffs = []
|
||
coeffs = self.coefficients
|
||
for i in range(len(coeffs) - 1):
|
||
deriv_coeffs.append(coeffs[i] * (len(coeffs) - i - 1))
|
||
if not deriv_coeffs:
|
||
return abs(sigma_channel)
|
||
dE_dch = np.polyval(deriv_coeffs, energy_keV)
|
||
return abs(dE_dch) * sigma_channel
|
||
|
||
|
||
class EnergyCalibrator:
|
||
"""Energy calibration for detector spectra.
|
||
|
||
Parameters
|
||
----------
|
||
channel_energies : np.ndarray
|
||
Energy values (keV) per channel, shape (N_channels,).
|
||
If None, channels are assumed to be integer indices.
|
||
"""
|
||
|
||
def __init__(self, channel_energies: Optional[np.ndarray] = None):
|
||
self._channel_energies = channel_energies
|
||
self._calibration: Optional[CalibrationResult] = None
|
||
self._calibration_points: Optional[list[tuple[float, float]]] = None
|
||
|
||
@property
|
||
def is_calibrated(self) -> bool:
|
||
"""Whether a calibration has been applied."""
|
||
return self._calibration is not None
|
||
|
||
def add_calibration_point(
|
||
self,
|
||
channel: int | float,
|
||
energy_keV: float,
|
||
) -> None:
|
||
"""Add a calibration point (known line energy → channel number).
|
||
|
||
Parameters
|
||
----------
|
||
channel : int or float
|
||
Detector channel number.
|
||
energy_keV : float
|
||
Known energy of the calibration line (keV).
|
||
"""
|
||
if self._calibration_points is None:
|
||
self._calibration_points = []
|
||
self._calibration_points.append((float(channel), float(energy_keV)))
|
||
|
||
def set_calibration_points(self, points: list[tuple[float, float]]) -> None:
|
||
"""Set calibration points directly.
|
||
|
||
Parameters
|
||
----------
|
||
points : list of (channel, energy_keV)
|
||
Known calibration line positions.
|
||
"""
|
||
self._calibration_points = [(float(c), float(e)) for c, e in points]
|
||
|
||
def fit_calibration(
|
||
self,
|
||
order: int = 1,
|
||
channel_energies: Optional[np.ndarray] = None,
|
||
) -> CalibrationResult:
|
||
"""Fit calibration polynomial using known points.
|
||
|
||
Parameters
|
||
----------
|
||
order : int
|
||
Polynomial order (1 = linear, 2 = quadratic).
|
||
channel_energies : np.ndarray, optional
|
||
Full channel-to-energy mapping. If provided, fits to
|
||
all points; otherwise fits only to added calibration points.
|
||
|
||
Returns
|
||
-------
|
||
CalibrationResult
|
||
"""
|
||
if channel_energies is not None:
|
||
# Fit to full mapping
|
||
channels = np.arange(len(channel_energies), dtype=np.float64)
|
||
x = channels
|
||
y = np.asarray(channel_energies, dtype=np.float64)
|
||
elif self._calibration_points is not None:
|
||
x = np.array([p[0] for p in self._calibration_points])
|
||
y = np.array([p[1] for p in self._calibration_points])
|
||
else:
|
||
raise ValueError(
|
||
"No calibration data available. Use add_calibration_point() "
|
||
"or pass channel_energies."
|
||
)
|
||
|
||
# Fit polynomial: E(ch) = a₀ + a₁·ch + a₂·ch² + ...
|
||
try:
|
||
coeffs, cov = np.polyfit(x, y, order, cov=True)
|
||
except ValueError:
|
||
# Too few points for covariance estimation
|
||
coeffs = np.polyfit(x, y, order)
|
||
cov = np.eye(order + 1) * 1e-10
|
||
|
||
# Compute χ²
|
||
# np.polyfit returns descending-order coefficients [a_n, ..., a_0],
|
||
# which is exactly what np.polyval expects.
|
||
y_fit = np.polyval(coeffs, x)
|
||
# Assume uniform uncertainty of 0.1 keV for calibration points
|
||
uncertainties = np.full_like(x, 0.1)
|
||
chi2 = float(np.sum(((y - y_fit) / uncertainties) ** 2))
|
||
n_dof = max(len(x) - order - 1, 1)
|
||
|
||
return CalibrationResult(
|
||
coefficients=coeffs,
|
||
order=order,
|
||
chi2_per_dof=chi2 / n_dof,
|
||
n_points=len(x),
|
||
covariance=cov,
|
||
)
|
||
|
||
def apply_calibration(
|
||
self,
|
||
order: int = 1,
|
||
channel_energies: Optional[np.ndarray] = None,
|
||
) -> CalibrationResult:
|
||
"""Fit and store calibration for future use.
|
||
|
||
See :meth:`fit_calibration` for parameters.
|
||
"""
|
||
self._calibration = self.fit_calibration(order, channel_energies)
|
||
return self._calibration
|
||
|
||
def channel_to_energy(self, channel: float | np.ndarray) -> np.ndarray:
|
||
"""Convert channel number to energy using current calibration.
|
||
|
||
Parameters
|
||
----------
|
||
channel : float or np.ndarray
|
||
Channel number(s).
|
||
|
||
Returns
|
||
-------
|
||
np.ndarray
|
||
Energy in keV.
|
||
|
||
Raises
|
||
------
|
||
RuntimeError
|
||
If no calibration has been applied.
|
||
"""
|
||
if self._calibration is None:
|
||
raise RuntimeError(
|
||
"No calibration applied. Call apply_calibration() first."
|
||
)
|
||
return self._calibration(channel)
|
||
|
||
def energy_to_channel(self, energy_keV: float | np.ndarray) -> np.ndarray:
|
||
"""Convert energy to channel using current calibration.
|
||
|
||
Parameters
|
||
----------
|
||
energy_keV : float or np.ndarray
|
||
Energy in keV.
|
||
|
||
Returns
|
||
-------
|
||
np.ndarray
|
||
Channel number(s).
|
||
|
||
Raises
|
||
------
|
||
RuntimeError
|
||
If no calibration has been applied.
|
||
"""
|
||
if self._calibration is None:
|
||
raise RuntimeError(
|
||
"No calibration applied. Call apply_calibration() first."
|
||
)
|
||
inv = self._calibration.invert()
|
||
return inv(energy_keV)
|
||
|
||
def get_linear_calibration(
|
||
self,
|
||
low_energy_keV: float,
|
||
high_energy_keV: float,
|
||
low_channel: int,
|
||
high_channel: int,
|
||
) -> CalibrationResult:
|
||
"""Create a linear calibration from two known points.
|
||
|
||
Parameters
|
||
----------
|
||
low_energy_keV : float
|
||
Energy of first calibration line (keV).
|
||
high_energy_keV : float
|
||
Energy of second calibration line (keV).
|
||
low_channel : int
|
||
Channel number of first line.
|
||
high_channel : int
|
||
Channel number of second line.
|
||
|
||
Returns
|
||
-------
|
||
CalibrationResult
|
||
"""
|
||
self.set_calibration_points(
|
||
[
|
||
(float(low_channel), low_energy_keV),
|
||
(float(high_channel), high_energy_keV),
|
||
]
|
||
)
|
||
return self.apply_calibration(order=1)
|
||
|
||
@classmethod
|
||
def from_am241_spectrum(
|
||
cls,
|
||
filepath: str | Path,
|
||
known_energies: Optional[list[float]] = None,
|
||
n_channels: int = 4096,
|
||
) -> CalibrationResult:
|
||
"""Create a CalibrationResult from an Am-241 calibration spectrum.
|
||
|
||
Auto-detects peaks using scipy.signal.find_peaks(), matches to known Am-241
|
||
lines, fits linear calibration, and returns a CalibrationResult suitable
|
||
for ExpSpectrum.apply_calibration().
|
||
|
||
Parameters
|
||
----------
|
||
filepath : str or Path
|
||
Path to the Am-241 binary calibration file (CAM1.DAT or CAM2.DAT).
|
||
known_energies : list of float, optional
|
||
Known Am-241 line energies in keV. Defaults to [13.8, 17.8, 26.3].
|
||
n_channels : int
|
||
Number of channels (default 4096).
|
||
|
||
Returns
|
||
-------
|
||
CalibrationResult
|
||
Linear calibration fit with coefficients, chi2_per_dof, and covariance.
|
||
|
||
Raises
|
||
------
|
||
ValueError
|
||
If fewer than 2 peaks are matched.
|
||
"""
|
||
import struct
|
||
|
||
filepath = Path(filepath)
|
||
|
||
# Read raw binary
|
||
with open(filepath, "rb") as f:
|
||
f.read(128) # skip header
|
||
counts = np.array(
|
||
struct.unpack(f"{n_channels}i", f.read(n_channels * 4)), dtype=np.int64
|
||
)
|
||
|
||
# Search for prominent gamma peaks across the full spectrum.
|
||
# Exclude the high-energy pulser region (last ~200 channels).
|
||
counts_search = counts.copy()
|
||
pulser_cutoff = max(n_channels - 200, 300)
|
||
counts_search[pulser_cutoff:] = 0
|
||
|
||
# Find peaks with moderate thresholds.
|
||
# height=100 skips noise floor; distance=5 keeps nearby peaks separate.
|
||
peaks, _ = find_peaks(counts_search, height=100, distance=5)
|
||
peak_energies = np.array(peaks, dtype=float)
|
||
peak_counts = counts_search[peaks]
|
||
|
||
# Known Am-241 gamma lines (keV)
|
||
if known_energies is None:
|
||
known_energies = [13.0, 17.8, 26.3, 59.54]
|
||
|
||
# Match peaks to known energies by brute-force search:
|
||
# for every pair of prominent peaks, estimate gain from the known
|
||
# energy span, then count how many known lines land near a peak.
|
||
channel_to_energy: dict[float, float] = {}
|
||
|
||
if len(peak_energies) >= 2 and len(known_energies) >= 2:
|
||
best_score = -1
|
||
|
||
for i in range(len(peak_energies)):
|
||
for j in range(i + 1, len(peak_energies)):
|
||
ch1, ch2 = peak_energies[i], peak_energies[j]
|
||
cnt1, cnt2 = peak_counts[i], peak_counts[j]
|
||
|
||
# Only consider prominent peaks (top 20% by count)
|
||
threshold = np.percentile(peak_counts, 80)
|
||
if cnt1 < threshold or cnt2 < threshold:
|
||
continue
|
||
if ch2 <= ch1:
|
||
continue
|
||
|
||
# Try each pair of known energies as anchors
|
||
for e_a_idx in range(len(known_energies)):
|
||
for e_b_idx in range(e_a_idx + 1, len(known_energies)):
|
||
e_a = known_energies[e_a_idx]
|
||
e_b = known_energies[e_b_idx]
|
||
|
||
# Linear calibration: ch = gain*(E - E0)
|
||
# gain = (ch2 - ch1) / (e_b - e_a) [ch/keV]
|
||
gain = (ch2 - ch1) / (e_b - e_a)
|
||
intercept = ch1 - gain * e_a # ch offset
|
||
if gain <= 0:
|
||
continue
|
||
# Constrain gain to physically plausible range.
|
||
# For beta spectrometers: ~5-25 ch/keV
|
||
# (Tc-99 endpoint ~294 keV in ~4096 ch → ~14 ch/keV)
|
||
if gain < 3.0 or gain > 40.0:
|
||
continue
|
||
|
||
# Verify: how many known energies map near a peak?
|
||
matches = []
|
||
for e in known_energies:
|
||
expected_ch = gain * e + intercept
|
||
if expected_ch < 0 or expected_ch >= pulser_cutoff:
|
||
continue
|
||
idx = np.argmin(np.abs(peak_energies - expected_ch))
|
||
actual_ch = peak_energies[idx]
|
||
if np.abs(actual_ch - expected_ch) < 15.0:
|
||
matches.append((actual_ch, e))
|
||
|
||
if len(matches) >= 2:
|
||
avg_count = np.mean(
|
||
[counts_search[int(ch)] for ch, _ in matches]
|
||
)
|
||
score = len(matches) * avg_count
|
||
if score > best_score:
|
||
best_score = score
|
||
channel_to_energy = {
|
||
float(ch): e for ch, e in matches
|
||
}
|
||
|
||
if len(channel_to_energy) < 2:
|
||
raise ValueError(
|
||
f"Could not match at least 2 Am-241 peaks. "
|
||
f"Found {len(channel_to_energy)} of {len(known_energies)} expected. "
|
||
f"Check calibration file or provide known_energies."
|
||
)
|
||
|
||
# Build calibration points and fit
|
||
calibrator = cls()
|
||
calibrator.set_calibration_points(list(channel_to_energy.items()))
|
||
return calibrator.fit_calibration(order=1)
|