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

Updated window finding so it makes more sense and can now use workspace as input

parent 3eab8d66
......@@ -19,4 +19,4 @@
from .wfm_processing import WFMProcessor, make_frame_shifts, get_frame_shifts
from .data_handling import load, save_fits_stack, normalize
from .get_wfm_windows import detect_peaks, check_peaks, get_wfm_windows
from .get_wfm_windows import get_wfm_windows
from __future__ import division, print_function
import numpy as np
from mantid.simpleapi import *
from scipy.ndimage import gaussian_filter1d
from scipy.signal import find_peaks
......@@ -96,9 +97,9 @@ def check_peaks(peaks, nwindows):
################################################################################
################################################################################
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"):
def get_wfm_windows(workspace=None, data=None, filename=None, nwindows=6,
bg_threshold=0.05, win_threshold=0.3, plot=False, gsmooth=2,
xrange=None, output_format="tof", rebin_parameters=None):
"""
Get the left and right edges of wave-frame multiplication windows.
......@@ -143,11 +144,29 @@ def get_wfm_windows(data=None, filename=None, nwindows=6, bg_threshold=0.05,
['15000,64,20000', '25000,64,30000', ...]
"""
if data is None:
if filename is not None:
data = np.loadtxt(filename)
else:
raise RuntimeError("Either data or filename must be defined!")
if (workspace is None) and (data is None) and (filename is None):
raise RuntimeError("Either workspace, data, or filename must be defined!")
if workspace is not None:
summed = SumSpectra(workspace)
if rebin_parameters is None:
summed_x = summed.extractX()[0]
tof_min = np.amin(summed_x)
tof_max = np.amax(summed_x)
nbins = 512
binsize = (tof_max-tof_min) / float(nbins)
rebin_parameters = '{},{},{}'.format(tof_min, binsize, tof_max)
summed = Rebin(summed, rebin_parameters, PreserveEvents=False)
nz = summed.blocksize()
data = np.zeros([nz,2])
summed_x = summed.extractX()[0]
summed_y = summed.extractY()[0]
data[:,0] = 0.5*(summed_x[1:]+summed_x[:-1])
data[:,1] = summed_y
DeleteWorkspace(summed)
elif filename is not None:
data = np.loadtxt(filename)
nx = np.shape(data)[0]
print(nx)
......@@ -155,8 +174,6 @@ def get_wfm_windows(data=None, filename=None, nwindows=6, bg_threshold=0.05,
y = data[:,1]
# Smooth the data with a gaussian filter
if gsmooth == 0:
gsmooth = max(int(nx/500),1)
if gsmooth is not None:
y = gaussian_filter1d(y, gsmooth)
......@@ -261,26 +278,29 @@ def get_wfm_windows(data=None, filename=None, nwindows=6, bg_threshold=0.05,
# The window edge is the first value that exceeds the a fraction of the mean:
# `y > win_threshold * mean`.
# Compute frame means once
y_minus_background = y - bg_value
frame_mean = []
for p in range(nwindows):
frame_mean.append(np.average(
y_minus_background[peaks[p]:peaks[p+1]][np.where(
y_minus_background[peaks[p]:peaks[p+1]] > 0.0)]))
print(frame_mean[-1])
# Define left and right window edges
ledges = [i_start]
redges = []
y_minus_background = y - bg_value
for p in range(1,len(peaks)-1):
# Towards the right ===================
rmean = np.average(y_minus_background[peaks[p]:peaks[p+1]])
# Find left edge iterating towards the right
# Towards the right: Find left edge iterating towards the right
for i in range(peaks[p]+1,peaks[p+1]):
if (y_minus_background[i] - y_minus_background[peaks[p]]) >= (win_threshold*(rmean-y_minus_background[peaks[p]])):
if y_minus_background[i] >= win_threshold*frame_mean[p]:
ledges.append(i)
break
# Towards the left =======================
lmean = np.average(y_minus_background[peaks[p-1]:peaks[p]])
# Find left edge iterating towards the right
# Towards the left: Find left edge iterating towards the right
for i in range(peaks[p]-1,peaks[p-1],-1):
if (y_minus_background[i] - y_minus_background[peaks[p]]) >= (win_threshold*(lmean-y_minus_background[peaks[p]])):
if y_minus_background[i] >= win_threshold*frame_mean[p-1]:
redges.append(i)
break
......@@ -368,62 +388,75 @@ if __name__ == "__main__":
parser = argparse.ArgumentParser(
description="Automatically detect wave-frame multiplication windows")
parser.add_argument("-w", "--workspace", default=None, dest="workspace",
help="A Mantid workspace containig the WFM data.")
parser.add_argument("-d", "--data", type=str, default=None, dest="data",
help="the 2D data array to process")
help="A 2D data array to process.")
parser.add_argument("-f", "--filename", type=str, default=None, dest="filename",
help="the name of the file to process")
parser.add_argument("-f", "--filename", type=str, default=None,
dest="filename",
help="The name of the file to process.")
parser.add_argument("-n", "--nwindows", type=int, default=6, dest="nwindows",
help="the number of windows to be found")
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",
help="threshold above which we are no longer in"
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).")
parser.add_argument("-w","--win-threshold", type=float, default=0.3, dest="win_threshold",
help="threshold to find window edge, as a percentage of"
"average signal inside window.")
parser.add_argument("-t", "--win-threshold", type=float, default=0.3,
dest="win_threshold",
help="Threshold to find window edge, as a percentage "
"of average signal inside window.")
parser.add_argument("-g", "--gsmooth", type=int, default=0, dest="gsmooth",
help="the width of the Gaussian kernel to smooth the"
"data. If 0 (the default), then a guess is made based"
"on the resolution of the data set. If None, then no"
help="The width of the Gaussian kernel to smooth the "
"data. If 0 (the default), then a guess is made based "
"on the resolution of the data set. If None, then no "
"smoothing is carried out.")
parser.add_argument("-x", "--xrange", default=None, dest="xrange",
help="an array or tuple containing 2 elements which are"
"the start and end TOF values to be considered for"
"analysis. The data outside these bounds is ignored."
"You can set one of the limits to None if you want to"
"have no limit on that end. For example:"
help="An array or tuple containing 2 elements which "
"are the start and end TOF values to be considered "
"for analysis. The data outside these bounds is "
"ignored. You can set one of the limits to None if "
"you want to have no limit on that end. For example: "
"xrange=(1000,55000) or xrange=[5000,None].")
parser.add_argument("-p", "--plot", action="store_true",
help="output the results to a plot")
parser.add_argument("-r", "--rebin-step-for-string-output", type=str, default=None,
dest="rebin_step_for_string_output",
help="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.")
parser.add_argument("-p", "--plot", type=str, default=None,
help="Name of the file for saving a figure showing "
"the frames.")
parser.add_argument("-o", "--output-format", type=str, default=None,
dest="output_format",
help="A string. Can be 3 different things:\n"
" - \"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, ...]\n"
" - \"indices\": returns the indices of the left and"
"right edges as arrays of values.\n"
" - \"rebin:64\": If in this format, 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', ...]")
options = parser.parse_args()
output = get_wfm_windows(data=options.data,
output = get_wfm_windows(workspace=options.workspace,
data=options.data,
filename=options.filename,
nwindows=options.nwindows,
bg_threshold=options.bg_threshold,
win_threshold=options.win_threshold,
gsmooth=options.gsmooth,
xrange=options.xrange,
rebin_step_for_string_output=options.rebin_step_for_string_output,
output_format=options.output_format,
plot=options.plot)
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