Commit 71b22409 authored by Neil Vaytet's avatar Neil Vaytet
Browse files

update to background finding

parent c959b978
......@@ -237,9 +237,9 @@ def check_peaks(peaks, nwindows):
################################################################################
################################################################################
def get_wfm_windows(data=None, filename=None, nwindows=6, bg_threshold=0.05,
def get_wfm_windows(data=None, filename=None, nwindows=6, bg_threshold=1.0e-4,
win_threshold=0.3, plot=False, gsmooth=0, xrange=None,
rebin_step_for_string_output=None):
output_format="tof"):
"""
Get the left and right edges of wave-frame multiplication windows.
......@@ -269,17 +269,19 @@ def get_wfm_windows(data=None, filename=None, nwindows=6, bg_threshold=0.05,
if you want to have no limit on that end.
For example: xrange=(1000,55000) or xrange=[5000,None]
- rebin_step_for_string_output: a number. If specified, the function will
return an array of strings that are to be
used as the input for the Rebin algorithm.
The number given as an argument will appear
as the bin step in the middle of each range.
For example, if the number supplied is 64,
the output will look like:
['15000,64,20000', '25000,64,30000', ...]
If unspecified, the function will return
the indices (NOT the Tof values) of the left
and right edges as arrays of values.
- output_format: a string. Can be 3 different things:
- "TOF" or "tof": returns an array of floats containing the
time-of-flight boundaries for the windows in the format
[frame1_left_edge, frame1_right_edge, frame2_left_edge, frame2_right_edge, ...]
- "indices": returns the indices of the left and right
edges as arrays of values.
- "rebin:64": If in this format, specified, the function
will return an array of strings that are to be used as
the input for the Rebin algorithm. The number given
after the colon will appear as the bin step in the
middle of each range. For example, if the number
supplied is 64, the output will look like:
['15000,64,20000', '25000,64,30000', ...]
"""
if data is None:
......@@ -292,7 +294,7 @@ def get_wfm_windows(data=None, filename=None, nwindows=6, bg_threshold=0.05,
x = data[:,0]
y = data[:,1]
print(x)
# Smooth the data with a gaussian filter
if gsmooth == 0:
gsmooth = max(int(nx/500),1)
......@@ -320,16 +322,34 @@ def get_wfm_windows(data=None, filename=None, nwindows=6, bg_threshold=0.05,
# Find min and max values
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")
# Find the leading and trailing edges; i.e. the leftmost and rightmost
# points that exceed the value `bg_threshold * (ymax - ymin)
# points that exceed the background value
i_start = 0
i_end = nx-1
for i in range(lbound,nx):
if y[i] > bg_threshold*(ymax - ymin):
if y[i] > bg_value:
i_start = i
break
for i in range(rbound,1,-1):
if y[i] > bg_threshold*(ymax - ymin):
if y[i] > bg_value:
i_end = i
break
......@@ -350,7 +370,7 @@ def get_wfm_windows(data=None, filename=None, nwindows=6, bg_threshold=0.05,
# Now filter out peaks that are between start and end
good_peaks = [i_start]
for p in peaks:
if (p > i_start+0.25*mpd) and (p < i_end-0.25*mpd):
if (p > i_start+0.5*mpd) and (p < i_end-0.5*mpd):
good_peaks.append(p)
good_peaks.append(i_end)
loop = check_peaks(good_peaks, nwindows)
......@@ -425,7 +445,7 @@ def get_wfm_windows(data=None, filename=None, nwindows=6, bg_threshold=0.05,
ax = fig.add_subplot(111)
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%5], alpha=0.5))
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)
......@@ -440,22 +460,29 @@ def get_wfm_windows(data=None, filename=None, nwindows=6, bg_threshold=0.05,
ax.plot(x[good_peaks[p]], y[good_peaks[p]], "o", color="r")
ax.plot(x[good_peaks[0]], y[good_peaks[0]], "o", color="deepskyblue", label="Leading edge")
ax.plot(x[good_peaks[-1]], y[good_peaks[-1]], "o", color="yellow", label="Trailing edge")
ax.set_xlabel("Time-of-flight")
ax.set_xlabel("Arrival time at detector")
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)
fig.savefig("wfm_auto_frames.pdf",bbox_inches="tight")
# If we want output in Rebin format, then construct the array of strings
if rebin_step_for_string_output is not None:
output = []
# Format the output
output_format = output_format.lower() # convert to lowercase
output = []
if output_format.startswith("rebin"):
# If we want output in Rebin format, then construct the array of strings
for i in range(len(ledges)):
output.append("{},{},{}".format(x[ledges[i]],rebin_step_for_string_output,x[redges[i]]))
return output
# If not, the indices of the edges are returned, NOT the Tof values
elif output_format == "tof":
# Return the tofs as float
for i in range(len(ledges)):
output.extend([x[ledges[i]],x[redges[i]]])
else:
return [ledges, redges]
# If not, the indices of the edges are returned
for i in range(len(ledges)):
output.extend([ledges[i],redges[i]])
return output
################################################################################
################################################################################
......
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