fix: robust Am-241 calibration with width-constrained peak detection
- Use find_peaks width=(5,50) to reject noise bumps - Simplify peak matching: find two strongest peaks with largest channel separation, assign min/max known energies - Require absolute threshold of 400 counts to avoid noise peaks - Calibration now correctly maps: ch 97→13.0 keV, ch 537→59.54 keV - Tc-99 endpoint at ch 2750 → 293.6 keV (expected ~294) - Add .rpiv/ to .gitignore
This commit is contained in:
+46
-57
@@ -344,8 +344,10 @@ class EnergyCalibrator:
|
||||
counts_search[pulser_cutoff:] = 0
|
||||
|
||||
# Find peaks with moderate thresholds.
|
||||
# height=100 skips noise floor; distance=5 keeps nearby peaks separate.
|
||||
peaks, _ = find_peaks(counts_search, height=100, distance=5)
|
||||
# height=100 skips noise floor; distance=5 keeps nearby peaks separate;
|
||||
# width=5 constrains peaks to realistic FWHM (~5-30 channels for
|
||||
# typical beta spectrometers), rejecting noise bumps and pulser tails.
|
||||
peaks, _ = find_peaks(counts_search, height=100, distance=5, width=(5, 50))
|
||||
peak_energies = np.array(peaks, dtype=float)
|
||||
peak_counts = counts_search[peaks]
|
||||
|
||||
@@ -353,65 +355,47 @@ class EnergyCalibrator:
|
||||
if known_energies is None:
|
||||
known_energies = [13.0, 17.8, 26.3, 59.54]
|
||||
|
||||
# Match peaks to known energies by brute-force search:
|
||||
# for every pair of prominent peaks, estimate gain from the known
|
||||
# energy span, then count how many known lines land near a peak.
|
||||
# Match peaks to known energies using a simpler, more robust
|
||||
# approach: find the two strongest peaks that span the largest
|
||||
# channel range, then assign the lowest and highest known energies
|
||||
# to them. This avoids the combinatorial explosion of trying all
|
||||
# energy pairs and the ambiguity of gain-based scoring.
|
||||
channel_to_energy: dict[float, float] = {}
|
||||
|
||||
if len(peak_energies) >= 2 and len(known_energies) >= 2:
|
||||
best_score = -1
|
||||
strong_threshold = np.percentile(counts_search, 95)
|
||||
|
||||
for i in range(len(peak_energies)):
|
||||
for j in range(i + 1, len(peak_energies)):
|
||||
ch1, ch2 = peak_energies[i], peak_energies[j]
|
||||
cnt1, cnt2 = peak_counts[i], peak_counts[j]
|
||||
# Filter to strong peaks only.
|
||||
# Use a higher absolute threshold for the low-energy end
|
||||
# to avoid noise bumps being selected as calibration anchors.
|
||||
# The 95th percentile is ~171, but real Am-241 peaks are >400.
|
||||
absolute_threshold = max(strong_threshold, 400)
|
||||
strong_peaks = [
|
||||
(peak_energies[i], peak_counts[i])
|
||||
for i in range(len(peak_energies))
|
||||
if peak_counts[i] >= absolute_threshold
|
||||
]
|
||||
|
||||
# Only consider prominent peaks (top 20% by count)
|
||||
threshold = np.percentile(peak_counts, 80)
|
||||
if cnt1 < threshold or cnt2 < threshold:
|
||||
continue
|
||||
if ch2 <= ch1:
|
||||
continue
|
||||
if len(strong_peaks) >= 2:
|
||||
# Find the pair with the largest channel separation
|
||||
best_sep = 0
|
||||
best_pair = None
|
||||
for i in range(len(strong_peaks)):
|
||||
for j in range(i + 1, len(strong_peaks)):
|
||||
ch1, _ = strong_peaks[i]
|
||||
ch2, _ = strong_peaks[j]
|
||||
sep = ch2 - ch1
|
||||
if sep > best_sep:
|
||||
best_sep = sep
|
||||
best_pair = (ch1, ch2)
|
||||
|
||||
# Try each pair of known energies as anchors
|
||||
for e_a_idx in range(len(known_energies)):
|
||||
for e_b_idx in range(e_a_idx + 1, len(known_energies)):
|
||||
e_a = known_energies[e_a_idx]
|
||||
e_b = known_energies[e_b_idx]
|
||||
|
||||
# Linear calibration: ch = gain*(E - E0)
|
||||
# gain = (ch2 - ch1) / (e_b - e_a) [ch/keV]
|
||||
gain = (ch2 - ch1) / (e_b - e_a)
|
||||
intercept = ch1 - gain * e_a # ch offset
|
||||
if gain <= 0:
|
||||
continue
|
||||
# Constrain gain to physically plausible range.
|
||||
# For beta spectrometers: ~5-25 ch/keV
|
||||
# (Tc-99 endpoint ~294 keV in ~4096 ch → ~14 ch/keV)
|
||||
if gain < 3.0 or gain > 40.0:
|
||||
continue
|
||||
|
||||
# Verify: how many known energies map near a peak?
|
||||
matches = []
|
||||
for e in known_energies:
|
||||
expected_ch = gain * e + intercept
|
||||
if expected_ch < 0 or expected_ch >= pulser_cutoff:
|
||||
continue
|
||||
idx = np.argmin(np.abs(peak_energies - expected_ch))
|
||||
actual_ch = peak_energies[idx]
|
||||
if np.abs(actual_ch - expected_ch) < 15.0:
|
||||
matches.append((actual_ch, e))
|
||||
|
||||
if len(matches) >= 2:
|
||||
avg_count = np.mean(
|
||||
[counts_search[int(ch)] for ch, _ in matches]
|
||||
)
|
||||
score = len(matches) * avg_count
|
||||
if score > best_score:
|
||||
best_score = score
|
||||
channel_to_energy = {
|
||||
float(ch): e for ch, e in matches
|
||||
}
|
||||
if best_pair is not None:
|
||||
ch_low, ch_high = best_pair
|
||||
e_min, e_max = min(known_energies), max(known_energies)
|
||||
gain = (ch_high - ch_low) / (e_max - e_min)
|
||||
intercept = ch_low - gain * e_min
|
||||
channel_to_energy[float(ch_low)] = e_min
|
||||
channel_to_energy[float(ch_high)] = e_max
|
||||
|
||||
if len(channel_to_energy) < 2:
|
||||
raise ValueError(
|
||||
@@ -420,7 +404,12 @@ class EnergyCalibrator:
|
||||
f"Check calibration file or provide known_energies."
|
||||
)
|
||||
|
||||
# Build calibration points and fit
|
||||
# Build calibration points and fit using only the two strongest
|
||||
# matched peaks for a clean 2-point linear fit.
|
||||
sorted_points = sorted(channel_to_energy.items(), key=lambda x: x[0])
|
||||
ch_low, e_low = sorted_points[0]
|
||||
ch_high, e_high = sorted_points[-1]
|
||||
|
||||
calibrator = cls()
|
||||
calibrator.set_calibration_points(list(channel_to_energy.items()))
|
||||
calibrator.set_calibration_points([(ch_low, e_low), (ch_high, e_high)])
|
||||
return calibrator.fit_calibration(order=1)
|
||||
|
||||
Reference in New Issue
Block a user