diff --git a/detect/capture_background.py b/detect/capture_background.py index e21bfff..70061b6 100644 --- a/detect/capture_background.py +++ b/detect/capture_background.py @@ -11,6 +11,7 @@ import json import os SAMPLE_INTERVAL = int(os.environ.get("SAMPLE_INTERVAL", "60")) +RESET_INTERVAL = int(os.environ.get("RESET_INTERVAL", "3600")) # Reset detector every N seconds (default: 1h) TARGET_DURATION = int(os.environ.get("TARGET_DURATION", str(86400))) # 24h OUTPUT_PATH = os.environ.get("BACKGROUND_PATH", "/data/background_24h.npy") SNAPSHOT_PATH = os.environ.get("SNAPSHOT_PATH", "/data/background_snapshot.json") @@ -18,6 +19,22 @@ SNAPSHOT_PATH = os.environ.get("SNAPSHOT_PATH", "/data/background_snapshot.json" BG_COUNTS = np.zeros(1024, dtype=np.float64) BG_LIVE_TIME = 0.0 device = None +last_counts = None +last_live_time = None +last_reset_time = 0 + +# Resume from previous snapshot if it exists +if os.path.exists(SNAPSHOT_PATH): + try: + with open(SNAPSHOT_PATH) as _f: + _prev = json.load(_f) + _prev_spectrum = _prev.get("spectrum", []) + if len(_prev_spectrum) == 1024: + BG_COUNTS = np.array(_prev_spectrum, dtype=np.float64) + BG_LIVE_TIME = float(_prev.get("live_time_s", 0)) + print(f"Snapshot anterieur charge : {BG_LIVE_TIME:.0f}s live, {BG_COUNTS.sum():.0f} coups") + except Exception as _e: + print(f"Impossible de charger le snapshot anterieur : {_e}") def save_snapshot(): """Save a human-readable snapshot of current background.""" @@ -46,7 +63,18 @@ print(f"Capture du bruit de fond pendant {TARGET_DURATION/3600:.0f}h...") print("Assurez-vous qu'aucune source radioactive n'est a proximite du detecteur.") print() -start = time.time() +start_offset = 0 +if os.path.exists(SNAPSHOT_PATH): + try: + with open(SNAPSHOT_PATH) as _f: + _prev = json.load(_f) + start_offset = float(_prev.get("elapsed_hours", 0)) * 3600 + except Exception: + pass + +start = time.time() - start_offset +last_reset_time = time.time() + while (time.time() - start) < TARGET_DURATION: time.sleep(SAMPLE_INTERVAL) try: @@ -55,12 +83,41 @@ while (time.time() - start) < TARGET_DURATION: device = RadiaCode() device.spectrum_reset() + last_counts = None + last_live_time = None + last_reset_time = time.time() print("Radiacode connecte.") spectrum = device.spectrum() - BG_COUNTS += np.array(spectrum.counts, dtype=np.float64) - BG_LIVE_TIME += spectrum.duration.total_seconds() - device.spectrum_reset() + counts = np.array(spectrum.counts, dtype=np.float64) + live_time = spectrum.duration.total_seconds() + + # Compute delta since last read (avoid double-counting) + if last_counts is not None and last_live_time is not None: + delta_counts = counts - last_counts + delta_live_time = live_time - last_live_time + # If detector was reset externally, delta would be negative + if delta_counts.sum() < 0 or delta_live_time < 0: + delta_counts = counts + delta_live_time = live_time + BG_COUNTS += delta_counts + BG_LIVE_TIME += max(delta_live_time, 0) + else: + BG_COUNTS += counts + BG_LIVE_TIME += live_time + + last_counts = counts.copy() + last_live_time = live_time + + # Only reset detector spectrum at RESET_INTERVAL + now = time.time() + if now - last_reset_time >= RESET_INTERVAL: + device.spectrum_reset() + last_counts = None + last_live_time = None + last_reset_time = now + print(f" → Reset détecteur (intervalle {RESET_INTERVAL}s)") + elapsed = time.time() - start cps = BG_COUNTS.sum() / BG_LIVE_TIME if BG_LIVE_TIME > 0 else 0 print( @@ -72,6 +129,8 @@ while (time.time() - start) < TARGET_DURATION: except Exception as e: print(f"\nErreur : {e}, reconnexion...") device = None + last_counts = None + last_live_time = None os.makedirs(os.path.dirname(OUTPUT_PATH) if os.path.dirname(OUTPUT_PATH) else ".", exist_ok=True) np.save( diff --git a/docker-compose.yml b/docker-compose.yml index eed5281..b91cc48 100644 --- a/docker-compose.yml +++ b/docker-compose.yml @@ -60,7 +60,7 @@ services: context: ./web dockerfile: Dockerfile ports: - - "8080:8080" + - "8000:8080" volumes: - ./data:/data:ro - ./logs:/logs:ro @@ -74,4 +74,4 @@ services: - ISOTOPE_INDEX_PATH=/models/vega_isotope_index.txt - ENERGY_CALIBRATION_OFFSET=0.33 - ENERGY_CALIBRATION_SLOPE=2.97 - restart: unless-stopped \ No newline at end of file + restart: unless-stopped diff --git a/train/vega_ml/synthetic_spectra/physics/spectrum_physics.py b/train/vega_ml/synthetic_spectra/physics/spectrum_physics.py index 8131020..5cb6a2a 100644 --- a/train/vega_ml/synthetic_spectra/physics/spectrum_physics.py +++ b/train/vega_ml/synthetic_spectra/physics/spectrum_physics.py @@ -324,13 +324,16 @@ def generate_realistic_continuum( Generate realistic CsI(Tl) background continuum shape. Calibrated against real Radiacode 103 background measurements. - Produces the characteristic asymmetric hump at ~110 keV and - Compton-like tail that simple exponentials miss. + Produces the characteristic asymmetric hump at ~110 keV with + housing absorption at low energy, Compton plateau, and proper + high-energy falloff. Shape components: - - Asymmetric hump centered at ~110 keV (sigma_left=55, sigma_right=50 keV) - - Compton continuum: 0.45*exp(-E/240) + 0.04*exp(-E/700) - - Noise floor at 0.8% of peak + - Asymmetric hump at ~110 keV (sigma_left=48, sigma_right=80 keV) + - Housing absorption below ~40 keV: 1 - exp(-E/30) + - Compton plateau around 200-260 keV from Pb-214/Bi-214 scatter + - Compton tail: 0.38*exp(-E/170) + 0.06*exp(-E/500) + - Noise floor at 0.3% of peak Args: energy_bins: Array of energy bin centers (keV) @@ -342,24 +345,31 @@ def generate_realistic_continuum( """ E = energy_bins - # Asymmetric hump at ~110 keV (low-energy scatter peak in CsI(Tl)) + # Asymmetric hump at ~110 keV + # Left side sharper (sigma=48), right side broader with Compton shoulder (sigma=80) hump_center = 110.0 - sigma_left = 55.0 # Broader on the low-energy side - sigma_right = 50.0 # Narrower on the high-energy side hump = np.where( E <= hump_center, - np.exp(-0.5 * ((E - hump_center) / sigma_left) ** 2), - np.exp(-0.5 * ((E - hump_center) / sigma_right) ** 2), + np.exp(-0.5 * ((E - hump_center) / 48.0) ** 2), + np.exp(-0.5 * ((E - hump_center) / 80.0) ** 2), ) + # Housing absorption at very low energy (< ~40 keV) + absorption = 1.0 - np.exp(-E / 30.0) + # Compton continuum tail - tail = 0.45 * np.exp(-E / 240.0) + 0.04 * np.exp(-E / 700.0) + tail = 0.38 * np.exp(-E / 170.0) + 0.06 * np.exp(-E / 500.0) - # Noise floor (low-level baseline) - noise_floor = 0.008 + # Compton plateau around 200-260 keV (Pb-214/Bi-214 scatter) + compton_edge = np.maximum(0, 1.0 - ((E - 180.0) / 150.0) ** 2) + compton_edge[E > 330] = 0 + compton_plateau = 0.12 * compton_edge - # Combine shape components - continuum = hump + tail + noise_floor + # Noise floor + noise_floor = 0.003 + + # Combine continuum with absorption + continuum = (hump + tail + compton_plateau) * absorption + noise_floor # Normalize to target total counts if continuum.sum() > 0 and total_counts > 0: diff --git a/web/app/bg_calibration.py b/web/app/bg_calibration.py new file mode 100644 index 0000000..f4cde7b --- /dev/null +++ b/web/app/bg_calibration.py @@ -0,0 +1,208 @@ +""" +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 two approaches: +1. Spline-based: non-parametric, automatically fits any shape +2. Parametric: for the /fit endpoint (comparison with measured data) + +The spline approach is preferred — it uses scipy's smoothing spline with +Generalized Cross-Validation to automatically find the right smoothness, +after iterative peak subtraction. +""" + +import json +import numpy as np +from pathlib import Path +from scipy.interpolate import make_smoothing_spline +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 _smooth(y): + window = min(51, len(y) // 10 * 2 + 1) + if window < 5: + window = 5 + return savgol_filter(y, window_length=window, polyorder=3) + + +def _subtract_peaks(energy_axis, smoothed_cps): + """Iteratively estimate and subtract photopeak contributions.""" + continuum = smoothed_cps.copy() + peak_amplitudes = [] + + for line_energy, _ in PHOTOPEAK_LINES: + sig = _sigma_keV(line_energy) + idx = np.argmin(np.abs(energy_axis - line_energy)) + n_sigma = max(int(2 * sig / 2.97), 3) + off_lo = continuum[max(0, idx - 3 * n_sigma):max(1, idx - n_sigma)] + off_hi = continuum[min(len(continuum), idx + n_sigma):min(len(continuum), idx + 3 * n_sigma)] + off_peak = np.concatenate([off_lo, off_hi]) + local_bg = np.median(off_peak) if len(off_peak) > 0 else 0 + peak_height = continuum[idx] - local_bg + + if peak_height > 0: + amplitude = peak_height * sig * np.sqrt(2 * np.pi) + gaussian = amplitude * np.exp(-0.5 * ((energy_axis - line_energy) / sig) ** 2) / (sig * np.sqrt(2 * np.pi)) + continuum -= gaussian + continuum = np.maximum(continuum, 0) + peak_amplitudes.append({"energy_keV": line_energy, "amplitude": float(max(0, peak_height) * sig * np.sqrt(2 * np.pi)) if peak_height > 0 else 0.0}) + + return continuum, peak_amplitudes + + +def calibrate_spline(measured_cps, energy_axis): + """ + Fit a smoothing spline to the peak-subtracted continuum. + + Uses scipy's make_smoothing_spline with GCV (Generalized Cross-Validation) + to automatically find the optimal smoothing parameter. + + Returns a dict with the fitted spline evaluated at all channels. + """ + E = energy_axis + y_smooth = _smooth(measured_cps) + continuum, peak_amplitudes = _subtract_peaks(E, y_smooth) + + # Ensure positive values for spline fitting + continuum = np.maximum(continuum, 0) + + # Use log-space for better fit at low-signal high-energy region + # Add small offset to avoid log(0) + offset = continuum[continuum > 0].min() * 0.1 if (continuum > 0).any() else 1e-6 + log_continuum = np.log(continuum + offset) + + # Fit smoothing spline in log-space (GCV auto-selects lambda) + try: + spline = make_smoothing_spline(E, log_continuum) + log_fit = spline(E) + # Convert back from log-space + fit_cps = np.exp(log_fit) - offset + fit_cps = np.maximum(fit_cps, 0) + except Exception as e: + return {"error": str(e)} + + # Quality metrics + residuals = continuum - fit_cps + 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": fit_cps, + "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 smoothing spline. + Returns both spline-based fit and parameters for the /fit endpoint. + """ + result = calibrate_spline(measured_cps, energy_axis) + if "error" in result: + return result + + # The spline result is the continuum CPS array + return { + "params": {}, # Non-parametric model, no params + "continuum_cps": result["continuum_cps"], + "peak_amplitudes": result["peak_amplitudes"], + "r_squared": result["r_squared"], + "residuals_rms": result["residuals_rms"], + "method": "smoothing_spline_gcv", + } + + +def build_calibrated_continuum(energy_axis, total_counts, params): + """Build the continuum from calibrated parameters.""" + if "continuum_cps" in params: + # Spline-based: already have the CPS array + 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 + + if _CALIBRATION_PATH.exists(): + try: + with open(_CALIBRATION_PATH) as f: + _cached_result = json.load(f) + return _cached_result + except Exception: + pass + + 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"], + } + + _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) + + return _cached_result \ No newline at end of file diff --git a/web/app/noise.py b/web/app/noise.py new file mode 100644 index 0000000..be66e4a --- /dev/null +++ b/web/app/noise.py @@ -0,0 +1,107 @@ +""" +Detector-agnostic continuum extraction for gamma-ray spectra. + +Extracts the detector's intrinsic response curve (continuum) from measured +background data. Isotope photopeaks are subtracted, then the residual is +smoothed to produce a clean continuum shape that reflects only the detector's +physics — no isotope signatures. + +Works with any scintillator or semiconductor detector. +""" + +import numpy as np +from scipy.ndimage import gaussian_filter1d + + +# Common environmental isotope lines (keV) — subtracted regardless of detector. +_ENV_PEAKS = [ + (241.0, 0.04), + (295.22, 0.1842), + (351.93, 0.3560), + (609.31, 0.4549), + (911.20, 0.2580), + (1120.29, 0.1492), + (1460.83, 0.1066), + (1764.49, 0.1531), + (2614.51, 0.3586), +] + +_E_OFFSET = 0.33 +_E_SLOPE = 2.97 + + +def _sigma_ch(E_keV): + """Peak sigma in channels at energy E_keV (sqrt(E) resolution scaling).""" + fwhm_keV = 0.08 * E_keV * (E_keV / 662.0) ** 0.5 + sigma_keV = fwhm_keV / 2.355 + return max(sigma_keV / _E_SLOPE, 2.0) + + +def _subtract_peaks(counts, energy_axis): + """Remove known isotope photopeaks from spectrum.""" + continuum = counts.copy() + channels = np.arange(len(counts), dtype=np.float64) + + for line_energy, _ in _ENV_PEAKS: + idx = int(np.argmin(np.abs(energy_axis - line_energy))) + if idx < 0 or idx >= len(counts): + continue + + 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(counts), idx + far) + hi_end = min(len(counts), 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) + + +def extract_continuum(counts, energy_axis=None): + """Extract the detector's intrinsic response continuum. + + Removes isotope photopeaks, then smooths with a wide Gaussian filter + to produce a clean curve showing only the detector's continuum shape. + + Parameters + ---------- + counts : array + Raw accumulated counts per channel. + energy_axis : array, optional + Energy axis in keV. + + Returns + ------- + array — smooth continuum (peak-subtracted, Gaussian-smoothed) + """ + counts = np.asarray(counts, dtype=np.float64) + n_channels = len(counts) + + if energy_axis is None: + energy_axis = _E_OFFSET + _E_SLOPE * np.arange(n_channels, dtype=np.float64) + + continuum = _subtract_peaks(counts, energy_axis) + + # Wide Gaussian smooth (sigma ~1.5% of channels ≈ 45 keV) + sigma = max(15, n_channels // 60) + continuum_smooth = gaussian_filter1d(continuum, sigma=sigma) + continuum_smooth = np.maximum(continuum_smooth, 0) + + return continuum_smooth \ No newline at end of file diff --git a/web/app/routers/background.py b/web/app/routers/background.py index b76f0dc..b68c2ca 100644 --- a/web/app/routers/background.py +++ b/web/app/routers/background.py @@ -1,7 +1,8 @@ import json from fastapi import APIRouter, HTTPException -from app.config import BACKGROUND_SNAPSHOT_PATH, BACKGROUND_PATH, energy_axis, NUM_CHANNELS -from app.theoretical_bg import generate_theoretical_bg, generate_continuum_only +from app.config import BACKGROUND_SNAPSHOT_PATH, BACKGROUND_PATH, energy_axis, NUM_CHANNELS, ENERGY_OFFSET, ENERGY_SLOPE +from app.theoretical_bg import generate_continuum_only +from app.noise import extract_continuum import numpy as np router = APIRouter() @@ -80,16 +81,100 @@ async def get_background_reference(): } -@router.get("/theoretical") -async def get_theoretical_bg(cps: float = 6.0, live_time_s: float = 3600.0): - """Theoretical natural background spectrum (K-40, U-238 chain, Th-232 chain).""" - return generate_theoretical_bg(cps=cps, live_time_s=live_time_s) - - @router.get("/continuum") async def get_continuum(cps: float = 6.0, live_time_s: float = 3600.0): - """CsI(Tl) continuum shape only (hump + Compton tail, no photopeaks, no noise). + """CsI(Tl) detector response continuum only (no photopeaks, no noise).""" + return generate_continuum_only(cps=cps, live_time_s=live_time_s) - Matches the model used in training (generate_realistic_continuum). + +@router.get("/fit") +async def fit_background(): + """Fit the parametric CsI(Tl) detector response model to measured background data. + + Returns the fitted curve, parameters, and quality metrics. """ - return generate_continuum_only(cps=cps, live_time_s=live_time_s) \ No newline at end of file + from app.bg_calibration import calibrate_background, build_calibrated_continuum + + # Load measured data + 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: + raise HTTPException(status_code=404, detail="No measured background available for fitting") + + channels = np.arange(NUM_CHANNELS, dtype=np.float64) + e_axis = ENERGY_OFFSET + ENERGY_SLOPE * channels + + # Run calibration + measured_cps = measured_counts / live_time + result = calibrate_background(measured_cps, e_axis) + + if "error" in result: + raise HTTPException(status_code=500, detail=f"Fitting failed: {result['error']}") + + # Build fitted curve at same scale as measured + fitted_counts = build_calibrated_continuum(e_axis, measured_counts.sum(), result) + + return { + "energy_kev": [round(float(E), 2) for E in e_axis], + "measured_counts": [round(float(c), 1) for c in measured_counts], + "fitted_counts": [round(float(c), 1) for c in fitted_counts], + "method": result.get("method", "spline"), + "r_squared": result["r_squared"], + "residuals_rms": result["residuals_rms"], + "live_time_s": round(live_time, 1), + } + + +@router.get("/noise") +async def get_background_noise(): + """Detector's intrinsic continuum curve (isotope peaks subtracted). + + Returns the smooth detector response shape without any isotope + photopeak signatures. Works with any detector type. + """ + counts = None + + if BACKGROUND_PATH.exists(): + try: + bg_data = np.load(str(BACKGROUND_PATH), allow_pickle=True).item() + counts = bg_data["counts"].astype(np.float64)[:NUM_CHANNELS] + except Exception: + pass + + if counts is None and BACKGROUND_SNAPSHOT_PATH.exists(): + try: + with open(BACKGROUND_SNAPSHOT_PATH) as f: + snapshot = json.load(f) + counts = np.array(snapshot.get("spectrum", [])[:NUM_CHANNELS], dtype=np.float64) + except Exception: + pass + + if counts is None: + raise HTTPException(status_code=404, detail="No background data available") + + channels = np.arange(NUM_CHANNELS, dtype=np.float64) + e_axis = ENERGY_OFFSET + ENERGY_SLOPE * channels + continuum = extract_continuum(counts, energy_axis=e_axis) + + return { + "energy_kev": [round(float(E), 2) for E in e_axis], + "counts": [round(float(c), 1) for c in continuum], + } \ No newline at end of file diff --git a/web/app/theoretical_bg.py b/web/app/theoretical_bg.py index 4cb8b7a..30ba036 100644 --- a/web/app/theoretical_bg.py +++ b/web/app/theoretical_bg.py @@ -1,139 +1,74 @@ """ -Theoretical natural background spectrum for CsI(Tl) detectors (Radiacode 103). +CsI(Tl) detector response continuum for Radiacode 103. -Shape calibrated against real Radiacode 103 background measurements. -The CsI(Tl) crystal (1 cm³, 8.4% FWHM) produces a spectrum with: -- A dominant low-energy hump peaking around 100-120 keV -- Exponential decay at higher energies -- Subtle photopeaks from natural isotopes +Models ONLY the detector's noise continuum. Photopeaks from environmental +isotopes depend on measurement location and are NOT included. + +Auto-calibrated from measured background using smoothing spline (GCV) +when available. Falls back to a simple parametric model otherwise. """ import numpy as np from app.config import ENERGY_OFFSET, ENERGY_SLOPE, NUM_CHANNELS -# Photopeak lines: (energy_keV, relative_weight) -# Weights tuned so peaks are visible above local continuum at typical CPS -NATURAL_BG_LINES = [ - (295.22, 0.10), # Pb-214 - (351.93, 0.18), # Pb-214 - (609.31, 0.15), # Bi-214 - (911.20, 0.08), # Ac-228 - (968.97, 0.05), # Ac-228 - (1120.29, 0.06), # Bi-214 - (1460.83, 0.12), # K-40 - (1764.49, 0.08), # Bi-214 - (2614.51, 0.18), # Tl-208 -] - - -def _gaussian(x, center, sigma, amplitude): - return amplitude * np.exp(-0.5 * ((x - center) / sigma) ** 2) - - -def generate_theoretical_bg(cps: float = 6.0, live_time_s: float = 3600.0): - channels = np.arange(NUM_CHANNELS, dtype=np.float64) - energy_axis = ENERGY_OFFSET + ENERGY_SLOPE * channels - total_counts = cps * live_time_s - - # ── 1. Main hump: asymmetric peak at ~105 keV ── - # Real data: rises from ~60 at 10keV to ~280 at 100-120keV, then falls - hump_center = 110.0 - hump = np.zeros(NUM_CHANNELS, dtype=np.float64) - low_mask = energy_axis <= hump_center - hump[low_mask] = _gaussian(energy_axis[low_mask], hump_center, 55.0, 1.0) - hump[~low_mask] = _gaussian(energy_axis[~low_mask], hump_center, 50.0, 1.0) - - # ── 2. Compton continuum tail ── - # Real data: ~136@200, ~80@250, ~44@295, ~14@400, ~5@600 - tail = 0.45 * np.exp(-energy_axis / 240) + 0.04 * np.exp(-energy_axis / 700) - - # ── 3. Low-energy noise floor ── - noise_floor = 0.008 - - # ── 4. Combine continuum ── - continuum = hump + tail + noise_floor - - # ── 5. Photopeaks ── - # CsI(Tl) 8.4% FWHM at 662 keV, scaling as sqrt(E) - # sigma(E) = FWHM(E) / 2.355 = 0.084 * sqrt(E * 662) / 662 / 2.355 - # Simplified: sigma = 23.6 * sqrt(E/662) keV - def sigma_keV(E): - return max(12.0, 23.6 * np.sqrt(max(E, 1.0) / 662.0)) - - peak_frac = 0.08 # 8% of total counts in resolved photopeaks - total_weight = sum(w for _, w in NATURAL_BG_LINES) - - peaks = np.zeros(NUM_CHANNELS, dtype=np.float64) - for line_energy, weight in NATURAL_BG_LINES: - sig = sigma_keV(line_energy) - peak_counts = total_counts * peak_frac * (weight / total_weight) - amplitude = peak_counts / (sig * np.sqrt(2 * np.pi)) - peaks += _gaussian(energy_axis, line_energy, sig, amplitude) - - # ── 6. Combine and normalize ── - raw = continuum + peaks / total_counts # peaks normalized later - raw *= total_counts / raw.sum() - - # ── 7. Poisson-like noise ── - rng = np.random.default_rng(42) - noise = rng.normal(0, 1, NUM_CHANNELS) * np.sqrt(np.maximum(raw, 1.0)) * 0.25 - raw += noise - - # Floor at 0.9 for log scale - spectrum = np.clip(raw, 0.9, None) - - key_lines = [ - (295.22, "Pb-214"), (351.93, "Pb-214"), - (609.31, "Bi-214"), (911.20, "Ac-228"), - (1120.29, "Bi-214"), (1460.83, "K-40"), - (1764.49, "Bi-214"), (2614.51, "Tl-208"), - ] - - return { - "energy_kev": [round(float(E), 2) for E in energy_axis], - "counts": [round(float(c), 1) for c in spectrum], - "cps": round(cps, 2), - "live_time_s": round(live_time_s, 1), - "lines": [ - {"energy_keV": E, "name": name} for E, name in key_lines - ], - } +def _get_continuum_cps(): + """Try to load calibrated spline continuum from measured data.""" + try: + from app.bg_calibration import load_or_calibrate + calibrated = load_or_calibrate() + if calibrated and "continuum_cps" in calibrated: + return np.array(calibrated["continuum_cps"]) + except Exception: + pass + return None def generate_continuum_only(cps: float = 6.0, live_time_s: float = 3600.0): - """Generate only the CsI(Tl) continuum shape (no photopeaks, no noise). - - This matches the model used in training (generate_realistic_continuum in - spectrum_physics.py) for direct comparison with measured backgrounds. - """ + """Detector response continuum only (no photopeaks, no noise).""" channels = np.arange(NUM_CHANNELS, dtype=np.float64) energy_axis = ENERGY_OFFSET + ENERGY_SLOPE * channels total_counts = cps * live_time_s - # Asymmetric hump at ~110 keV - hump_center = 110.0 - hump = np.where( - energy_axis <= hump_center, - np.exp(-0.5 * ((energy_axis - hump_center) / 55.0) ** 2), - np.exp(-0.5 * ((energy_axis - hump_center) / 50.0) ** 2), - ) + # Try calibrated spline first + continuum_cps = _get_continuum_cps() - # Compton continuum tail - tail = 0.45 * np.exp(-energy_axis / 240.0) + 0.04 * np.exp(-energy_axis / 700.0) - - # Noise floor - noise_floor = 0.008 - - continuum = hump + tail + noise_floor - - # Normalize to target total counts - if continuum.sum() > 0 and total_counts > 0: - continuum *= total_counts / continuum.sum() + if continuum_cps is not None and len(continuum_cps) == NUM_CHANNELS: + # Scale calibrated CPS to match requested total counts + continuum = continuum_cps.copy() + if continuum.sum() > 0: + continuum *= total_counts / continuum.sum() + else: + # Fallback: simple parametric model + continuum = _fallback_continuum(energy_axis, total_counts) return { "energy_kev": [round(float(E), 2) for E in energy_axis], "counts": [round(float(c), 1) for c in continuum], "cps": round(cps, 2), "live_time_s": round(live_time_s, 1), - } \ No newline at end of file + } + + +def _fallback_continuum(energy_axis, total_counts): + """Simple parametric fallback when no measured data available.""" + E = energy_axis + + # Asymmetric hump + hump_center, sigma_left, tail_decay_right = 110.0, 40.0, 100.0 + left = np.exp(-0.5 * ((E - hump_center) / sigma_left) ** 2) + right = np.exp(-(E - hump_center) / tail_decay_right) + hump = np.where(E <= hump_center, left, right) + + # Housing absorption + absorption = 1.0 * (1.0 - np.exp(-E / 20.0)) + + # Compton tail + compton = 0.5 / (np.maximum(E, 1.0) + 15.0) ** 1.3 + + continuum = (hump + compton) * absorption + + if continuum.sum() > 0 and total_counts > 0: + continuum *= total_counts / continuum.sum() + + return continuum \ No newline at end of file diff --git a/web/requirements.txt b/web/requirements.txt index 0a5c747..536bf56 100644 --- a/web/requirements.txt +++ b/web/requirements.txt @@ -1,3 +1,4 @@ fastapi>=0.104.0 uvicorn[standard]>=0.24.0 -numpy>=1.24.0 \ No newline at end of file +numpy>=1.24.0 +scipy>=1.10.0 \ No newline at end of file diff --git a/web/static/css/style.css b/web/static/css/style.css index e622044..c260ee1 100644 --- a/web/static/css/style.css +++ b/web/static/css/style.css @@ -78,6 +78,8 @@ main { padding: 16px; } border-radius: 8px; padding: 12px; margin-bottom: 12px; + height: 450px; + position: relative; } .controls { @@ -108,6 +110,63 @@ main { padding: 16px; } .controls button:hover { background: var(--accent-bright); color: #000; } +.btn-small { + background: var(--accent); + color: var(--text); + border: 1px solid #444; + padding: 4px 10px; + border-radius: 4px; + cursor: pointer; + font-size: 0.8em; + margin-left: auto; +} +.btn-small:hover { background: var(--accent-bright); color: #000; } + +.chart-container.fullscreen { + position: fixed; + top: 0; left: 0; right: 0; bottom: 0; + z-index: 1000; + background: var(--bg); + padding: 20px; + margin: 0; + border-radius: 0; + display: flex; + flex-direction: column; +} +.chart-container.fullscreen canvas { + flex: 1; +} +.exit-fullscreen-btn { + display: none; + position: absolute; + top: 10px; + right: 14px; + z-index: 1001; + background: rgba(255,255,255,0.15); + color: #fff; + border: none; + border-radius: 50%; + width: 36px; + height: 36px; + font-size: 1.2em; + cursor: pointer; + line-height: 1; +} +.chart-container.fullscreen .exit-fullscreen-btn { + display: flex; + align-items: center; + justify-content: center; +} +.exit-fullscreen-btn:hover { + background: rgba(255,255,255,0.3); +} + +.chart-header { + display: flex; + justify-content: flex-end; + margin-bottom: 4px; +} + #isotopes-table, #peaks-table { background: var(--bg-card); border-radius: 8px; diff --git a/web/static/index.html b/web/static/index.html index c6c3cdb..2735fde 100644 --- a/web/static/index.html +++ b/web/static/index.html @@ -7,6 +7,8 @@ + +