diff --git a/becquerel/core/peakfinder.py b/becquerel/core/peakfinder.py index a4f56286..a7f95dd4 100644 --- a/becquerel/core/peakfinder.py +++ b/becquerel/core/peakfinder.py @@ -302,11 +302,12 @@ def add_peak(self, xpeak): def plot(self, facecolor="red", linecolor="red", alpha=0.5, peaks=True): """Plot the peak signal-to-noise ratios calculated using the kernel.""" bin_edges = self.spectrum.bin_edges_raw + bin_centers = self.spectrum.bin_centers_raw if facecolor is not None: - plt.fill_between(bin_edges[:-1], self.snr, 0, color=facecolor, alpha=alpha) + plt.fill_between(bin_centers, self.snr, 0, color=facecolor, alpha=alpha) if linecolor is not None: - plt.plot(bin_edges[:-1], self.snr, "-", color=linecolor) + plt.plot(bin_centers, self.snr, "-", color=linecolor) if peaks: for cent, snr, fwhm in zip(self.centroids, self.snrs, self.fwhms): plt.plot([cent] * 2, [0, snr], "b-", lw=1.5) @@ -314,7 +315,7 @@ def plot(self, facecolor="red", linecolor="red", alpha=0.5, peaks=True): plt.plot( [cent - fwhm / 2, cent + fwhm / 2], [snr / 2] * 2, "b-", lw=1.5 ) - plt.xlim(0, bin_edges.max()) + plt.xlim(0, bin_centers.max()) plt.ylim(0) plt.xlabel("x") plt.ylabel("SNR") @@ -389,15 +390,9 @@ def find_peaks(self, xmin=None, xmax=None, min_snr=2, max_num=40): max_num = int(max_num) if max_num < 1: raise PeakFinderError("Must keep at least 1 peak, not {}".format(max_num)) - # calculate the first derivative and second derivatives of the SNR - d1 = (self.snr[2:] - self.snr[:-2]) / 2 - d1 = np.append(0, d1) - d1 = np.append(d1, 0) - d2 = self.snr[2:] - 2 * self.snr[1:-1] + self.snr[:-2] - d2 = np.append(0, d2) - d2 = np.append(d2, 0) + # find maxima - peak = (d1[2:] < 0) & (d1[:-2] > 0) & (d2[1:-1] < 0) + peak = (self.snr[:-2] < self.snr[1:-1]) & (self.snr[1:-1] >= self.snr[2:]) peak = np.append(False, peak) peak = np.append(peak, False) # select peaks using SNR and centroid criteria @@ -415,3 +410,4 @@ def find_peaks(self, xmin=None, xmax=None, min_snr=2, max_num=40): self.backgrounds = self.backgrounds[-max_num:] # sort by centroid self.sort_by(self.centroids) + self.peak = peak diff --git a/tests/autocal_test.py b/tests/autocal_test.py index 407840e2..9d94a7f3 100644 --- a/tests/autocal_test.py +++ b/tests/autocal_test.py @@ -87,13 +87,45 @@ def test_peakfinder(): kernel = bq.GaussianPeakFilter(500, 50, fwhm_at_0=10) finder = bq.PeakFinder(spec1, kernel) finder.find_peak(500, min_snr=3.0) - assert np.isclose(finder.channels[0], 485.5) + assert np.isclose(finder.centroids[0], 485.5) finder.reset() finder.find_peaks() - assert len(finder.channels) == 9 + assert len(finder.centroids) == 9 finder.reset() finder.find_peaks(min_snr=0.5, xmin=50, xmax=1000, max_num=10) - assert len(finder.channels) == 10 + assert len(finder.centroids) == 10 + + +def test_peakfinder_double_peaks(): + """Test that the new peakfinder logic correctly gives only one peak, even + when min_sep is not specified, while old logic gives two peaks. + """ + + # Toy data: gaussian lineshape + x = np.arange(31) + y = 1.5 * np.exp(-0.5 * (x[:-1] - 15) ** 2 / (2 * 3 ** 2)) + spec = bq.Spectrum(bin_edges_raw=x, counts=y) + kernel = bq.GaussianPeakFilter(15, 7.5, fwhm_at_0=7.5) + finder = bq.PeakFinder(spec, kernel) + + # OLD code from PeakFinder.find_peaks + d1 = (finder.snr[2:] - finder.snr[:-2]) / 2 + d1 = np.append(0, d1) + d1 = np.append(d1, 0) + d2 = finder.snr[2:] - 2 * finder.snr[1:-1] + finder.snr[:-2] + d2 = np.append(0, d2) + d2 = np.append(d2, 0) + peak_old = (d1[2:] < 0) & (d1[:-2] > 0) & (d2[1:-1] < 0) + peak_old = np.append(False, peak_old) + peak_old = np.append(peak_old, False) + assert peak_old.sum() == 2 + + # NEW peakfind + finder.find_peaks(min_snr=1.0) + assert len(finder.centroids) == 1 + assert finder.centroids[0] in finder.spectrum.bin_centers_raw[peak_old] + peak_new = np.isclose(finder.spectrum.bin_centers_raw, finder.centroids[0]) + assert peak_new.sum() == 1 def test_peakfinder_exceptions(): @@ -192,7 +224,7 @@ def test_autocal_spec1(): finder = bq.PeakFinder(spec1, kernel) cal = bq.AutoCalibrator(finder) finder.find_peaks(min_snr=1, xmin=50) - assert len(cal.peakfinder.channels) == 10 + assert len(cal.peakfinder.centroids) == 10 cal.fit( REQUIRED, optional=OPTIONAL, @@ -210,7 +242,7 @@ def test_autocal_one_line(): finder = bq.PeakFinder(spec1, kernel) cal = bq.AutoCalibrator(finder) finder.find_peaks(min_snr=8, xmin=50) - assert len(cal.peakfinder.channels) == 1 + assert len(cal.peakfinder.centroids) == 1 cal.fit( [1460.82], gain_range=[2.5, 4.0], @@ -230,7 +262,7 @@ def test_autocal_exceptions(): cal = bq.AutoCalibrator(finder) # only one peak found, but multiple energies required finder.find_peaks(min_snr=10, xmin=50) - assert len(cal.peakfinder.channels) == 1 + assert len(cal.peakfinder.centroids) == 1 with pytest.raises(bq.AutoCalibratorError): cal.fit( REQUIRED, @@ -240,7 +272,7 @@ def test_autocal_exceptions(): # multiple peaks found, but only one energy given finder.reset() finder.find_peaks(min_snr=1, xmin=50) - assert len(cal.peakfinder.channels) == 10 + assert len(cal.peakfinder.centroids) == 10 with pytest.raises(bq.AutoCalibratorError): cal.fit( [1460.82], @@ -250,7 +282,7 @@ def test_autocal_exceptions(): # more energies required than peaks found finder.reset() finder.find_peaks(min_snr=4, xmin=50) - assert len(cal.peakfinder.channels) == 4 + assert len(cal.peakfinder.centroids) == 4 with pytest.raises(bq.AutoCalibratorError): cal.fit( [238.63, 351.93, 609.32, 1460.82, 2614.3], @@ -265,7 +297,7 @@ def test_autocal_no_fit(): finder = bq.PeakFinder(spec1, kernel) cal = bq.AutoCalibrator(finder) finder.find_peaks(min_snr=2, xmin=50) - assert len(cal.peakfinder.channels) == 8 + assert len(cal.peakfinder.centroids) == 8 with pytest.raises(bq.AutoCalibratorError): cal.fit( REQUIRED, @@ -282,7 +314,7 @@ def test_autocal_plot(): finder = bq.PeakFinder(spec1, kernel) cal = bq.AutoCalibrator(finder) finder.find_peaks(min_snr=1, xmin=50) - assert len(cal.peakfinder.channels) == 10 + assert len(cal.peakfinder.centroids) == 10 cal.fit( REQUIRED, optional=OPTIONAL, @@ -299,15 +331,15 @@ def test_autocal_spec2(): kernel = bq.GaussianPeakFilter(3700, 10, 5) finder = bq.PeakFinder(spec2, kernel) cal = bq.AutoCalibrator(finder) - cal.peakfinder.find_peaks(min_snr=15, xmin=1000) - assert len(cal.peakfinder.channels) == 12 + cal.peakfinder.find_peaks(min_snr=20, xmin=1000) + assert len(cal.peakfinder.centroids) == 9 cal.fit( REQUIRED, optional=OPTIONAL, gain_range=[0.1, 0.6], de_max=5.0, ) - assert len(cal.fit_channels) == 6 + assert len(cal.fit_channels) == 3 assert np.isclose(cal.gain, 0.3785, rtol=1e-2) @@ -317,7 +349,7 @@ def test_autocal_spec3(): finder = bq.PeakFinder(spec3, kernel) cal = bq.AutoCalibrator(finder) cal.peakfinder.find_peaks(min_snr=3, xmin=100) - assert len(cal.peakfinder.channels) == 8 + assert len(cal.peakfinder.centroids) == 8 # this fit succeeds but misidentifies the lines cal.fit( [609.32, 1460.82], @@ -344,7 +376,7 @@ def test_autocal_spec4(): finder = bq.PeakFinder(spec4, kernel) cal = bq.AutoCalibrator(finder) cal.peakfinder.find_peaks(min_snr=3, xmin=100) - assert len(cal.peakfinder.channels) == 7 + assert len(cal.peakfinder.centroids) == 7 cal.fit( [356.0129, 661.657, 1460.82], optional=[911.20, 1120.294, 1620.50, 1764.49, 2118.514, 2614.3],