diff --git a/detect/radiacode_monitor.py b/detect/radiacode_monitor.py
index 792d890..1ffbdb8 100644
--- a/detect/radiacode_monitor.py
+++ b/detect/radiacode_monitor.py
@@ -30,6 +30,56 @@ CPS_LOG_PATH = os.environ.get("CPS_LOG_PATH", "/data/cps_log.jsonl")
ENERGY_OFFSET = float(os.environ.get("ENERGY_CALIBRATION_OFFSET", "0.33"))
ENERGY_SLOPE = float(os.environ.get("ENERGY_CALIBRATION_SLOPE", "2.97"))
+# CsI(Tl) non-linear response correction
+# CsI(Tl) produces more light per keV at low energies, shifting peaks to higher
+# apparent energies. Model: E_apparent = E_true * (1 + alpha * exp(-E_true/beta))
+# Calibrated from Am-241 (59.5 keV appears at ~71.6 keV) and K-40 (correct at 1460.8 keV).
+CSI_NONLINEAR_ALPHA = float(os.environ.get("CSI_NONLINEAR_ALPHA", "0.37"))
+CSI_NONLINEAR_BETA = float(os.environ.get("CSI_NONLINEAR_BETA", "100.0"))
+
+
+def correct_csilinear_energy(spectrum_rate, num_channels=1023):
+ """Apply inverse CsI(Tl) non-linear response correction to spectrum channels.
+
+ CsI(Tl) has non-proportional scintillation response at low energies,
+ causing peaks to appear at higher channels than their true energy position.
+ This function remaps channels so that peaks appear at their theoretical
+ energy positions, matching what the model was trained on.
+
+ For each output channel j (true energy position), we find the input
+ channel i (apparent energy position) where the detector actually placed
+ counts for that true energy.
+
+ Args:
+ spectrum_rate: Array of 1023 channel count rates
+ num_channels: Number of channels
+
+ Returns:
+ Corrected spectrum with peaks at theoretical energy positions
+ """
+ alpha = CSI_NONLINEAR_ALPHA
+ beta = CSI_NONLINEAR_BETA
+
+ # For each output channel j, compute the apparent energy where
+ # counts for true energy E_true(j) actually appear
+ output_channels = np.arange(num_channels, dtype=np.float64)
+ e_true = ENERGY_OFFSET + ENERGY_SLOPE * output_channels
+
+ # Forward model: E_apparent = E_true * (1 + alpha * exp(-E_true / beta))
+ e_apparent = e_true * (1 + alpha * np.exp(-e_true / beta))
+
+ # Input channel where the detector placed counts for this true energy
+ source_channels = (e_apparent - ENERGY_OFFSET) / ENERGY_SLOPE
+ source_channels = np.clip(source_channels, 0, num_channels - 1.001)
+
+ # Linear interpolation from source channels
+ lower = np.floor(source_channels).astype(int)
+ upper = np.minimum(lower + 1, num_channels - 1)
+ frac = source_channels - lower
+
+ corrected = spectrum_rate[lower] * (1 - frac) + spectrum_rate[upper] * frac
+ return corrected
+
# Logging
logging.basicConfig(
level=logging.INFO,
@@ -46,10 +96,10 @@ class RadiacodeMonitor:
def __init__(self):
# Charger le modèle PyTorch
device_str = os.environ.get("VEGA_DEVICE", "cpu")
- self.device = torch.device(device_str)
+ self.torch_device = torch.device(device_str)
- log.info(f"Chargement du modèle depuis {MODEL_PATH} sur {self.device}...")
- checkpoint = torch.load(MODEL_PATH, map_location=self.device, weights_only=False)
+ log.info(f"Chargement du modèle depuis {MODEL_PATH} sur {self.torch_device}...")
+ checkpoint = torch.load(MODEL_PATH, map_location=self.torch_device, weights_only=False)
# Importer VegaModel (depuis le volume monté)
vega_ml_path = os.environ.get("VEGA_ML_PATH", "/models/vega_ml")
@@ -86,42 +136,62 @@ class RadiacodeMonitor:
else:
log.warning(f"Pas de fichier background : {BACKGROUND_PATH}")
+ # Connexion persistante au Radiacode
+ self._rc = None
+ self.reconnect_backoff = 0
+
# Compteurs cumulés
self.cumulated_counts = np.zeros(1024, dtype=np.float64)
self.cumulated_live_time = 0.0
self.last_report_date = None
self.connected = False
- def try_connect(self):
- """Tente de se connecter au Radiacode. Retourne le device ou None."""
+ def _connect(self):
+ """Tente d'établir une connexion persistante au Radiacode."""
try:
from radiacode import RadiaCode
- device = RadiaCode()
- log.info("Radiacode connecté")
+ self._rc = RadiaCode()
self.connected = True
- return device
+ self.reconnect_backoff = 0
+ log.info("Radiacode connecté")
+ return True
except Exception as e:
- log.debug(f"Détecteur non disponible : {e}")
+ self._rc = None
self.connected = False
- return None
+ self.reconnect_backoff = min(self.reconnect_backoff + 1, 10)
+ log.debug(f"Détecteur non disponible (retry dans {self.reconnect_backoff} cycles) : {e}")
+ return False
+
+ def _disconnect(self):
+ """Ferme la connexion au Radiacode."""
+ if self._rc is not None:
+ try:
+ del self._rc
+ except Exception:
+ pass
+ self._rc = None
+ self.connected = False
def sample_once(self):
"""Échantillonne une fois. Retourne True si succès."""
- device = None
- try:
- device = self.try_connect()
- if device is None:
+ # Établir la connexion si nécessaire
+ if self._rc is None:
+ if self.reconnect_backoff > 0:
+ self.reconnect_backoff -= 1
+ return False
+ if not self._connect():
return False
- spectrum = device.spectrum()
+ try:
+ spectrum = self._rc.spectrum()
counts = np.array(spectrum.counts, dtype=np.float64)
live_time = spectrum.duration.total_seconds()
if live_time > 0 and counts.sum() > 0:
self.cumulated_counts += counts
self.cumulated_live_time += live_time
- device.spectrum_reset()
+ self._rc.spectrum_reset()
log.info(
f"Échantillon : {counts.sum():.0f} coups en {live_time:.1f}s "
f"(cumul : {self.cumulated_live_time/3600:.1f}h)"
@@ -129,14 +199,9 @@ class RadiacodeMonitor:
return True
return False
except Exception as e:
- log.warning(f"Erreur échantillonnage : {e}")
+ log.warning(f"Erreur échantillonnage, reconnexion : {e}")
+ self._disconnect()
return False
- finally:
- if device:
- try:
- del device
- except Exception:
- pass
def save_state(self):
"""Ecrit l'etat actuel du moniteur dans un fichier JSON atomique."""
@@ -147,10 +212,10 @@ class RadiacodeMonitor:
if self.cumulated_live_time > 0:
rate = self.cumulated_counts / self.cumulated_live_time
if self.bg_counts is not None and self.bg_live_time is not None:
- bg_rate = self.bg_counts / self.bg_live_time
- net_rate = np.clip(rate - bg_rate, 0, None)
+ bg_rate = self.bg_counts[:1023] / self.bg_live_time
+ net_rate = np.clip(rate[:1023] - bg_rate, 0, None)
else:
- net_rate = rate
+ net_rate = rate[:1023]
isotopes = self.run_inference(net_rate)
state = {
@@ -190,11 +255,15 @@ class RadiacodeMonitor:
def run_inference(self, spectrum_rate):
"""Exécute l'inférence PyTorch sur le spectre cumulé."""
if spectrum_rate.max() > 0:
- normalized = spectrum_rate / spectrum_rate.max()
+ # Apply CsI(Tl) non-linear correction so peaks appear
+ # at theoretical energy positions (matching training data)
+ corrected = correct_csilinear_energy(spectrum_rate)
+ log_spectrum = np.log1p(np.maximum(corrected, 0))
+ normalized = log_spectrum / log_spectrum.max()
else:
return []
- tensor = torch.tensor(normalized, dtype=torch.float32).unsqueeze(0).to(self.device)
+ tensor = torch.tensor(normalized, dtype=torch.float32).unsqueeze(0).to(self.torch_device)
with torch.no_grad():
logits, activities = self.model(tensor)
@@ -227,10 +296,10 @@ class RadiacodeMonitor:
rate = self.cumulated_counts / self.cumulated_live_time
if self.bg_counts is not None and self.bg_live_time is not None:
- bg_rate = self.bg_counts / self.bg_live_time
- net_rate = np.clip(rate - bg_rate, 0, None)
+ bg_rate = self.bg_counts[:1023] / self.bg_live_time
+ net_rate = np.clip(rate[:1023] - bg_rate, 0, None)
else:
- net_rate = rate
+ net_rate = rate[:1023]
results = self.run_inference(net_rate)
@@ -280,7 +349,7 @@ class RadiacodeMonitor:
log.info("Radiacode 103 — Moniteur d'isotopes")
log.info("=" * 50)
log.info(f"Modèle : {MODEL_PATH}")
- log.info(f"Device : {self.device}")
+ log.info(f"Device : {self.torch_device}")
log.info(f"Isotopes : {self.isotope_index.num_isotopes}")
log.info(
f"Background : {'chargé' if self.bg_counts is not None else 'non disponible'}"
diff --git a/docker-compose.yml b/docker-compose.yml
index b91cc48..31d75a4 100644
--- a/docker-compose.yml
+++ b/docker-compose.yml
@@ -20,11 +20,11 @@ services:
- MODEL_DIR=/models
- NUM_SAMPLES=50000
- EPOCHS=100
- - BATCH_SIZE=64
+ - BATCH_SIZE=32
- LEARNING_RATE=0.001
- DETECTOR=radiacode_103
- - MIN_DURATION=43200
- - MAX_DURATION=86400
+ - MIN_DURATION=30
+ - MAX_DURATION=300
- SEED=42
- MEASURED_BACKGROUND_PATH=/data/background_24h.npy
restart: "no"
@@ -53,6 +53,8 @@ services:
- REPORT_HOUR=0
- MIN_LIVE_TIME=3600
- THRESHOLD=0.5
+ - CSI_NONLINEAR_ALPHA=0.37
+ - CSI_NONLINEAR_BETA=100.0
restart: unless-stopped
web:
@@ -74,4 +76,6 @@ services:
- ISOTOPE_INDEX_PATH=/models/vega_isotope_index.txt
- ENERGY_CALIBRATION_OFFSET=0.33
- ENERGY_CALIBRATION_SLOPE=2.97
+ - CSI_NONLINEAR_ALPHA=0.37
+ - CSI_NONLINEAR_BETA=100.0
restart: unless-stopped
diff --git a/test_detection.py b/test_detection.py
new file mode 100644
index 0000000..48b11cd
--- /dev/null
+++ b/test_detection.py
@@ -0,0 +1,202 @@
+#!/usr/bin/env python3
+"""Test isotope detection: vega_best vs vega_final, with/without background subtraction."""
+
+import os
+import sys
+import json
+import numpy as np
+import torch
+from pathlib import Path
+
+# Paths (container-mounted)
+MODELS_DIR = Path(os.environ.get("MODELS_DIR", "/models"))
+DATA_DIR = Path(os.environ.get("DATA_DIR", "/data"))
+VEGA_ML_PATH = Path(os.environ.get("VEGA_ML_PATH", "/models/vega_ml"))
+
+# Add vega_ml to path
+sys.path.insert(0, str(VEGA_ML_PATH))
+from training.vega.model import VegaModel, VegaConfig
+from training.vega.isotope_index import IsotopeIndex
+
+# Energy calibration
+ENERGY_OFFSET = 0.33
+ENERGY_SLOPE = 2.97
+THRESHOLD = 0.5
+
+# CsI(Tl) non-linear response correction
+CSI_NONLINEAR_ALPHA = 0.37
+CSI_NONLINEAR_BETA = 100.0
+
+
+def correct_csilinear_energy(spectrum_rate, num_channels=1023):
+ """Apply inverse CsI(Tl) non-linear response correction.
+
+ Remaps channels so peaks appear at theoretical energy positions
+ (matching training data), correcting for the detector's non-proportional
+ scintillation response that shifts low-energy peaks upward.
+ """
+ alpha = CSI_NONLINEAR_ALPHA
+ beta = CSI_NONLINEAR_BETA
+
+ output_channels = np.arange(num_channels, dtype=np.float64)
+ e_true = ENERGY_OFFSET + ENERGY_SLOPE * output_channels
+
+ # Forward model: E_apparent = E_true * (1 + alpha * exp(-E_true / beta))
+ e_apparent = e_true * (1 + alpha * np.exp(-e_true / beta))
+
+ # Input channel where detector placed counts for this true energy
+ source_channels = (e_apparent - ENERGY_OFFSET) / ENERGY_SLOPE
+ source_channels = np.clip(source_channels, 0, num_channels - 1.001)
+
+ lower = np.floor(source_channels).astype(int)
+ upper = np.minimum(lower + 1, num_channels - 1)
+ frac = source_channels - lower
+
+ corrected = spectrum_rate[lower] * (1 - frac) + spectrum_rate[upper] * frac
+ return corrected
+
+
+def load_model(model_path):
+ """Load a VegaModel checkpoint."""
+ device = torch.device("cpu")
+ checkpoint = torch.load(model_path, map_location=device, weights_only=False)
+ config = VegaConfig(**checkpoint["model_config"])
+ model = VegaModel(config)
+ model.load_state_dict(checkpoint["model_state_dict"])
+ model.eval()
+ return model, config
+
+
+def run_inference(model, config, isotope_index, spectrum_rate, threshold=THRESHOLD,
+ apply_correction=True):
+ """Run inference on a spectrum rate array (1023 channels)."""
+ if spectrum_rate.max() <= 0:
+ return []
+
+ # Apply CsI(Tl) non-linear correction so peaks match training data positions
+ if apply_correction:
+ spectrum_rate = correct_csilinear_energy(spectrum_rate)
+
+ log_spectrum = np.log1p(np.maximum(spectrum_rate, 0))
+ max_val = log_spectrum.max()
+ normalized = log_spectrum / max_val if max_val > 0 else log_spectrum
+ tensor = torch.tensor(normalized, dtype=torch.float32).unsqueeze(0)
+
+ with torch.no_grad():
+ logits, activities = model(tensor)
+
+ probs = torch.sigmoid(logits).numpy()[0]
+ activities = activities.numpy()[0] * config.max_activity_bq
+
+ results = []
+ for i in range(len(probs)):
+ if probs[i] >= threshold:
+ results.append({
+ "isotope": isotope_index.index_to_name(i),
+ "probability": float(probs[i]),
+ "activity_bq": float(activities[i]),
+ })
+ return sorted(results, key=lambda x: -x["probability"])
+
+
+def main():
+ # Load isotope index
+ isotope_index = IsotopeIndex.load(MODELS_DIR / "vega_isotope_index.txt")
+ print(f"Isotope index: {isotope_index.num_isotopes} isotopes\n")
+
+ # Load monitor state (real spectrum from detector)
+ with open(DATA_DIR / "monitor_state.json") as f:
+ state = json.load(f)
+
+ counts = np.array(state["counts"], dtype=np.float64)
+ live_time = state["cumulated_live_time_s"]
+ print(f"Spectre reel : {live_time:.0f}s live time, {counts.sum():.0f} coups, {len(counts)} canaux")
+ print(f"CPS : {state['cps']:.2f}")
+
+ # Load background
+ bg_data = np.load(DATA_DIR / "background_24h.npy", allow_pickle=True).item()
+ bg_counts = bg_data["counts"].astype(np.float64)
+ bg_live_time = float(bg_data["duration"])
+ print(f"Background : {bg_live_time/3600:.1f}h, {bg_counts.sum():.0f} coups\n")
+
+ # Prepare spectra
+ rate = counts[:1023] / live_time
+ bg_rate = bg_counts[:1023] / bg_live_time
+ net_rate = np.clip(rate - bg_rate, 0, None)
+
+ # Apply CsI correction to show peak positions
+ corrected_rate = correct_csilinear_energy(rate)
+ corrected_net = correct_csilinear_energy(net_rate)
+
+ print("=" * 70)
+ print(f" Sans correction CsI:")
+ print(f" Canal max (brut) : {rate.argmax():>4d} ({ENERGY_OFFSET + ENERGY_SLOPE * rate.argmax():.1f} keV)")
+ print(f" Canal max (net) : {net_rate.argmax():>4d} ({ENERGY_OFFSET + ENERGY_SLOPE * net_rate.argmax():.1f} keV)")
+ print(f" Avec correction CsI:")
+ print(f" Canal max (brut) : {corrected_rate.argmax():>4d} ({ENERGY_OFFSET + ENERGY_SLOPE * corrected_rate.argmax():.1f} keV)")
+ print(f" Canal max (net) : {corrected_net.argmax():>4d} ({ENERGY_OFFSET + ENERGY_SLOPE * corrected_net.argmax():.1f} keV)")
+ print(f" Rate max (brut) : {rate.max():.2f} cps")
+ print(f" Rate max (net) : {net_rate.max():.2f} cps")
+ print("=" * 70)
+
+ # Am-241 should be at 59.5 keV → ch ~20
+ print(f"\n Am-241 region (59.5 keV) apres correction CsI:")
+ for ch in range(16, 26):
+ e = ENERGY_OFFSET + ENERGY_SLOPE * ch
+ print(f" ch {ch:3d} ({e:5.1f} keV): brut={corrected_rate[ch]:.5f} net={corrected_net[ch]:.5f}")
+
+ # Load both models
+ models = {
+ "vega_best": load_model(MODELS_DIR / "vega_best.pt"),
+ "vega_final": load_model(MODELS_DIR / "vega_final.pt"),
+ }
+
+ scenarios = {
+ "brut (sans soustraction)": rate,
+ "net (avec soustraction bg)": net_rate,
+ }
+
+ for model_name, (model, config) in models.items():
+ print(f"\n{'─' * 70}")
+ print(f" Modele : {model_name}")
+ print(f"{'─' * 70}")
+ for scenario_name, spectrum in scenarios.items():
+ print(f"\n Scenario : {scenario_name}")
+ results = run_inference(model, config, isotope_index, spectrum)
+ if results:
+ print(f" {'Isotope':>10s} {'Probabilite':>12s} {'Activite (Bq)':>15s}")
+ print(f" {'─'*10} {'─'*12} {'─'*15}")
+ for r in results:
+ print(f" {r['isotope']:>10s} {r['probability']*100:>11.1f}% {r['activity_bq']:>15.1f}")
+ else:
+ print(f" Aucun isotope detecte (seuil = {THRESHOLD})")
+
+ # Also show top-10 probabilities below threshold for context
+ print(f"\n{'═' * 70}")
+ print(" Top-10 probabilites (tous scenarios, meme sous le seuil)")
+ print(f"{'═' * 70}")
+ for model_name, (model, config) in models.items():
+ print(f"\n Modele : {model_name}")
+ for scenario_name, spectrum in scenarios.items():
+ if spectrum.max() <= 0:
+ continue
+ # Apply CsI correction before inference
+ corrected = correct_csilinear_energy(spectrum)
+ log_spectrum = np.log1p(np.maximum(corrected, 0))
+ max_val = log_spectrum.max()
+ normalized = log_spectrum / max_val if max_val > 0 else log_spectrum
+ tensor = torch.tensor(normalized, dtype=torch.float32).unsqueeze(0)
+ with torch.no_grad():
+ logits, _ = model(tensor)
+ probs = torch.sigmoid(logits).numpy()[0]
+ top10 = np.argsort(probs)[::-1][:10]
+ print(f"\n {scenario_name} :")
+ for idx in top10:
+ name = isotope_index.index_to_name(idx)
+ prob = probs[idx]
+ marker = " *" if prob >= THRESHOLD else ""
+ print(f" {name:>10s} : {prob*100:6.2f}%{marker}")
+
+
+if __name__ == "__main__":
+ main()
\ No newline at end of file
diff --git a/train/entrypoint.sh b/train/entrypoint.sh
index 7b2c3e4..ae6d206 100755
--- a/train/entrypoint.sh
+++ b/train/entrypoint.sh
@@ -5,11 +5,11 @@ DATA_DIR="${DATA_DIR:-/data/synthetic}"
MODEL_DIR="${MODEL_DIR:-/models}"
NUM_SAMPLES="${NUM_SAMPLES:-50000}"
EPOCHS="${EPOCHS:-100}"
-BATCH_SIZE="${BATCH_SIZE:-64}"
+BATCH_SIZE="${BATCH_SIZE:-32}"
LEARNING_RATE="${LEARNING_RATE:-0.001}"
DETECTOR="${DETECTOR:-radiacode_103}"
-MIN_DURATION="${MIN_DURATION:-43200}"
-MAX_DURATION="${MAX_DURATION:-86400}"
+MIN_DURATION="${MIN_DURATION:-30}"
+MAX_DURATION="${MAX_DURATION:-300}"
SEED="${SEED:-42}"
MEASURED_BACKGROUND_PATH="${MEASURED_BACKGROUND_PATH:-}"
@@ -20,7 +20,7 @@ echo " Data dir : $DATA_DIR"
echo " Model dir : $MODEL_DIR"
echo " Samples : $NUM_SAMPLES"
echo " Detector : $DETECTOR"
-echo " Duration : $MIN_DURATION-$MAX_DURATION s"
+echo " Duration : $MIN_DURATION-$MAX_DURATION s"
echo " Epochs : $EPOCHS"
echo " Batch size : $BATCH_SIZE"
echo " Learning rate: $LEARNING_RATE"
diff --git a/train/vega_ml/synthetic_spectra/config.py b/train/vega_ml/synthetic_spectra/config.py
index 6a91e2f..f73d392 100644
--- a/train/vega_ml/synthetic_spectra/config.py
+++ b/train/vega_ml/synthetic_spectra/config.py
@@ -3,99 +3,77 @@ Detector Configuration Module
Contains configuration parameters for Radiacode gamma spectrometers
and other detector settings.
+
+Energy calibration matches the real Radiacode 103:
+ E(keV) = 0.33 + 2.97 * channel_index
+Uses 1023 channels (channel 1023 is overflow, excluded).
"""
-from dataclasses import dataclass, field
-from typing import Dict, Optional
+from dataclasses import dataclass
+from typing import Dict
import numpy as np
@dataclass
class DetectorConfig:
"""Configuration for a gamma spectrometer detector."""
-
- name: str
- # Energy range in keV
- energy_min_kev: float = 20.0
- energy_max_kev: float = 3000.0
-
- # Number of channels
- num_channels: int = 1024
- # Some devices/software workflows treat channel 0 as unreliable/noisy.
- # This project models "usable" channels by skipping the first raw channel.
- skip_first_channel: bool = True
-
+ name: str
+ # Energy calibration: E = calibration_offset + calibration_slope * channel
+ # Must match the real detector calibration used in inference.
+ calibration_offset_kev: float = 0.33
+ calibration_slope_kev: float = 2.97
+
+ # Number of usable channels (1023 for Radiacode, channel 1023 is overflow)
+ num_channels: int = 1023
+
# FWHM at 662 keV (Cs-137 reference) as fraction
fwhm_at_662: float = 0.084 # 8.4%
fwhm_uncertainty: float = 0.003 # ±0.3%
-
+
# Detector crystal type
crystal_type: str = "CsI(Tl)"
-
+
# Sensitivity: counts per second at 1 μSv/h for Cs-137
sensitivity_cps_per_usvh: float = 30.0
-
+
# Detector volume in cm³
detector_volume_cm3: float = 1.0
-
- def get_channel_width_kev(self) -> float:
- """Get the width of each channel in keV."""
- return (self.energy_max_kev - self.energy_min_kev) / self.num_channels
-
- def get_energy_bins(self) -> np.ndarray:
- """Get array of energy bin centers (keV) for the modeled usable channels."""
- channel_width = self.get_channel_width_kev()
- # Raw device channels are assumed to be 0..num_channels-1 with centers:
- # E_center(k) = E_min + (k + 0.5) * channel_width
- # If we skip the first raw channel (k=0), we model usable channels k=1..num_channels-1.
- start_raw_channel = 1 if self.skip_first_channel else 0
- raw_channels = np.arange(start_raw_channel, self.num_channels, dtype=np.float64)
- return self.energy_min_kev + (raw_channels + 0.5) * channel_width
-
+ def get_energy_bins(self) -> np.ndarray:
+ """Get array of energy bin centers (keV) matching the real detector calibration."""
+ channels = np.arange(self.num_channels, dtype=np.float64)
+ return self.calibration_offset_kev + self.calibration_slope_kev * channels
+
def get_fwhm_at_energy(self, energy_kev: float) -> float:
"""
Calculate FWHM at a given energy.
-
+
For scintillators, FWHM scales approximately as sqrt(E).
- FWHM(E) = FWHM_662 * sqrt(662/E) * E / 662 = FWHM_662 * sqrt(E/662)
+ FWHM(E) = FWHM_662 * sqrt(E/662)
"""
- return self.fwhm_at_662 * np.sqrt(662.0 / energy_kev) * energy_kev
-
+ return self.fwhm_at_662 * np.sqrt(energy_kev / 662.0) * 662.0
+
def get_sigma_at_energy(self, energy_kev: float) -> float:
- """
- Get Gaussian sigma at a given energy.
- sigma = FWHM / (2 * sqrt(2 * ln(2))) ≈ FWHM / 2.355
- """
+ """Get Gaussian sigma at a given energy."""
fwhm = self.get_fwhm_at_energy(energy_kev)
return fwhm / 2.355
-
+
def energy_to_channel(self, energy_kev: float) -> int:
- """Convert energy in keV to modeled usable channel index."""
- channel_width = self.get_channel_width_kev()
- raw_channel = int((energy_kev - self.energy_min_kev) / channel_width)
- if self.skip_first_channel:
- channel = raw_channel - 1
- max_channel = self.num_channels - 2
- else:
- channel = raw_channel
- max_channel = self.num_channels - 1
- return max(0, min(max_channel, channel))
+ """Convert energy in keV to channel index."""
+ channel = int((energy_kev - self.calibration_offset_kev) / self.calibration_slope_kev)
+ return max(0, min(self.num_channels - 1, channel))
def channel_to_energy(self, channel: int) -> float:
- """Convert modeled usable channel index to energy bin center (keV)."""
- channel_width = self.get_channel_width_kev()
- raw_channel = channel + (1 if self.skip_first_channel else 0)
- raw_channel = max(0, min(self.num_channels - 1, int(raw_channel)))
- return self.energy_min_kev + (raw_channel + 0.5) * channel_width
+ """Convert channel index to energy in keV."""
+ return self.calibration_offset_kev + self.calibration_slope_kev * channel
# Pre-defined configurations for Radiacode devices
RADIACODE_CONFIGS: Dict[str, DetectorConfig] = {
"radiacode_101": DetectorConfig(
name="Radiacode 101",
- fwhm_at_662=0.095, # 9.5% (original model, similar to 102)
+ fwhm_at_662=0.095, # 9.5%
fwhm_uncertainty=0.004,
crystal_type="CsI(Tl)",
sensitivity_cps_per_usvh=30.0,
@@ -119,8 +97,7 @@ RADIACODE_CONFIGS: Dict[str, DetectorConfig] = {
),
"radiacode_103g": DetectorConfig(
name="Radiacode 103G",
- energy_min_kev=25.0, # Tech spec lists 0.025…3 MeV
- fwhm_at_662=0.074, # 7.4% (GAGG crystal - better resolution)
+ fwhm_at_662=0.074, # 7.4% (GAGG crystal)
fwhm_uncertainty=0.003,
crystal_type="GAGG(Ce)",
sensitivity_cps_per_usvh=40.0,
@@ -131,12 +108,12 @@ RADIACODE_CONFIGS: Dict[str, DetectorConfig] = {
fwhm_at_662=0.084, # 8.4%
fwhm_uncertainty=0.003,
crystal_type="CsI(Tl)",
- sensitivity_cps_per_usvh=77.0, # Higher sensitivity
- detector_volume_cm3=2.5, # Larger crystal
+ sensitivity_cps_per_usvh=77.0,
+ detector_volume_cm3=2.5,
),
}
def get_default_config() -> DetectorConfig:
"""Get the default detector configuration (Radiacode 103)."""
- return RADIACODE_CONFIGS["radiacode_103"]
+ return RADIACODE_CONFIGS["radiacode_103"]
\ No newline at end of file
diff --git a/train/vega_ml/synthetic_spectra/generate_spectra.py b/train/vega_ml/synthetic_spectra/generate_spectra.py
index 8df1215..b5a6c79 100644
--- a/train/vega_ml/synthetic_spectra/generate_spectra.py
+++ b/train/vega_ml/synthetic_spectra/generate_spectra.py
@@ -128,19 +128,21 @@ def generate_training_batch(
num_samples: int,
output_dir: Path,
detector_name: str = "radiacode_103",
- duration_range: tuple = (60, 300),
+ duration_range: tuple = (30, 300),
activity_range: tuple = (1.0, 100.0),
- single_isotope_fraction: float = 0.4,
- dual_isotope_fraction: float = 0.3,
- multi_isotope_fraction: float = 0.2,
+ single_isotope_fraction: float = 0.3,
+ dual_isotope_fraction: float = 0.2,
+ multi_isotope_fraction: float = 0.15,
background_only_fraction: float = 0.1,
+ low_signal_fraction: float = 0.15,
+ subtracted_fraction: float = 0.1,
save_png: bool = False,
random_seed: int = None,
measured_background_path: str = None,
) -> list:
"""
Generate a batch of training samples with various configurations.
-
+
Args:
num_samples: Total number of samples to generate
output_dir: Output directory for spectra and labels
@@ -151,9 +153,11 @@ def generate_training_batch(
dual_isotope_fraction: Fraction of two-isotope samples
multi_isotope_fraction: Fraction of 3+ isotope samples
background_only_fraction: Fraction of background-only samples
+ low_signal_fraction: Fraction of low-activity samples (0.01-5 Bq)
+ subtracted_fraction: Fraction of background-subtracted samples
save_png: Whether to also save PNG images
random_seed: Random seed for reproducibility
-
+
Returns:
List of generated spectra
"""
@@ -181,11 +185,13 @@ def generate_training_batch(
n_dual = int(num_samples * dual_isotope_fraction)
n_multi = int(num_samples * multi_isotope_fraction)
n_background = int(num_samples * background_only_fraction)
-
+ n_low_signal = int(num_samples * low_signal_fraction)
+ n_subtracted = int(num_samples * subtracted_fraction)
+
# Adjust to ensure we hit exactly num_samples
- remaining = num_samples - (n_single + n_dual + n_multi + n_background)
+ remaining = num_samples - (n_single + n_dual + n_multi + n_background + n_low_signal + n_subtracted)
n_single += remaining
-
+
total_generated = 0
print(f"\nGenerating {num_samples} synthetic spectra:")
@@ -193,6 +199,8 @@ def generate_training_batch(
print(f" - Dual isotope: {n_dual}")
print(f" - Multi isotope (3+): {n_multi}")
print(f" - Background only: {n_background}")
+ print(f" - Low signal (0.01-5 Bq): {n_low_signal}")
+ print(f" - Background-subtracted: {n_subtracted}")
print()
sample_num = 0
@@ -314,6 +322,77 @@ def generate_training_batch(
sample_num += 1
+ # Generate low-signal samples (weak sources, 0.01-5 Bq)
+ print("Generating low-signal samples...")
+ for i in range(n_low_signal):
+ isotope = np.random.choice(isotope_pool)
+ activity = np.random.uniform(0.01, 5.0)
+ duration = np.random.uniform(*duration_range)
+
+ spectrum = generate_single_isotope_sample(
+ generator,
+ isotope,
+ activity,
+ duration,
+ detector_name=detector_name,
+ include_background=True,
+ measured_background_path=measured_background_path,
+ )
+
+ save_spectrum(
+ spectrum,
+ spectra_dir,
+ save_image=True,
+ image_format='npy'
+ )
+ del spectrum
+
+ sample_num += 1
+
+ if sample_num % 100 == 0:
+ print(f" Generated {sample_num}/{num_samples} samples...")
+
+ # Generate background-subtracted samples (simulates inference pipeline)
+ print("Generating background-subtracted samples...")
+ for i in range(n_subtracted):
+ num_iso = np.random.choice([1, 2, 3], p=[0.5, 0.3, 0.2])
+ isotopes = np.random.choice(isotope_pool, size=num_iso, replace=False)
+ activities = [np.random.uniform(0.1, 50.0) for _ in range(num_iso)]
+ duration = np.random.uniform(*duration_range)
+
+ sources = [
+ IsotopeSource(
+ isotope_name=name,
+ activity_bq=activity,
+ include_daughters=True
+ )
+ for name, activity in zip(isotopes, activities)
+ ]
+
+ config = SpectrumConfig(
+ duration_seconds=duration,
+ sources=sources,
+ include_background=True,
+ subtract_background=True,
+ detector_name=detector_name,
+ measured_background_path=measured_background_path,
+ )
+
+ spectrum = generator.generate_spectrum(config)
+
+ save_spectrum(
+ spectrum,
+ spectra_dir,
+ save_image=True,
+ image_format='npy'
+ )
+ del spectrum
+
+ sample_num += 1
+
+ if sample_num % 100 == 0:
+ print(f" Generated {sample_num}/{num_samples} samples...")
+
total_generated = sample_num
print(f"\nGenerated {total_generated} samples total")
diff --git a/train/vega_ml/synthetic_spectra/generator.py b/train/vega_ml/synthetic_spectra/generator.py
index 62170e0..1f54300 100644
--- a/train/vega_ml/synthetic_spectra/generator.py
+++ b/train/vega_ml/synthetic_spectra/generator.py
@@ -49,14 +49,14 @@ class IsotopeSource:
@dataclass
class SpectrumConfig:
"""Configuration for a single spectrum generation."""
-
+
# Time parameters
duration_seconds: float = 60.0
time_interval_seconds: float = 1.0 # Each row in the spectrogram
-
+
# Sources to include
sources: List[IsotopeSource] = field(default_factory=list)
-
+
# Background options
include_background: bool = True
background_cps: float = 5.0
@@ -64,18 +64,25 @@ class SpectrumConfig:
include_radon: bool = True
include_thorium: bool = True
measured_background_path: Optional[str] = None
-
+
+ # Background subtraction simulation
+ # When True, generates a second independent background realization
+ # and subtracts it from the spectrum, then clips negatives to 0.
+ # This simulates what happens at inference time (measured bg subtraction).
+ subtract_background: bool = False
+
# Detector configuration
detector_name: str = "radiacode_103"
-
+
# Noise options
apply_poisson: bool = True
apply_electronic: bool = False
electronic_noise_sigma: float = 0.5
-
- # Normalization
+
+ # Normalization — "log1p" preserves relative signal levels,
+ # works well after background subtraction where many channels are ~0.
normalize: bool = True
- normalization_method: str = "max" # max, sum, log, sqrt
+ normalization_method: str = "log1p" # max, sum, log, sqrt, log1p
@dataclass
@@ -272,7 +279,7 @@ class SpectrumGenerator:
all_source_isotopes.extend(src_iso)
all_background_isotopes.extend(bg_iso)
- # Apply noise
+ # Apply noise before any subtraction (Poisson noise on raw counts)
if config.apply_poisson:
spectrum = apply_poisson_noise(spectrum)
@@ -282,6 +289,24 @@ class SpectrumGenerator:
config.electronic_noise_sigma
)
+ # Simulate background subtraction (matches inference pipeline)
+ if config.subtract_background and config.include_background:
+ # Generate an independent background realization
+ bg_spectrum2, _ = generate_environmental_background(
+ self.energy_bins,
+ config.duration_seconds,
+ background_cps=config.background_cps,
+ include_k40=config.include_k40,
+ include_radon=config.include_radon,
+ include_thorium=config.include_thorium,
+ detector_config=self.detector_config,
+ measured_background_path=config.measured_background_path,
+ )
+ if config.apply_poisson:
+ bg_spectrum2 = apply_poisson_noise(bg_spectrum2)
+ # Subtract and clip — same as inference: net = clip(rate - bg_rate, 0, inf)
+ spectrum = np.maximum(spectrum - bg_spectrum2, 0)
+
# Normalize if requested
if config.normalize:
spectrum = normalize_spectrum(spectrum, config.normalization_method)
diff --git a/train/vega_ml/synthetic_spectra/physics/spectrum_physics.py b/train/vega_ml/synthetic_spectra/physics/spectrum_physics.py
index 5cb6a2a..231a33e 100644
--- a/train/vega_ml/synthetic_spectra/physics/spectrum_physics.py
+++ b/train/vega_ml/synthetic_spectra/physics/spectrum_physics.py
@@ -184,38 +184,148 @@ def calculate_expected_counts(
return expected
+def _k_escape_fraction(energy_kev: float, detector_config: Optional[DetectorConfig] = None) -> float:
+ """
+ Calculate K-escape peak fraction for CsI(Tl) detector.
+
+ For iodine K-shell (binding energy ~33.2 keV), when a gamma photon
+ interacts with the K-shell, there's a chance the K X-ray escapes the
+ crystal, producing a peak at E - E_Ka (~28.5 keV for I K-alpha).
+
+ The escape fraction decreases with energy as the photoelectric cross-section
+ ratio (K-shell / total) decreases.
+
+ Args:
+ energy_kev: Gamma energy in keV
+ detector_config: Detector configuration
+
+ Returns:
+ Fraction of photopeak counts that appear in the K-escape peak
+ """
+ if energy_kev <= 33.2:
+ return 0.0
+
+ # K-shell binding energy for iodine
+ k_binding = 33.2 # keV
+
+ # K-escape fraction for CsI(Tl) detector
+ # Based on measured data: ~35% at 60 keV, ~15% at 150 keV, ~5% at 662 keV
+ # Model as: fraction = A * (1 - exp(-E/B)) where A and B are fit parameters
+ # Fitted to typical CsI K-escape measurements
+ fraction = 0.40 * (1.0 - np.exp(-(energy_kev - k_binding) / 80.0))
+
+ return float(np.clip(fraction, 0.0, 0.45))
+
+
+def _asymmetric_peak(
+ energy_bins: np.ndarray,
+ peak_energy: float,
+ sigma: float,
+ amplitude: float,
+ tail_fraction: float = 0.0,
+ tail_sigma_ratio: float = 3.0
+) -> np.ndarray:
+ """
+ Generate an asymmetric peak using an exponentially-modified Gaussian.
+
+ For scintillation detectors at low energies, incomplete charge collection
+ creates a low-energy tail. The tail fraction increases at lower energies.
+
+ Args:
+ energy_bins: Array of energy bin centers (keV)
+ peak_energy: Center energy of peak (keV)
+ sigma: Gaussian sigma (keV)
+ amplitude: Total peak area (counts)
+ tail_fraction: Fraction of peak area in low-energy tail (0-0.5)
+ tail_sigma_ratio: Ratio of tail sigma to peak sigma
+
+ Returns:
+ Array of counts in each bin
+ """
+ # Main Gaussian component
+ main_peak = gaussian_peak(energy_bins, peak_energy, sigma, amplitude * (1 - tail_fraction))
+
+ if tail_fraction <= 0:
+ return main_peak
+
+ # Low-energy tail: Gaussian shifted to lower energy with broader width
+ tail_sigma = sigma * tail_sigma_ratio
+ tail_energy = peak_energy - 2.0 * sigma # Tail centered 2 sigma below peak
+ tail_peak = gaussian_peak(energy_bins, tail_energy, tail_sigma, amplitude * tail_fraction)
+
+ return main_peak + tail_peak
+
+
def generate_peak_spectrum(
energy_bins: np.ndarray,
peak_params: PeakParameters,
detector_config: Optional[DetectorConfig] = None
) -> np.ndarray:
"""
- Generate a single gamma peak with detector response.
-
+ Generate a single gamma peak with realistic CsI(Tl) detector response.
+
+ Includes:
+ - Asymmetric peak shape (low-energy tail from incomplete charge collection)
+ - K-escape peak (Iodine K-shell X-ray escape at E - 28.5 keV)
+ - Energy-dependent resolution
+
+ Note: Peaks are placed at theoretical gamma energies. The non-linear
+ CsI(Tl) response correction is applied in the inference pipeline
+ (radiacode_monitor.py), not here, to keep training data detector-independent.
+
Args:
- energy_bins: Array of energy bin centers (keV)
+ energy_bins: Array of energy bin centers (keV) matching detector calibration
peak_params: Peak parameters
detector_config: Detector configuration
-
+
Returns:
Array of expected counts in each bin (not yet Poisson sampled)
"""
if detector_config is None:
detector_config = get_default_config()
-
+
# Calculate expected counts
- amplitude = calculate_expected_counts(peak_params, detector_config)
-
- if amplitude <= 0:
+ total_amplitude = calculate_expected_counts(peak_params, detector_config)
+
+ if total_amplitude <= 0:
return np.zeros_like(energy_bins)
-
+
# Calculate peak width
fwhm_kev = calculate_fwhm(peak_params.energy_kev, detector_config.fwhm_at_662)
sigma = fwhm_to_sigma(fwhm_kev)
-
- # Generate Gaussian peak
- peak = gaussian_peak(energy_bins, peak_params.energy_kev, sigma, amplitude)
-
+
+ # Low-energy tail fraction: increases at lower energies due to
+ # incomplete charge collection in CsI(Tl)
+ if peak_params.energy_kev < 200:
+ tail_frac = 0.15 * (1.0 - peak_params.energy_kev / 200.0)
+ else:
+ tail_frac = 0.0
+
+ # Generate main peak (asymmetric)
+ peak = _asymmetric_peak(
+ energy_bins, peak_params.energy_kev, sigma,
+ total_amplitude, tail_fraction=tail_frac
+ )
+
+ # K-escape peak for CsI(Tl)
+ escape_frac = _k_escape_fraction(peak_params.energy_kev, detector_config)
+ if escape_frac > 0:
+ escape_energy = peak_params.energy_kev - 28.5 # I K-alpha at 28.5 keV
+ if escape_energy > 20: # Only if above detection threshold
+ escape_amplitude = total_amplitude * escape_frac
+ # Reduce main peak amplitude
+ peak = peak * (1 - escape_frac)
+
+ # Escape peak has slightly broader resolution
+ escape_fwhm = calculate_fwhm(escape_energy, detector_config.fwhm_at_662)
+ escape_sigma = fwhm_to_sigma(escape_fwhm) * 1.3
+
+ escape_peak = _asymmetric_peak(
+ energy_bins, escape_energy, escape_sigma,
+ escape_amplitude, tail_fraction=0.25
+ )
+ peak = peak + escape_peak
+
return peak
@@ -636,11 +746,11 @@ def apply_electronic_noise(
def normalize_spectrum(
spectrum: np.ndarray,
- method: str = "max"
+ method: str = "log1p"
) -> np.ndarray:
"""
Normalize a spectrum for ML training.
-
+
Args:
spectrum: Raw count spectrum
method: Normalization method
@@ -648,7 +758,8 @@ def normalize_spectrum(
- "sum": Divide by total counts (probability distribution)
- "log": Log transform then max normalize
- "sqrt": Square root transform then max normalize
-
+ - "log1p": log(1+x) then max normalize (best for bg-subtracted spectra)
+
Returns:
Normalized spectrum
"""
@@ -657,7 +768,7 @@ def normalize_spectrum(
if max_val > 0:
return spectrum / max_val
return spectrum
-
+
elif method == "sum":
total = spectrum.sum()
if total > 0:
@@ -678,6 +789,13 @@ def normalize_spectrum(
if max_val > 0:
return sqrt_spec / max_val
return sqrt_spec
-
+
+ elif method == "log1p":
+ log_spec = np.log1p(np.maximum(spectrum, 0))
+ max_val = log_spec.max()
+ if max_val > 0:
+ return log_spec / max_val
+ return log_spec
+
else:
raise ValueError(f"Unknown normalization method: {method}")
diff --git a/train/vega_ml/training/vega/__init__.py b/train/vega_ml/training/vega/__init__.py
index 93de62c..389a9d8 100644
--- a/train/vega_ml/training/vega/__init__.py
+++ b/train/vega_ml/training/vega/__init__.py
@@ -14,7 +14,12 @@ Features:
from .model import VegaModel, VegaConfig
from .dataset import SpectrumDataset, create_data_loaders
-from .train import train_vega, VegaTrainer
+
+def __getattr__(name):
+ if name in ('train_vega', 'VegaTrainer'):
+ from .train import train_vega, VegaTrainer
+ return locals()[name]
+ raise AttributeError(f"module {__name__!r} has no attribute {name!r}")
__all__ = [
'VegaModel',
diff --git a/train/vega_ml/training/vega/dataset.py b/train/vega_ml/training/vega/dataset.py
index 6b2cf66..c6303a2 100644
--- a/train/vega_ml/training/vega/dataset.py
+++ b/train/vega_ml/training/vega/dataset.py
@@ -31,24 +31,38 @@ class SpectrumSample:
detector: str
+def normalize_log1p(spectrum: np.ndarray) -> np.ndarray:
+ """Log1p normalization: log(1 + x) / max(log(1 + x)).
+
+ Preserves relative signal levels across channels, works well when
+ many channels are zero (e.g. after background subtraction).
+ """
+ log_spec = np.log1p(np.maximum(spectrum, 0))
+ max_val = log_spec.max()
+ if max_val > 0:
+ return log_spec / max_val
+ return log_spec
+
+
class SpectrumDataset(Dataset):
"""
PyTorch Dataset for synthetic gamma spectra.
-
+
Loads spectra from numpy files and their labels from JSON files.
Supports both individual JSON files per sample (efficient for large datasets)
and combined labels.json (legacy format).
-
+
Converts to tensors suitable for the Vega model.
"""
-
+
def __init__(
self,
data_dir: Path,
isotope_index: Optional[IsotopeIndex] = None,
max_activity_bq: float = 1000.0,
collapse_time: bool = True,
- transform=None
+ transform=None,
+ normalization: str = "log1p"
):
"""
Initialize the dataset.
@@ -66,6 +80,7 @@ class SpectrumDataset(Dataset):
self.max_activity_bq = max_activity_bq
self.collapse_time = collapse_time
self.transform = transform
+ self.normalization = normalization
# Detect label format and load sample list
self.use_individual_labels = self._detect_label_format()
@@ -156,7 +171,15 @@ class SpectrumDataset(Dataset):
if self.collapse_time and spectrum.ndim == 2:
# Average across time intervals to get single spectrum
spectrum = spectrum.mean(axis=0)
-
+
+ # Normalize spectrum
+ if self.normalization == "log1p":
+ spectrum = normalize_log1p(spectrum)
+ elif self.normalization == "max":
+ max_val = spectrum.max()
+ if max_val > 0:
+ spectrum = spectrum / max_val
+
# Convert to tensor
spectrum_tensor = torch.tensor(spectrum, dtype=torch.float32)
diff --git a/web/app/bg_calibration.py b/web/app/bg_calibration.py
index f4cde7b..7fc3367 100644
--- a/web/app/bg_calibration.py
+++ b/web/app/bg_calibration.py
@@ -4,19 +4,15 @@ 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.
+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.interpolate import make_smoothing_spline
+from scipy.ndimage import gaussian_filter1d
from scipy.signal import savgol_filter
from app.config import ENERGY_OFFSET, ENERGY_SLOPE, NUM_CHANNELS
@@ -31,77 +27,76 @@ 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 _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, smoothed_cps):
- """Iteratively estimate and subtract photopeak contributions."""
- continuum = smoothed_cps.copy()
- peak_amplitudes = []
+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:
- 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
+ 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:
- 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))
+ gaussian = peak_height * np.exp(-0.5 * ((channels - idx) / sig) ** 2)
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
+ return np.maximum(continuum, 0), [{"energy_keV": e, "amplitude": 0.0} for e, _ in PHOTOPEAK_LINES]
def calibrate_spline(measured_cps, energy_axis):
"""
- Fit a smoothing spline to the peak-subtracted continuum.
+ Fit continuum using peak subtraction + Gaussian smoothing.
- 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.
+ 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.
"""
- E = energy_axis
- y_smooth = _smooth(measured_cps)
- continuum, peak_amplitudes = _subtract_peaks(E, y_smooth)
+ # 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)
- # Ensure positive values for spline fitting
- continuum = np.maximum(continuum, 0)
+ # Step 2: Subtract known photopeaks
+ continuum, peak_amplitudes = _subtract_peaks(energy_axis, y_smooth)
- # 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)}
+ # 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 - fit_cps
+ 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": fit_cps,
+ "continuum_cps": continuum_smooth,
"peak_amplitudes": peak_amplitudes,
"r_squared": float(r_squared),
"residuals_rms": float(np.sqrt(np.mean(residuals ** 2))),
@@ -109,29 +104,24 @@ def calibrate_spline(measured_cps, energy_axis):
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.
- """
+ """Fit the continuum model using peak subtraction + Gaussian smoothing."""
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
+ "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",
+ "method": "peak_subtract_gaussian",
}
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()
@@ -151,13 +141,20 @@ def load_or_calibrate():
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_result = json.load(f)
- return _cached_result
+ cached = json.load(f)
+ # Invalidate if method changed
+ if cached.get("method") != "peak_subtract_gaussian":
+ cached = None
except Exception:
- pass
+ cached = None
+
+ if cached and "continuum_cps" in cached:
+ _cached_result = cached
+ return _cached_result
from app.config import BACKGROUND_PATH, BACKGROUND_SNAPSHOT_PATH
@@ -199,10 +196,14 @@ def load_or_calibrate():
"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)
+ # 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
\ No newline at end of file
diff --git a/web/app/config.py b/web/app/config.py
index 30eb039..197be85 100644
--- a/web/app/config.py
+++ b/web/app/config.py
@@ -1,4 +1,5 @@
import os
+import numpy as np
from pathlib import Path
STATE_PATH = Path(os.environ.get("STATE_PATH", "/data/monitor_state.json"))
@@ -11,8 +12,47 @@ ISOTOPE_INDEX_PATH = Path(os.environ.get("ISOTOPE_INDEX_PATH", "/models/vega_iso
ENERGY_OFFSET = float(os.environ.get("ENERGY_CALIBRATION_OFFSET", "0.33"))
ENERGY_SLOPE = float(os.environ.get("ENERGY_CALIBRATION_SLOPE", "2.97"))
NUM_CHANNELS = 1023 # Last channel (1023) is overflow bin, excluded from display
+ENERGY_MIN = 30.0 # keV - detector lower limit
+ENERGY_MAX = 3000.0 # keV - detector upper limit (3 MeV)
+
+# CsI(Tl) non-linear response correction parameters
+# Matches the detector's non-proportional scintillation response
+CSI_NONLINEAR_ALPHA = float(os.environ.get("CSI_NONLINEAR_ALPHA", "0.37"))
+CSI_NONLINEAR_BETA = float(os.environ.get("CSI_NONLINEAR_BETA", "100.0"))
+
+
+def correct_csi_nonlinear(spectrum, num_channels=1023):
+ """Apply inverse CsI(Tl) non-linear response correction.
+
+ Remaps spectrum channels so peaks appear at their theoretical energy
+ positions, correcting for the detector's non-proportional scintillation
+ response that shifts low-energy peaks upward.
+ """
+ alpha = CSI_NONLINEAR_ALPHA
+ beta = CSI_NONLINEAR_BETA
+ output_channels = np.arange(num_channels, dtype=np.float64)
+ e_true = ENERGY_OFFSET + ENERGY_SLOPE * output_channels
+ e_apparent = e_true * (1 + alpha * np.exp(-e_true / beta))
+ source_channels = (e_apparent - ENERGY_OFFSET) / ENERGY_SLOPE
+ source_channels = np.clip(source_channels, 0, num_channels - 1.001)
+ lower = np.floor(source_channels).astype(int)
+ upper = np.minimum(lower + 1, num_channels - 1)
+ frac = source_channels - lower
+ return spectrum[lower] * (1 - frac) + spectrum[upper] * frac
def energy_axis():
- """Generate energy axis in keV from channel numbers."""
- return [round(ENERGY_OFFSET + ENERGY_SLOPE * i, 2) for i in range(NUM_CHANNELS)]
\ No newline at end of file
+ """Generate energy axis in keV from channel numbers, clipped to detector range."""
+ axis = [round(ENERGY_OFFSET + ENERGY_SLOPE * i, 2) for i in range(NUM_CHANNELS)]
+ return [e for e in axis if ENERGY_MIN <= e <= ENERGY_MAX]
+
+
+def energy_mask():
+ """Return boolean mask of channels within detector energy range."""
+ return [ENERGY_MIN <= ENERGY_OFFSET + ENERGY_SLOPE * i <= ENERGY_MAX for i in range(NUM_CHANNELS)]
+
+
+def clip_to_range(arr):
+ """Clip array to detector energy range using energy mask."""
+ mask = energy_mask()
+ return [arr[i] for i in range(len(arr)) if mask[i]]
\ No newline at end of file
diff --git a/web/app/routers/background.py b/web/app/routers/background.py
index b68c2ca..5f3f4e8 100644
--- a/web/app/routers/background.py
+++ b/web/app/routers/background.py
@@ -1,6 +1,6 @@
import json
from fastapi import APIRouter, HTTPException
-from app.config import BACKGROUND_SNAPSHOT_PATH, BACKGROUND_PATH, energy_axis, NUM_CHANNELS, ENERGY_OFFSET, ENERGY_SLOPE
+from app.config import BACKGROUND_SNAPSHOT_PATH, BACKGROUND_PATH, energy_axis, NUM_CHANNELS, ENERGY_OFFSET, ENERGY_SLOPE, clip_to_range
from app.theoretical_bg import generate_continuum_only
from app.noise import extract_continuum
import numpy as np
@@ -25,8 +25,9 @@ def _load_reference():
return None
try:
bg_data = np.load(str(BACKGROUND_PATH), allow_pickle=True).item()
+ raw_counts = [round(float(c), 1) for c in bg_data["counts"][:NUM_CHANNELS]]
return {
- "counts": [round(float(c), 1) for c in bg_data["counts"][:NUM_CHANNELS]],
+ "counts": clip_to_range(raw_counts),
"live_time_s": round(float(bg_data["duration"]), 1),
}
except Exception:
@@ -54,11 +55,12 @@ async def get_background_spectrum():
"""Live background spectrum (from snapshot) with energy axis."""
snapshot = _load_snapshot()
live_time = snapshot.get("live_time_s", 0)
+ raw_spectrum = snapshot.get("spectrum", [0] * 1024)[:NUM_CHANNELS]
return {
- "channels": list(range(NUM_CHANNELS)),
+ "channels": clip_to_range(list(range(NUM_CHANNELS))),
"energy_kev": energy_axis(),
- "counts": snapshot.get("spectrum", [0] * 1024)[:NUM_CHANNELS],
+ "counts": clip_to_range(raw_spectrum),
"live_time_s": live_time,
"cps": snapshot.get("cps", 0),
"top_peaks": snapshot.get("top_peaks", []),
@@ -74,7 +76,7 @@ async def get_background_reference():
raise HTTPException(status_code=404, detail="No 24h reference background available")
return {
- "channels": list(range(NUM_CHANNELS)),
+ "channels": clip_to_range(list(range(NUM_CHANNELS))),
"energy_kev": energy_axis(),
"counts": ref["counts"],
"live_time_s": ref["live_time_s"],
@@ -84,7 +86,10 @@ async def get_background_reference():
@router.get("/continuum")
async def get_continuum(cps: float = 6.0, live_time_s: float = 3600.0):
"""CsI(Tl) detector response continuum only (no photopeaks, no noise)."""
- return generate_continuum_only(cps=cps, live_time_s=live_time_s)
+ raw = generate_continuum_only(cps=cps, live_time_s=live_time_s)
+ raw["energy_kev"] = clip_to_range(raw["energy_kev"])
+ raw["counts"] = clip_to_range(raw["counts"])
+ return raw
@router.get("/fit")
@@ -132,10 +137,14 @@ async def fit_background():
# Build fitted curve at same scale as measured
fitted_counts = build_calibrated_continuum(e_axis, measured_counts.sum(), result)
+ e_list = [round(float(E), 2) for E in e_axis]
+ m_list = [round(float(c), 1) for c in measured_counts]
+ f_list = [round(float(c), 1) for c in fitted_counts]
+
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],
+ "energy_kev": clip_to_range(e_list),
+ "measured_counts": clip_to_range(m_list),
+ "fitted_counts": clip_to_range(f_list),
"method": result.get("method", "spline"),
"r_squared": result["r_squared"],
"residuals_rms": result["residuals_rms"],
@@ -174,7 +183,9 @@ async def get_background_noise():
e_axis = ENERGY_OFFSET + ENERGY_SLOPE * channels
continuum = extract_continuum(counts, energy_axis=e_axis)
+ e_list = [round(float(E), 2) for E in e_axis]
+ c_list = [round(float(c), 1) for c in continuum]
return {
- "energy_kev": [round(float(E), 2) for E in e_axis],
- "counts": [round(float(c), 1) for c in continuum],
+ "energy_kev": clip_to_range(e_list),
+ "counts": clip_to_range(c_list),
}
\ No newline at end of file
diff --git a/web/app/routers/spectrum.py b/web/app/routers/spectrum.py
index 1026d30..7c5aaba 100644
--- a/web/app/routers/spectrum.py
+++ b/web/app/routers/spectrum.py
@@ -1,6 +1,7 @@
import json
from fastapi import APIRouter, HTTPException
-from app.config import STATE_PATH, BACKGROUND_PATH, energy_axis, NUM_CHANNELS
+from app.config import (STATE_PATH, BACKGROUND_PATH, energy_axis, NUM_CHANNELS,
+ clip_to_range, correct_csi_nonlinear)
import numpy as np
router = APIRouter()
@@ -8,7 +9,7 @@ router = APIRouter()
@router.get("/current")
async def get_current_spectrum():
- """Current accumulated spectrum with energy axis."""
+ """Current accumulated spectrum with energy axis (CsI-corrected)."""
if not STATE_PATH.exists():
raise HTTPException(status_code=503, detail="Monitor not started yet")
@@ -18,6 +19,9 @@ async def get_current_spectrum():
except (json.JSONDecodeError, OSError):
raise HTTPException(status_code=503, detail="Monitor state file corrupt")
+ raw_counts = state.get("counts", [0] * 1024)[:NUM_CHANNELS]
+ # Apply CsI correction so peaks appear at theoretical energy positions
+ corrected_counts = correct_csi_nonlinear(np.array(raw_counts, dtype=np.float64))
return {
"timestamp": state.get("timestamp", ""),
"connected": state.get("connected", False),
@@ -27,15 +31,15 @@ async def get_current_spectrum():
"total_counts": state.get("total_counts", 0),
"background_subtracted": state.get("background_subtracted", False),
"isotopes_detected": state.get("isotopes_detected", []),
- "channels": list(range(NUM_CHANNELS)),
+ "channels": clip_to_range(list(range(NUM_CHANNELS))),
"energy_kev": energy_axis(),
- "counts": state.get("counts", [0] * 1024)[:NUM_CHANNELS],
+ "counts": clip_to_range([round(float(c), 1) for c in corrected_counts]),
}
@router.get("/difference")
async def get_difference_spectrum():
- """Background-subtracted spectrum (net signal)."""
+ """Background-subtracted spectrum (net signal, CsI-corrected)."""
if not STATE_PATH.exists():
raise HTTPException(status_code=503, detail="Monitor not started yet")
@@ -59,18 +63,21 @@ async def get_difference_spectrum():
bg_live_time = float(bg_data["duration"])
bg_rate = bg_counts / bg_live_time
net_rate = np.clip(rate - bg_rate, 0, None)
- net_counts = net_rate * live_time
+ # Apply CsI correction to net spectrum
+ corrected_net = correct_csi_nonlinear(net_rate)
+ net_counts = corrected_net * live_time
bg_available = True
else:
net_counts = counts
bg_available = False
+ net_counts_list = [round(float(c), 1) for c in net_counts]
return {
"timestamp": state.get("timestamp", ""),
"cumulated_live_time_s": live_time,
"background_subtracted": bg_available,
- "channels": list(range(NUM_CHANNELS)),
+ "channels": clip_to_range(list(range(NUM_CHANNELS))),
"energy_kev": energy_axis(),
- "counts": [round(float(c), 1) for c in net_counts],
- "raw_counts": state.get("counts", [])[:NUM_CHANNELS],
+ "counts": clip_to_range(net_counts_list),
+ "raw_counts": clip_to_range(state.get("counts", [])[:NUM_CHANNELS]),
}
\ No newline at end of file
diff --git a/web/app/theoretical_bg.py b/web/app/theoretical_bg.py
index 30ba036..a5ab85e 100644
--- a/web/app/theoretical_bg.py
+++ b/web/app/theoretical_bg.py
@@ -30,7 +30,7 @@ def generate_continuum_only(cps: float = 6.0, live_time_s: float = 3600.0):
energy_axis = ENERGY_OFFSET + ENERGY_SLOPE * channels
total_counts = cps * live_time_s
- # Try calibrated spline first
+ # Try calibrated spline first (fits measured background)
continuum_cps = _get_continuum_cps()
if continuum_cps is not None and len(continuum_cps) == NUM_CHANNELS:
@@ -39,7 +39,7 @@ def generate_continuum_only(cps: float = 6.0, live_time_s: float = 3600.0):
if continuum.sum() > 0:
continuum *= total_counts / continuum.sum()
else:
- # Fallback: simple parametric model
+ # Fallback: parametric model
continuum = _fallback_continuum(energy_axis, total_counts)
return {
diff --git a/web/static/css/style.css b/web/static/css/style.css
index eb660ea..a52bbaf 100644
--- a/web/static/css/style.css
+++ b/web/static/css/style.css
@@ -68,18 +68,29 @@ nav a {
nav a:hover { background: rgba(255,255,255,0.1); }
nav a.active { color: var(--accent-bright); border-bottom: 2px solid var(--accent-bright); }
-main { padding: 16px; }
+main { padding: 12px 0; }
.tab-content { display: none; }
.tab-content.active { display: block; }
+.controls, .bg-stats, #isotopes-table, #peaks-table, .history-item, .chart-header {
+ margin-left: 12px;
+ margin-right: 12px;
+}
+
.chart-container {
background: var(--bg-card);
border-radius: 8px;
padding: 12px;
- margin-bottom: 12px;
+ margin: 0 12px 12px 12px;
height: 450px;
+ width: calc(100% - 24px);
position: relative;
+ overflow: hidden;
+}
+
+.chart-container canvas {
+ display: block;
}
.controls {
diff --git a/web/static/index.html b/web/static/index.html
index 688fc24..4591d99 100644
--- a/web/static/index.html
+++ b/web/static/index.html
@@ -4,7 +4,7 @@
Radiacode 103 — Dashboard
-
+
@@ -83,12 +83,12 @@
-
+
-
+
-
-
+
+