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

improvement to background search and plotting

parent a09201b2
......@@ -237,7 +237,7 @@ def check_peaks(peaks, nwindows):
################################################################################
################################################################################
def get_wfm_windows(data=None, filename=None, nwindows=6, bg_threshold=1.5,
def get_wfm_windows(data=None, filename=None, nwindows=6, bg_threshold=0.05,
win_threshold=0.3, plot=False, gsmooth=0, xrange=None,
output_format="tof"):
......@@ -326,7 +326,7 @@ def get_wfm_windows(data=None, filename=None, nwindows=6, bg_threshold=1.5,
# 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
bg_value = bg_raw + bg_threshold * (ymax - bg_raw)
hmin = np.amin(hist)
hmax = np.amax(hist)
......@@ -370,7 +370,10 @@ def get_wfm_windows(data=None, filename=None, nwindows=6, bg_threshold=1.5,
# time-of-flight.
mpd = int(0.75 * float(i_end - i_start) / nwindows)
loop = True
max_loops = 10
nloops = 0
while loop:
nloops += 1
print("The minimum peak distance (mpd) is:",mpd)
# Find valleys using `detect_peaks` function from Marcos Duarte
peaks = detect_peaks(y, mpd=mpd, valley=True)
......@@ -408,6 +411,8 @@ def get_wfm_windows(data=None, filename=None, nwindows=6, bg_threshold=1.5,
if mpd < 1:
print("The peak detection has failed with too few or too many peaks")
loop = False
if nloops > max_loops:
loop = False
# Now for each valley, iterate to one side starting from the valley center and
# find the window edge. We start from the first valley, which is the second
......@@ -422,22 +427,23 @@ def get_wfm_windows(data=None, filename=None, nwindows=6, bg_threshold=1.5,
# Define left and right window edges
ledges = [i_start]
redges = []
y_minus_background = y - bg_value
for p in range(1,len(good_peaks)-1):
# Towards the right ===================
rmean = np.average(y[good_peaks[p]:good_peaks[p+1]])
rmean = np.average(y_minus_background[good_peaks[p]:good_peaks[p+1]])
# Find left edge iterating towards the right
for i in range(good_peaks[p]+1,good_peaks[p+1]):
if (y[i] - y[good_peaks[p]]) >= (win_threshold*(rmean-y[good_peaks[p]])):
if (y_minus_background[i] - y_minus_background[good_peaks[p]]) >= (win_threshold*(rmean-y_minus_background[good_peaks[p]])):
ledges.append(i)
break
# Towards the left =======================
lmean = np.average(y[good_peaks[p-1]:good_peaks[p]])
lmean = np.average(y_minus_background[good_peaks[p-1]:good_peaks[p]])
# Find left edge iterating towards the right
for i in range(good_peaks[p]-1,good_peaks[p-1],-1):
if (y[i] - y[good_peaks[p]]) >= (win_threshold*(lmean-y[good_peaks[p]])):
if (y_minus_background[i] - y_minus_background[good_peaks[p]]) >= (win_threshold*(lmean-y_minus_background[good_peaks[p]])):
redges.append(i)
break
......@@ -460,9 +466,9 @@ def get_wfm_windows(data=None, filename=None, nwindows=6, bg_threshold=1.5,
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))
for p in range(len(good_peaks)-1):
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)
rmean = np.average(y_minus_background[good_peaks[p]:good_peaks[p+1]])
ax.plot([x[good_peaks[p]],x[good_peaks[p+1]]],[rmean+bg_value,rmean+bg_value],color="lime", lw=2)
ax.plot([x[good_peaks[p]],x[good_peaks[p+1]]],[win_threshold*rmean+bg_value,win_threshold*rmean+bg_value], color="sienna", lw=2)
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_value, bg_value], "--", color="pink", lw=1, label="Background threshold")
......@@ -486,7 +492,11 @@ def get_wfm_windows(data=None, filename=None, nwindows=6, bg_threshold=1.5,
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")
if isinstance(plot, str):
fname = plot
else:
fname = "wfm_auto_frames.pdf"
fig.savefig(fname,bbox_inches="tight")
# Format the output
output_format = output_format.lower() # convert to lowercase
......@@ -526,7 +536,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=1.5, dest="bg_threshold",
parser.add_argument("-b","--bg-threshold", type=float, default=0.05, 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