From db300bb1caaa5eb6f364907ddd900117cd260e28 Mon Sep 17 00:00:00 2001
From: Guenter Quast <guenter.quast@online.de>
Date: Sat, 13 Apr 2024 14:52:23 +0200
Subject: [PATCH] nicer output graphics, kafe2 output otpional
---
examples/findPeaks.py | 89 +++++++++++++++++++++++++++++--------------
1 file changed, 60 insertions(+), 29 deletions(-)
diff --git a/examples/findPeaks.py b/examples/findPeaks.py
index 7b0b676..80a88d1 100755
--- a/examples/findPeaks.py
+++ b/examples/findPeaks.py
@@ -56,6 +56,7 @@ parser.add_argument('-t','--threshold',
help = "threshold, i.e. minimum valid channel nmber (51)",
type=float, default=51)
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')
args = parser.parse_args()
@@ -69,6 +70,10 @@ fit_range_factor = args.fit_range_factor # fit range = fit_range_factor * fwhm
min_channel = args.threshold # threshold for min. valid channel number
verbose = args.verbose
plot = not args.noplot
+kafe2_plots = args.kafe2_plots
+
+# factor to convert Gaussian Sigma to full-width-half-maximum
+sig2fwhm = 2.3548
# read data
hst = np.loadtxt(filename, dtype=np.uint32)
@@ -97,6 +102,7 @@ hdata.label = "spectrum"
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 )")
@@ -132,45 +138,70 @@ for i, p in enumerate(peaks):
if verbose:
fit.report() # optionally, report fit results
-# Ergebnisse anzeigen und speichern:
+# show and save results
print("{} Peak {}: {:.2f} ± {:.2g} ( {:.2f} {:.1f}% {:.1f}% )".format(i, p,
fit.parameter_values[1],
fit.parameter_errors[1],
fit.parameter_values[2],
100*fit.parameter_values[2]/fit.parameter_values[1],
- 235.48*fit.parameter_values[2]/fit.parameter_values[1]
+ 100*sig2fwhm*fit.parameter_values[2]/fit.parameter_values[1]
))
fit_results = np.append(fit_results, fit)
- plot_ranges[i][0] = p-wid
- plot_ranges[i][1] = p+wid+1
+ plot_ranges[i][0] = max(0, p - wid)
+ plot_ranges[i][1] = min(len(hst), p + wid+1)
+ plot_ranges2[i][0] = max(0, p - 2*wid)
+ plot_ranges2[i][1] = min(len(hst), p + 2*wid+1)
# plot results
-sprctrum_fig = plt.figure('Spectrum', figsize=(10, 5))
-xhst = np.linspace(0, len(hst), len(hst), endpoint=False) +0.5
-plt.errorbar(xhst, hst, yerr=np.sqrt(hst),
- fmt='.', color='grey', alpha=0.1)
-for i, fit in enumerate(fit_results):
- xplt = np.linspace(plot_ranges[i][0],plot_ranges[i][1],
- 10*int((plot_ranges[i][1]-plot_ranges[i][0])))
- plt.plot(xplt, gauss_plus_bkg(xplt, *fit.parameter_values),
- linestyle = 'solid', linewidth=3)
- mu = fit.parameter_values[1]
- sig = fit.parameter_values[2]
- mx = gauss_plus_bkg(mu, *fit.parameter_values)
- h = fit.parameter_values[0]/np.sqrt(2*np.pi)/sig
- plt.vlines(mu, mx - h, mx,
- linewidth=3)
- plt.vlines(mu, 0, mx - h,
- linewidth=2, linestyle='dashed')
-
-plt.show()
-
+fig = plt.figure('Spectrum', figsize=(12, 10))
+fig.suptitle("Gamma spectrum " + filename)
+fig.subplots_adjust(left=0.12, bottom=0.1, right=0.95, top=0.95, wspace=None, hspace=0.1)
+ax0=fig.add_subplot(211)
+ax1=fig.add_subplot(212)
+ax0.set_ylabel('Entries per Channel')
+ax0.grid(linestyle='dotted', which='both')
+ax0.set_xticklabels([])
+ax1.set_ylabel('Entries per Channel')
+ax1.set_xlabel('Channel #')
+ax1.grid(linestyle='dotted', which='both')
-# kafe2 plots
if plot:
- fit_results = np.append(fit_results, hdata)
- plot_ranges[len(peaks+1)][1] = len(hst)
- kafe2.plot(fits = fit_results, x_range = list(map(tuple, plot_ranges)),
+ xhst = np.linspace(0, len(hst), len(hst), endpoint=False) +0.5
+ ax0.errorbar(xhst, hst, yerr=np.sqrt(hst), fmt='.', color='grey', markersize=2, linewidth=2)
+ ax1.errorbar(xhst, hst, yerr=np.sqrt(hst), fmt='.', color='grey', alpha=0.1)
+ colors = ('cadetblue', 'orange', 'olive', 'orchid', 'turquoise', 'tomato', 'green', 'pink', 'salmon', 'yellowgreen' )
+ for i, fit in enumerate(fit_results):
+ # plot fitted peak in fit range
+ colr = colors[i%10]
+ 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)
+ # 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),
+ linestyle = 'dotted', linewidth=2, color = colr)
+ # show fitted peak properties
+ mu = fit.parameter_values[1]
+ sig = fit.parameter_values[2]
+ fwhm = sig2fwhm * sig
+ mx = gauss_plus_bkg(mu, *fit.parameter_values)
+ h = fit.parameter_values[0]/np.sqrt(2*np.pi)/sig
+ ax1.vlines(mu, mx - h, mx,
+ linewidth=3, color = 'goldenrod')
+ 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')
+
+ if kafe2_plots:
+ # kafe2 plots
+ fit_results = np.append(fit_results, hdata)
+ plot_ranges[len(peaks+1)][1] = len(hst)
+ kafe2.plot(fits = fit_results, x_range = list(map(tuple, plot_ranges)),
x_label = "Channel number", y_label = "Counts",
font_scale = 0.5, save=False)
-
+ else:
+ plt.show()
+
--
GitLab