""" CsI(Tl) detector response continuum calibration for Radiacode 103. Models ONLY the detector's noise continuum. Photopeaks from environmental isotopes depend on measurement location and are NOT part of this model. Uses iterative peak subtraction followed by Gaussian smoothing to produce a clean continuum shape. This approach tracks the measured background closely at all energies, unlike log-space splines which collapse in low-signal regions. """ import json import numpy as np from pathlib import Path from scipy.ndimage import gaussian_filter1d from scipy.signal import savgol_filter from app.config import ENERGY_OFFSET, ENERGY_SLOPE, NUM_CHANNELS PHOTOPEAK_LINES = [ (295.22, 0.1842), (351.93, 0.3560), (609.31, 0.4549), (911.20, 0.2580), (968.97, 0.1580), (1120.29, 0.1492), (1460.83, 0.1066), (1764.49, 0.1531), (2614.51, 0.3586), ] def _sigma_keV(E): return max(12.0, 23.6 * np.sqrt(max(E, 1.0) / 662.0)) def _sigma_ch(E_keV): fwhm_keV = 0.08 * E_keV * (E_keV / 662.0) ** 0.5 sigma_keV = fwhm_keV / 2.355 return max(sigma_keV / ENERGY_SLOPE, 2.0) def _subtract_peaks(energy_axis, spectrum): """Remove known isotope photopeaks from spectrum.""" continuum = spectrum.copy() channels = np.arange(len(spectrum), dtype=np.float64) for line_energy, _ in PHOTOPEAK_LINES: idx = int(np.argmin(np.abs(energy_axis - line_energy))) sig = _sigma_ch(line_energy) far = int(5 * sig) + 3 lo_start = max(0, idx - far - int(3 * sig)) lo_end = max(0, idx - far) hi_start = min(len(spectrum), idx + far) hi_end = min(len(spectrum), idx + far + int(3 * sig)) baseline_regions = [] if lo_end > lo_start: baseline_regions.append(continuum[lo_start:lo_end]) if hi_end > hi_start: baseline_regions.append(continuum[hi_start:hi_end]) if not baseline_regions: continue local_bg = float(np.median(np.concatenate(baseline_regions))) peak_height = continuum[idx] - local_bg if peak_height > 0: gaussian = peak_height * np.exp(-0.5 * ((channels - idx) / sig) ** 2) continuum -= gaussian return np.maximum(continuum, 0), [{"energy_keV": e, "amplitude": 0.0} for e, _ in PHOTOPEAK_LINES] def calibrate_spline(measured_cps, energy_axis): """ Fit continuum using peak subtraction + Gaussian smoothing. Uses scipy's gaussian_filter1d after iterative peak subtraction, producing a smooth continuum that tracks the measured background closely at all energies including the high-energy tail. """ # Step 1: Smooth to reduce statistical noise window = min(51, len(measured_cps) // 10 * 2 + 1) if window < 5: window = 5 y_smooth = savgol_filter(measured_cps, window_length=window, polyorder=3) # Step 2: Subtract known photopeaks continuum, peak_amplitudes = _subtract_peaks(energy_axis, y_smooth) # Step 3: Gaussian smooth for final continuum shape sigma = max(15, len(continuum) // 60) continuum_smooth = gaussian_filter1d(continuum, sigma=sigma) continuum_smooth = np.maximum(continuum_smooth, 0) # Quality metrics residuals = continuum - continuum_smooth ss_res = np.sum(residuals ** 2) ss_tot = np.sum((continuum - continuum.mean()) ** 2) r_squared = 1.0 - ss_res / ss_tot if ss_tot > 0 else 0 return { "continuum_cps": continuum_smooth, "peak_amplitudes": peak_amplitudes, "r_squared": float(r_squared), "residuals_rms": float(np.sqrt(np.mean(residuals ** 2))), } def calibrate_background(measured_cps, energy_axis): """Fit the continuum model using peak subtraction + Gaussian smoothing.""" result = calibrate_spline(measured_cps, energy_axis) if "error" in result: return result return { "params": {}, "continuum_cps": result["continuum_cps"], "peak_amplitudes": result["peak_amplitudes"], "r_squared": result["r_squared"], "residuals_rms": result["residuals_rms"], "method": "peak_subtract_gaussian", } def build_calibrated_continuum(energy_axis, total_counts, params): """Build the continuum from calibrated parameters.""" if "continuum_cps" in params: cps = np.array(params["continuum_cps"]) if cps.sum() > 0: return cps * total_counts / cps.sum() return cps return np.zeros(len(energy_axis)) # Cached calibration _cached_result = None _CALIBRATION_PATH = Path("/data/bg_calibration.json") def load_or_calibrate(): """Load cached calibration or fit from measured data.""" global _cached_result if _cached_result is not None: return _cached_result # Try loading from cache file first (read-only volume is fine for reads) if _CALIBRATION_PATH.exists(): try: with open(_CALIBRATION_PATH) as f: cached = json.load(f) # Invalidate if method changed if cached.get("method") != "peak_subtract_gaussian": cached = None except Exception: cached = None if cached and "continuum_cps" in cached: _cached_result = cached return _cached_result from app.config import BACKGROUND_PATH, BACKGROUND_SNAPSHOT_PATH measured_counts = None live_time = 0 if BACKGROUND_PATH.exists(): try: bg_data = np.load(str(BACKGROUND_PATH), allow_pickle=True).item() measured_counts = bg_data["counts"].astype(np.float64)[:NUM_CHANNELS] live_time = float(bg_data["duration"]) except Exception: pass if measured_counts is None and BACKGROUND_SNAPSHOT_PATH.exists(): try: with open(BACKGROUND_SNAPSHOT_PATH) as f: snapshot = json.load(f) measured_counts = np.array(snapshot.get("spectrum", [])[:NUM_CHANNELS], dtype=np.float64) live_time = float(snapshot.get("live_time_s", 0)) except Exception: pass if measured_counts is None or live_time < 600: return None channels = np.arange(NUM_CHANNELS, dtype=np.float64) e_axis = ENERGY_OFFSET + ENERGY_SLOPE * channels measured_cps = measured_counts / live_time result = calibrate_background(measured_cps, e_axis) if "error" in result: return None _cached_result = { "continuum_cps": [round(float(c), 6) for c in result["continuum_cps"]], "method": result["method"], "r_squared": result["r_squared"], } # Write cache if volume is writable (may fail on read-only mounts) try: _CALIBRATION_PATH.parent.mkdir(parents=True, exist_ok=True) tmp = _CALIBRATION_PATH.with_suffix(".tmp") with open(tmp, "w") as f: json.dump(_cached_result, f, indent=2) tmp.replace(_CALIBRATION_PATH) except OSError: pass # Read-only volume — in-memory cache is sufficient return _cached_result