radiacode/web/app/bg_calibration.py

"""
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