diff --git a/exp_data/calibration.py b/exp_data/calibration.py index 0fc853a..a726f97 100644 --- a/exp_data/calibration.py +++ b/exp_data/calibration.py @@ -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) diff --git a/notebooks/output.root b/notebooks/output.root index 209df08..d91e06e 100644 Binary files a/notebooks/output.root and b/notebooks/output.root differ