Commit a09201b2 authored by Neil Vaytet's avatar Neil Vaytet
Browse files

Fixed automatic peak dismissal and updated background search method

parent c8a50d1a
......@@ -237,7 +237,7 @@ def check_peaks(peaks, nwindows):
################################################################################
################################################################################
def get_wfm_windows(data=None, filename=None, nwindows=6, bg_threshold=1.0e-4,
def get_wfm_windows(data=None, filename=None, nwindows=6, bg_threshold=1.5,
win_threshold=0.3, plot=False, gsmooth=0, xrange=None,
output_format="tof"):
......@@ -323,22 +323,29 @@ def get_wfm_windows(data=None, filename=None, nwindows=6, bg_threshold=1.0e-4,
ymin = np.amin(y[lbound:rbound+1])
ymax = np.amax(y[lbound:rbound+1])
# Determine background threshold
i_bg = None
y_normed = y / np.sum(y)
sorted = np.argsort(y_normed)
grad = np.gradient(y_normed[sorted])
cumsum = np.cumsum(grad)
for i in range(len(cumsum)):
if cumsum[i] > bg_threshold:
i_bg = i
break
if i_bg is not None:
bg_value = y[sorted][i_bg]
print("Background value is: {}".format(bg_value))
else:
bg_value = 0.05*(ymax - ymin)
print("Warning: background value was not found in the conventional way")
# Determine background threshold by histogramming data
hist, edges = np.histogram(y[lbound:rbound+1], bins=50)
bg_raw = edges[np.argmax(hist) + 1]
bg_value = bg_raw * bg_threshold
hmin = np.amin(hist)
hmax = np.amax(hist)
# # Determine background threshold by looking for strong cumulative gradients
# i_bg = None
# y_normed = y / np.sum(y)
# sorted = np.argsort(y_normed)
# grad = np.gradient(y_normed[sorted])
# cumsum = np.cumsum(grad)
# for i in range(len(cumsum)):
# if cumsum[i] > bg_threshold:
# i_bg = i
# break
# if i_bg is not None:
# bg_value = y[sorted][i_bg]
# print("Background value is: {}".format(bg_value))
# else:
# bg_value = 0.05*(ymax - ymin)
# print("Warning: background value was not found in the conventional way")
# Find the leading and trailing edges; i.e. the leftmost and rightmost
# points that exceed the background value
......@@ -375,19 +382,19 @@ def get_wfm_windows(data=None, filename=None, nwindows=6, bg_threshold=1.0e-4,
good_peaks.append(i_end)
loop = check_peaks(good_peaks, nwindows)
# Try to remove peaks of there are too many
# Try to remove peaks if there are too many
if (len(good_peaks)-2) > (nwindows - 1):
print("Too many peaks found. Attempting automatic dismissal")
# if the value at the bottom of the valley is over the window
# threshold, then this is not a real valley
# if the value at the bottom of the valley is over the lowest window
# average, then this is not a real valley
temp_peaks = good_peaks[:]
for p in range(1,len(good_peaks)-1):
rmean = np.average(y[good_peaks[p]:good_peaks[p+1]])
lmean = np.average(y[good_peaks[p-1]:good_peaks[p]])
mean = 0.5*(lmean + rmean)
if y[good_peaks[p]] > (win_threshold*mean):
mean = min(lmean,rmean)
if y[good_peaks[p]] > mean:
temp_peaks[p] = -100
print("Removed peak number {} at x,y = {},{}".format(p,x[temp_peaks[p]],y[temp_peaks[p]]))
print("Removed peak number {} at x,y = {},{}".format(p,x[good_peaks[p]],y[good_peaks[p]]))
good_peaks = []
for p in range(len(temp_peaks)):
if temp_peaks[p] != -100:
......@@ -445,7 +452,10 @@ def get_wfm_windows(data=None, filename=None, nwindows=6, bg_threshold=1.0e-4,
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle
fig = plt.figure()
ax = fig.add_subplot(111)
ratio = 0.35
sizex = 20.0
fig.set_size_inches(sizex, ratio*sizex)
ax = fig.add_subplot(121)
colors = ["r","g","b","magenta","cyan","orange"]
for i in range(len(ledges)):
ax.add_patch(Rectangle((x[ledges[i]], ymin), (x[redges[i]]-x[ledges[i]]), (ymax-ymin), facecolor=colors[i%6], alpha=0.5))
......@@ -453,10 +463,9 @@ def get_wfm_windows(data=None, filename=None, nwindows=6, bg_threshold=1.0e-4,
rmean = np.average(y[good_peaks[p]:good_peaks[p+1]])
ax.plot([x[good_peaks[p]],x[good_peaks[p+1]]],[rmean,rmean],color="lime", lw=2)
ax.plot([x[good_peaks[p]],x[good_peaks[p+1]]],[win_threshold*rmean,win_threshold*rmean], color="sienna", lw=2)
bg = bg_threshold*(ymax - ymin)
ax.plot(x, data[:,1], color="k", lw=2, label="Raw data")
ax.plot(x, y, color="lightgrey", lw=1, label="Smoothed data")
ax.plot([np.amin(x),np.amax(x)], [bg, bg], "--", color="pink", lw=1, label="Background threshold")
ax.plot([np.amin(x),np.amax(x)], [bg_value, bg_value], "--", color="pink", lw=1, label="Background threshold")
ax.plot([x[good_peaks[0]],x[good_peaks[1]]], [-ymax,-ymax], color="lime", label="Window mean", lw=2)
ax.plot([x[good_peaks[0]],x[good_peaks[1]]], [-ymax,-ymax], color="sienna", label="Window mean", lw=2)
for p in range(1,len(good_peaks)-1):
......@@ -467,7 +476,16 @@ def get_wfm_windows(data=None, filename=None, nwindows=6, bg_threshold=1.0e-4,
ax.set_ylabel("Amplitude")
ax.set_ylim([0.0,1.05*np.amax(data[:,1])])
ax.plot(x[good_peaks[0]], -ymax, "o", color="r",label="Valleys")
ax.legend(loc=(0,1.02),ncol=4, fontsize=7)
ax.legend(loc=(0,1.02),ncol=4, fontsize=10)
ax2 = fig.add_subplot(122)
ax2.hist(y[lbound:rbound+1], bins=edges, orientation='horizontal', label="Amplitude frequency", facecolor='None')
ax2.set_xlim([hmin,1.05*hmax])
ax2.set_ylim([0.0,1.05*np.amax(data[:,1])])
ax2.plot([hmin, hmax], [bg_raw, bg_raw], "--", color="b", lw=3, label="Background estimation")
ax2.plot([hmin, hmax], [bg_value, bg_value], "--", color="pink", lw=1, label="Background threshold")
ax2.set_xlabel("Counts in background search")
ax2.set_ylabel("Amplitude")
ax2.legend(loc=(0,1.02),ncol=4, fontsize=10)
fig.savefig("wfm_auto_frames.pdf",bbox_inches="tight")
# Format the output
......@@ -508,7 +526,7 @@ if __name__ == "__main__":
parser.add_argument("-n", "--nwindows", type=int, default=6, dest="nwindows",
help="the number of windows to be found")
parser.add_argument("-b","--bg-threshold", type=float, default=0.05, dest="bg_threshold",
parser.add_argument("-b","--bg-threshold", type=float, default=1.5, dest="bg_threshold",
help="threshold above which we are no longer in"
"background signal, as a percentage of (ymax - ymin).")
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment