Feb 2022
This commit is contained in:
@@ -9,7 +9,7 @@ initialized = False
|
||||
|
||||
|
||||
def initialize(params):
|
||||
global initialized, buffer, device, step_length, edge_type, refinement, dark_event, fel_on_event, use_dark, calib, use_filter, filter_window
|
||||
global initialized, buffer_savgol, device, step_length, edge_type, refinement, dark_event, fel_on_event, use_dark, calib, use_filter, filter_window, buffer
|
||||
|
||||
device = params["device"]
|
||||
step_length = params["step_length"]
|
||||
@@ -17,11 +17,12 @@ def initialize(params):
|
||||
refinement = params["refinement"]
|
||||
dark_event = params["dark_event"]
|
||||
fel_on_event = params["fel_on_event"]
|
||||
buffer = deque(maxlen=params["buffer_length"])
|
||||
buffer_savgol = deque(maxlen=params["buffer_length"])
|
||||
use_dark = params["use_dark"]
|
||||
calib = params["calib"]
|
||||
filter_window = params["filter_window"]
|
||||
use_filter = params['filter']
|
||||
# use_filter = params['filter']
|
||||
buffer = deque(maxlen=params["buffer_length"])
|
||||
initialized = True
|
||||
|
||||
|
||||
@@ -66,25 +67,28 @@ def find_edge(data, step_length=50, edge_type="falling", refinement=1):
|
||||
def process(data, pulse_id, timestamp, params):
|
||||
if not initialized:
|
||||
initialize(params)
|
||||
output = {}
|
||||
|
||||
# Read stream inputs
|
||||
prof_sig = data[params["prof_sig"]]
|
||||
if use_filter:
|
||||
prof_sig = savgol_filter(prof_sig,filter_window,3)
|
||||
prof_sig_savgol = savgol_filter(prof_sig, filter_window, 3)
|
||||
events = data[params["events"]]
|
||||
|
||||
if prof_sig.ndim == 1:
|
||||
prof_sig = prof_sig[np.newaxis, :]
|
||||
|
||||
if events[dark_event] and use_dark:
|
||||
buffer.append(prof_sig)
|
||||
|
||||
if prof_sig_savgol.ndim == 1:
|
||||
prof_sig_savgol = prof_sig_savgol[np.newaxis, :]
|
||||
|
||||
if events[dark_event] and use_dark:
|
||||
buffer_savgol.append(prof_sig_savgol)
|
||||
edge_results = {"edge_pos": np.nan, "xcorr": np.nan, "xcorr_ampl": np.nan, "signal":np.nan}
|
||||
else:
|
||||
if events[fel_on_event] and buffer:
|
||||
prof_sig = prof_sig / np.mean(buffer, axis=0)
|
||||
edge_results = find_edge(prof_sig, step_length, edge_type, refinement)
|
||||
if events[fel_on_event] and buffer_savgol:
|
||||
prof_sig_norm = prof_sig_savgol / np.mean(buffer_savgol, axis=0)
|
||||
edge_results = find_edge(prof_sig_norm, step_length, edge_type, refinement)
|
||||
elif events[fel_on_event] and not use_dark:
|
||||
edge_results = find_edge(prof_sig, step_length, edge_type, refinement)
|
||||
edge_results = find_edge(prof_sig_savgol, step_length, edge_type, refinement)
|
||||
else:
|
||||
edge_results = {"edge_pos": np.nan, "xcorr": np.nan, "xcorr_ampl": np.nan, "signal":np.nan}
|
||||
|
||||
@@ -92,9 +96,27 @@ def process(data, pulse_id, timestamp, params):
|
||||
edge_results["arrival_time"] = edge_results["edge_pos"] * calib
|
||||
|
||||
# Set bs outputs
|
||||
output = {}
|
||||
for key, value in edge_results.items():
|
||||
output[f"{device}:{key}"] = value
|
||||
|
||||
return output
|
||||
output[f"{device}:raw_wf"] = prof_sig
|
||||
output[f"{device}:raw_wf_savgol"] = prof_sig_savgol
|
||||
|
||||
if events[dark_event]:
|
||||
output[f"{device}:dark_wf"] = prof_sig
|
||||
output[f"{device}:dark_wf_savgol"] = prof_sig_savgol
|
||||
else:
|
||||
output[f"{device}:dark_wf"] = np.nan
|
||||
output[f"{device}:dark_wf_savgol"] = np.nan
|
||||
|
||||
if buffer:
|
||||
output[f"{device}:avg_dark_wf"] = np.mean(buffer, axis=0)
|
||||
else:
|
||||
output[f"{device}:avg_dark_wf"] = np.nan
|
||||
|
||||
if buffer_savgol:
|
||||
output[f"{device}:avg_dark_wf_savgol"] = np.mean(buffer_savgol, axis=0)
|
||||
else:
|
||||
output[f"{device}:avg_dark_wf_savgol"] = np.nan
|
||||
|
||||
return output
|
||||
|
||||
@@ -8,6 +8,5 @@ def process_image(image, pulse_id, timestamp, x_axis, y_axis, parameters, bsdata
|
||||
channels = ["intensity","x_center_of_mass","x_fwhm","x_rms","x_fit_amplitude", "x_fit_mean","x_fit_offset","x_fit_standard_deviation","x_profile"]
|
||||
prefix = parameters["camera_name"]
|
||||
for c in channels:
|
||||
ret[prefix+":"+c] = r[c]
|
||||
ret[prefix+":"+c] = r[c]
|
||||
return ret
|
||||
|
||||
|
||||
@@ -6,7 +6,7 @@ from cam_server.utils import create_thread_pvs, epics_lock
|
||||
|
||||
import json
|
||||
|
||||
import numpy
|
||||
import numpy as np
|
||||
import scipy.signal
|
||||
import scipy.optimize
|
||||
import numba
|
||||
@@ -20,50 +20,47 @@ output_pv, center_pv, fwhm_pv, ymin_pv, ymax_pv, axis_pv = None, None, None, Non
|
||||
roi = [0, 0]
|
||||
initialized = False
|
||||
sent_pid = -1
|
||||
nrows = 1
|
||||
axis = None
|
||||
|
||||
|
||||
@numba.njit(parallel=True)
|
||||
@numba.njit(parallel=False)
|
||||
def get_spectrum(image, background):
|
||||
y = image.shape[0]
|
||||
x = image.shape[1]
|
||||
|
||||
profile = numpy.zeros(x, dtype=numpy.uint32)
|
||||
profile = np.zeros(x, dtype=np.float64)
|
||||
|
||||
for i in numba.prange(y):
|
||||
for j in range(x):
|
||||
v = image[i, j]
|
||||
b = background[i, j]
|
||||
if v > b:
|
||||
v -= b
|
||||
else:
|
||||
v = 0
|
||||
profile[j] += v
|
||||
profile[j] += image[i, j] - background[i, j]
|
||||
return profile
|
||||
|
||||
|
||||
def initialize(parameters):
|
||||
def initialize(params):
|
||||
global ymin_pv, ymax_pv, axis_pv, output_pv, center_pv, fwhm_pv
|
||||
global channel_names
|
||||
epics_pv_name_prefix = parameters["camera_name"]
|
||||
output_pv_name = epics_pv_name_prefix + ":SPECTRUM_Y"
|
||||
center_pv_name = epics_pv_name_prefix + ":SPECTRUM_CENTER"
|
||||
fwhm_pv_name = epics_pv_name_prefix + ":SPECTRUM_FWHM"
|
||||
ymin_pv_name = epics_pv_name_prefix + ":SPC_ROI_YMIN"
|
||||
ymax_pv_name = epics_pv_name_prefix + ":SPC_ROI_YMAX"
|
||||
axis_pv_name = epics_pv_name_prefix + ":SPECTRUM_X"
|
||||
channel_names = [output_pv_name, center_pv_name, fwhm_pv_name, ymin_pv_name, ymax_pv_name, axis_pv_name]
|
||||
global channel_names, spectra_buffer
|
||||
camera_name = params["camera_name"]
|
||||
output_pv_name = camera_name + ":SPECTRUM_Y"
|
||||
center_pv_name = camera_name + ":SPECTRUM_CENTER"
|
||||
fwhm_pv_name = camera_name + ":SPECTRUM_FWHM"
|
||||
ymin_pv_name = camera_name + ":SPC_ROI_YMIN"
|
||||
ymax_pv_name = camera_name + ":SPC_ROI_YMAX"
|
||||
axis_pv_name = camera_name + ":SPECTRUM_X"
|
||||
com_pv_name = camera_name + ":SPECTRUM_COM"
|
||||
std_pv_name = camera_name + ":SPECTRUM_STD"
|
||||
channel_names = [output_pv_name, center_pv_name, fwhm_pv_name, ymin_pv_name, ymax_pv_name, axis_pv_name, com_pv_name, std_pv_name]
|
||||
|
||||
|
||||
def process_image(image, pulse_id, timestamp, x_axis, y_axis, parameters, bsdata=None, background=None):
|
||||
global roi, initialized, sent_pid
|
||||
global roi, initialized, sent_pid, nrows, axis
|
||||
global channel_names
|
||||
|
||||
|
||||
if not initialized:
|
||||
initialize(parameters)
|
||||
initialized = True
|
||||
[output_pv, center_pv, fwhm_pv, ymin_pv, ymax_pv, axis_pv] = create_thread_pvs(channel_names)
|
||||
[output_pv, center_pv, fwhm_pv, ymin_pv, ymax_pv, axis_pv, com_pv, std_pv] = create_thread_pvs(channel_names)
|
||||
processed_data = dict()
|
||||
epics_pv_name_prefix = parameters["camera_name"]
|
||||
camera_name = parameters["camera_name"]
|
||||
|
||||
if ymin_pv and ymin_pv.connected:
|
||||
roi[0] = ymin_pv.value
|
||||
@@ -85,12 +82,13 @@ def process_image(image, pulse_id, timestamp, x_axis, y_axis, parameters, bsdata
|
||||
# match the energy axis to image width
|
||||
axis = axis[:image.shape[1]]
|
||||
|
||||
processing_image = image
|
||||
processing_image = image.astype(np.float32) - np.float32(parameters["pixel_bkg"])
|
||||
nrows, ncols = processing_image.shape
|
||||
|
||||
# validate background data if passive mode (background subtraction handled here)
|
||||
background_image = parameters.pop('background_data', None)
|
||||
if isinstance(background_image, numpy.ndarray):
|
||||
if isinstance(background_image, np.ndarray):
|
||||
background_image = background_image.astype(np.float32)
|
||||
if background_image.shape != processing_image.shape:
|
||||
_logger.info("Invalid background shape: %s instead of %s" % (
|
||||
str(background_image.shape), str(processing_image.shape)))
|
||||
@@ -98,7 +96,7 @@ def process_image(image, pulse_id, timestamp, x_axis, y_axis, parameters, bsdata
|
||||
else:
|
||||
background_image = None
|
||||
|
||||
processed_data[epics_pv_name_prefix + ":processing_parameters"] = json.dumps(
|
||||
processed_data[camera_name + ":processing_parameters"] = json.dumps(
|
||||
{"roi": roi, "background": None if (background_image is None) else parameters.get('image_background')})
|
||||
|
||||
# crop the image in y direction
|
||||
@@ -113,7 +111,7 @@ def process_image(image, pulse_id, timestamp, x_axis, y_axis, parameters, bsdata
|
||||
if background_image is not None:
|
||||
spectrum = get_spectrum(processing_image, background_image)
|
||||
else:
|
||||
spectrum = processing_image.sum(0, 'uint32')
|
||||
spectrum = np.sum(processing_image, axis=0)
|
||||
|
||||
# smooth the spectrum with savgol filter with 51 window size and 3rd order polynomial
|
||||
smoothed_spectrum = scipy.signal.savgol_filter(spectrum, 51, 3)
|
||||
@@ -129,26 +127,38 @@ def process_image(image, pulse_id, timestamp, x_axis, y_axis, parameters, bsdata
|
||||
offset, amplitude, center, sigma = functions.gauss_fit_psss(smoothed_spectrum[::2], axis[::2],
|
||||
offset=minimum, amplitude=amplitude, skip=skip, maxfev=20)
|
||||
|
||||
smoothed_spectrum_normed = smoothed_spectrum / np.sum(smoothed_spectrum)
|
||||
spectrum_com = np.sum(axis * smoothed_spectrum_normed)
|
||||
spectrum_std = np.sqrt(np.sum((axis - spectrum_com) ** 2 * smoothed_spectrum_normed))
|
||||
|
||||
# outputs
|
||||
processed_data[epics_pv_name_prefix + ":SPECTRUM_Y"] = spectrum
|
||||
processed_data[epics_pv_name_prefix + ":SPECTRUM_X"] = axis
|
||||
processed_data[epics_pv_name_prefix + ":SPECTRUM_CENTER"] = numpy.float64(center)
|
||||
processed_data[epics_pv_name_prefix + ":SPECTRUM_FWHM"] = numpy.float64(2.355 * sigma)
|
||||
processed_data[camera_name + ":SPECTRUM_Y"] = spectrum
|
||||
processed_data[camera_name + ":SPECTRUM_X"] = axis
|
||||
processed_data[camera_name + ":SPECTRUM_CENTER"] = np.float64(center)
|
||||
processed_data[camera_name + ":SPECTRUM_FWHM"] = np.float64(2.355 * sigma)
|
||||
processed_data[camera_name + ":SPECTRUM_COM"] = spectrum_com
|
||||
processed_data[camera_name + ":SPECTRUM_STD"] = spectrum_std
|
||||
|
||||
|
||||
if epics_lock.acquire(False):
|
||||
try:
|
||||
if pulse_id > sent_pid:
|
||||
sent_pid = pulse_id
|
||||
if output_pv and output_pv.connected:
|
||||
output_pv.put(processed_data[epics_pv_name_prefix + ":SPECTRUM_Y"])
|
||||
#_logger.debug("caput on %s for pulse_id %s", output_pv, pulse_id)
|
||||
try:
|
||||
if pulse_id > sent_pid:
|
||||
sent_pid = pulse_id
|
||||
if output_pv and output_pv.connected:
|
||||
output_pv.put(processed_data[camera_name + ":SPECTRUM_Y"])
|
||||
|
||||
if center_pv and center_pv.connected:
|
||||
center_pv.put(processed_data[epics_pv_name_prefix + ":SPECTRUM_CENTER"])
|
||||
if center_pv and center_pv.connected:
|
||||
center_pv.put(processed_data[camera_name + ":SPECTRUM_CENTER"])
|
||||
|
||||
if fwhm_pv and fwhm_pv.connected:
|
||||
fwhm_pv.put(processed_data[epics_pv_name_prefix + ":SPECTRUM_FWHM"])
|
||||
finally:
|
||||
epics_lock.release()
|
||||
if fwhm_pv and fwhm_pv.connected:
|
||||
fwhm_pv.put(processed_data[camera_name + ":SPECTRUM_FWHM"])
|
||||
|
||||
if com_pv and com_pv.connected:
|
||||
com_pv.put(processed_data[camera_name + ":SPECTRUM_COM"])
|
||||
|
||||
if std_pv and std_pv.connected:
|
||||
std_pv.put(processed_data[camera_name + ":SPECTRUM_STD"])
|
||||
finally:
|
||||
epics_lock.release()
|
||||
|
||||
return processed_data
|
||||
|
||||
|
||||
Reference in New Issue
Block a user