From 66bcd1533aeff72fe7ce403a9f4a45213ff63b1f Mon Sep 17 00:00:00 2001
From: Guenter Quast <guenter.quast@kit.edu>
Date: Sun, 14 Apr 2024 22:16:27 +0200
Subject: [PATCH] improved graphics
---
examples/findPeaks.py | 67 +++++++++++++++++++++----------------------
1 file changed, 33 insertions(+), 34 deletions(-)
diff --git a/examples/findPeaks.py b/examples/findPeaks.py
index 8333da8..0c81fbd 100755
--- a/examples/findPeaks.py
+++ b/examples/findPeaks.py
@@ -3,10 +3,9 @@
algorithm: perform a tow-stage determination of peak posisions:
- 1. scipi.signal.find_peaks() is run on a smoothed histogram
- of channel counts.
- 1. the output is used as a starting point for histogram fits
- with kafe2 based on the binned neg-log-likelhood cost function.
+ 1. scipi.signal.find_peaks() is run on a smoothened histogram of channel counts.
+ 2. the characteristica of identified peaks are used as starting points for
+ histogram fits with kafe2 based on the binned neg-log-likelhood method.
"""
import argparse
@@ -20,7 +19,6 @@ from scipy.signal import find_peaks
from PhyPraKit import meanFilter
from kafe2 import HistFit, HistContainer, Plot, plot
-
#define fit function
def gauss_plus_bkg(x, Ns=1000, mu=1000, sig=50., Nb=100., s=0., mn=0, mx=1.) :
"""Gaussian shape on linear background;
@@ -39,7 +37,7 @@ def gauss_plus_bkg(x, Ns=1000, mu=1000, sig=50., Nb=100., s=0., mn=0, mx=1.) :
B = (1 + (x-mu)*s) / (s/2*(mx**2 - mn**2) + (1-s*mu)*(mx - mn))
return Ns*S + Nb*B
-# parse input
+# define command line arguments ...
parser = argparse.ArgumentParser(description=__doc__)
parser.add_argument('filename', nargs='?', default='Co60.hst')
parser.add_argument('-p','--prominence',
@@ -60,7 +58,7 @@ parser.add_argument('-t','--threshold',
parser.add_argument('-n','--noplot', action='store_true')
parser.add_argument('-k','--kafe2_plots', action='store_true')
parser.add_argument('-v','--verbose', action='store_true')
-
+# ... and parse command line input
args = parser.parse_args()
filename=args.filename
# some constants for peak-finder
@@ -82,54 +80,54 @@ hst = np.loadtxt(filename, dtype=np.uint32)
bin_edges= np.linspace(0, len(hst), len(hst)+1, endpoint=True)
# find maxima with scipy.signal.find_peaks()
-# smoothen data to reduce statistical noise before running peak finder
+# firts, smoothen data to reduce statistical noise
hst_s = meanFilter(hst, max(1, min_width//4))
+# search for peaks
peaks, peak_props = find_peaks(hst_s,
prominence=min_prominence,
width=min_width, rel_height=rel_height,
wlen=350)
-print(len(peaks), "peaks found: ", peaks)
-print("prominences: ", peak_props["prominences"])
-print("left_bases: ", peak_props["left_bases"])
-print("right_bases: ", peak_props["right_bases"])
-print("fwhms: ", peak_props["widths"])
+if verbose:
+ print(len(peaks), "peaks found: ", peaks)
+ print("prominences: ", peak_props["prominences"])
+ print("left_bases: ", peak_props["left_bases"])
+ print("right_bases: ", peak_props["right_bases"])
+ print("fwhms: ", peak_props["widths"])
+
prominences = peak_props["prominences"]
left_bases = peak_props["left_bases"]
right_bases = peak_props["right_bases"]
fwhms = peak_props["widths"]
# perform fit for precise determination
-# put data in kafe2 HistContainer
+# put data in kafe2 HistContainer
hdata = HistContainer(bin_edges=bin_edges)
hdata.set_bins(hst)
-hdata.label = "spectrum"
-
+hdata.label = "spectrum data"
+# initialize arrays for output
fit_results=[]
plot_ranges = np.zeros([len(peaks)+1, 2])
plot_ranges2 = np.zeros([len(peaks)+1, 2])
-
print("\n"+"*==* fit results:")
-print(" mu ± d_mu ( sig sig/mu FWHM/mu )")
+print(" mu ± d_mu ( sig sig/mu FWHM/mu )")
+# run fits for all identified peaks
for i, p in enumerate(peaks):
wid = int(fit_range_factor*fwhms[i])
base = (hst[left_bases[i]]+hst[right_bases[i]])/2.
- # print(p, wid, base)
-
mn = max(min_channel, int(p-wid))
mx = int(p+wid)
bins = range(mn, mx+1)
counts = hst[mn : mx]
- label = "peak "+str(i)
hfit_data = HistContainer(bin_edges = bins)
hfit_data.label="peak " + str(i)
hfit_data.set_bins(counts)
#Fit
- fit = HistFit(data = hfit_data, model_function = gauss_plus_bkg, density=False)
- # !!! Startwerte setzen !!!
+ fit = HistFit(data = hfit_data, model_function = gauss_plus_bkg, density=False)
+ # !!! set initial values
fit.set_parameter_values(mu=p, sig=wid, Nb=base*(mx-mn))
fit.fix_parameter("mn", mn)
fit.fix_parameter("mx", mx)
- # einige Parameter auf "vernünftige" Bereiche festlegen
+ # limit parameters to reasonable values
fit.limit_parameter("sig", 0., None)
fit.limit_parameter("Ns", 0., None)
#fit.limit_parameter("Nb", 0., None) # biases uncertainty if at zero boundary
@@ -143,7 +141,7 @@ for i, p in enumerate(peaks):
fit.report() # optionally, report fit results
# show and save results
- print("{} Peak {}: {:.2f} ± {:.2g} ( {:.2f} {:.1f}% {:.1f}% )".format(i, p,
+ print("{} peak@{}: {:.2f} ± {:.2g} ( {:.2f} {:.1f}% {:.1f}% )".format(i+1, p,
fit.parameter_values[1],
fit.parameter_errors[1],
fit.parameter_values[2],
@@ -172,17 +170,18 @@ ax1.grid(linestyle='dotted', which='both')
if plot:
# show spectrum and result of find_peaks
xhst = np.linspace(0, len(hst), len(hst), endpoint=False) +0.5
- ax0.plot(xhst, hst_s, 'b-', linewidth=3, label = 'smoothed channel counts')
- ax0.errorbar(xhst, hst, yerr=np.sqrt(hst), label='channel counts',
- fmt='.', color='grey', markersize=2, linewidth=2, alpha=0.5)
- ax0.plot(peaks, hst_s[peaks], '*', color='red', markersize=10, label='result of find_peaks()')
+ ax0.plot(xhst, hst_s, 'b-', linewidth=1, zorder=2, label = 'smoothed channel counts')
+ ax0.errorbar(xhst, hst, yerr=np.sqrt(hst), zorder=1, label='channel counts',
+ fmt='.', color='grey', markersize=2, linewidth=2, alpha=0.5, )
+ ax0.plot(peaks, hst_s[peaks], '*', color='red', markersize=10, zorder=3,
+ label='result of find_peaks()')
ax0.legend(loc = 'best')
# show fitted peaks
- ax1.errorbar(xhst, hst, yerr=np.sqrt(hst), fmt='.', color='grey', alpha=0.1,
+ ax1.errorbar(xhst, hst, yerr=np.sqrt(hst), fmt='.', color='grey', alpha=0.25, zorder=1,
label = 'channel counts')
# select colors for peaks
- pcolors = ('cadetblue', 'orange', 'olive', 'orchid', 'turquoise',
+ pcolors = ('steelblue', 'darkorange', 'green', 'orchid', 'turquoise',
'tomato', 'green', 'pink', 'salmon', 'yellowgreen' )
for i, fit in enumerate(fit_results):
# plot fitted peak in fit range
@@ -190,12 +189,12 @@ if plot:
xplt = np.linspace(plot_ranges[i][0], plot_ranges[i][1],
10*int((plot_ranges[i][1]-plot_ranges[i][0])))
ax1.plot(xplt, gauss_plus_bkg(xplt, *fit.parameter_values),
- linestyle = 'solid', linewidth=3, color= colr,
+ linestyle = 'solid', linewidth=3, color= colr, zorder=2,
label='peak ' + str(i+1) + '@' + str(int(10*fit.parameter_values[1])/10.) )
# plot fitted peak near fit region
xplt2 = np.linspace(plot_ranges2[i][0], plot_ranges2[i][1],
10*int((plot_ranges[i][1]-plot_ranges[i][0])))
- ax1.plot(xplt2, gauss_plus_bkg(xplt2, *fit.parameter_values),
+ ax1.plot(xplt2, gauss_plus_bkg(xplt2, *fit.parameter_values), zorder=2,
linestyle = 'dotted', linewidth=2, color = colr)
# show fitted peak properties
mu = fit.parameter_values[1]
@@ -208,7 +207,7 @@ if plot:
ax1.vlines(mu, 0, mx - h,
linewidth=1, linestyle='dashed', color='goldenrod')
ax1.hlines(mx-h/2, mu-fwhm/2, mu+fwhm/2,
- linewidth=2, color='gold')
+ linewidth=2, color='goldenrod')
ax1.legend(loc = 'best')
if kafe2_plots:
--
GitLab