From f5ad640ddcaad8ebc4be74ed61946e71949433f3 Mon Sep 17 00:00:00 2001
From: Cammille Carinan <cammille.carinan@xfel.eu>
Date: Mon, 7 Feb 2022 16:18:26 +0100
Subject: [PATCH] Add hRIXS functions
---
src/toolbox_scs/constants.py | 18 ++++
src/toolbox_scs/detectors/hrixs.py | 133 +++++++++++++++++++++++++++++
2 files changed, 151 insertions(+)
create mode 100644 src/toolbox_scs/detectors/hrixs.py
diff --git a/src/toolbox_scs/constants.py b/src/toolbox_scs/constants.py
index bc412c2..8c70dd5 100644
--- a/src/toolbox_scs/constants.py
+++ b/src/toolbox_scs/constants.py
@@ -66,6 +66,18 @@ mnemonics = {
"M2BEND": ({'source': 'SA3_XTD10_MIRR-2/MOTOR/BENDER',
'key': 'actualPosition.value',
'dim': None},),
+ "hRIXS_det": ({'source': 'SCS_HRIXS_DET/CAM/CAMERA:daqOutput',
+ 'key': 'data.image.pixels',
+ 'dim': ['x', 'y']},),
+ "chem_X": ({'source': 'SCS_CHEM_JET/MOTOR/MANA_X',
+ 'key': 'actualPosition.value',
+ 'dim': None},),
+ "chem_Y": ({'source': 'SCS_CHEM_JET/MOTOR/MANA_Y',
+ 'key': 'actualPosition.value',
+ 'dim': None},),
+ "chem_Z": ({'source': 'SCS_CHEM_JET/MOTOR/MANA_Z',
+ 'key': 'actualPosition.value',
+ 'dim': None},),
"VSLIT": ({'source': 'SA3_XTD10_VSLIT/MDL/BLADE',
'key': 'actualGap.value',
'dim': None},),
@@ -75,6 +87,12 @@ mnemonics = {
"HSLIT": ({'source': 'SCS_XTD10_HSLIT/MDL/BLADE',
'key': 'actualGap.value',
'dim': None},),
+ "H_BLADE_LEFT": ({'source': 'SCS_XTD10_HSLIT/MOTOR/BLADE_LEFT',
+ 'key': 'actualPosition.value',
+ 'dim': None},),
+ "H_BLADE_RIGHT": ({'source': 'SCS_XTD10_HSLIT/MOTOR/BLADE_RIGHT',
+ 'key': 'actualPosition.value',
+ 'dim': None},),
"transmission": ({'source': 'SA3_XTD10_VAC/MDL/GATT_TRANSMISSION_MONITOR',
'key': 'Estimated_Tr.value',
'dim': None},
diff --git a/src/toolbox_scs/detectors/hrixs.py b/src/toolbox_scs/detectors/hrixs.py
new file mode 100644
index 0000000..a66483d
--- /dev/null
+++ b/src/toolbox_scs/detectors/hrixs.py
@@ -0,0 +1,133 @@
+import numpy as np
+from scipy.optimize import curve_fit
+from scipy.signal import fftconvolve
+
+
+# -----------------------------------------------------------------------------
+# Curvature
+
+def find_curvature(image, frangex=None, frangey=None,
+ deg=2, axis=1, offset=100):
+ # Resolve arguments
+ x_range = (0, image.shape[1])
+ if frangex is not None:
+ x_range = (max(frangex[0], x_range[0]), min(frangex[1], x_range[1]))
+ y_range = (0, image.shape[0])
+ if frangex is not None:
+ y_range = (max(frangey[0], y_range[0]), min(frangey[1], y_range[1]))
+
+ axis_range = y_range if axis == 1 else x_range
+ axis_dim = image.shape[axis - 1]
+
+ # Get kernel
+ integral = image[slice(*y_range), slice(*x_range)].mean(axis=axis)
+ roi = np.ones([axis_range[1] - axis_range[0], axis_dim])
+ ref = roi * integral[:, np.newaxis]
+
+ # Get sliced image
+ slice_ = [slice(None), slice(None)]
+ slice_[axis - 1] = slice(max(axis_range[0] - offset, 0),
+ min(axis_range[1] + offset, axis_dim))
+ sliced = image[tuple(slice_)]
+ if axis == 0:
+ sliced = sliced.T
+
+ # Get curvature factor from cross correlation
+ crosscorr = fftconvolve(sliced,
+ ref[::-1, :],
+ axes=0, )
+ shifts = np.argmax(crosscorr, axis=0)
+ curv = np.polyfit(np.arange(axis_dim), shifts, deg=deg)
+ return curv[:-1][::-1]
+
+
+def correct_curvature(image, factor=None, axis=1):
+ if factor is None:
+ return
+
+ if axis == 1:
+ image = image.T
+
+ ydim, xdim = image.shape
+ x = np.arange(xdim + 1)
+ y = np.arange(ydim + 1)
+ xx, yy = np.meshgrid(x[:-1] + 0.5, y[:-1] + 0.5)
+ xxn = xx - factor[0] * yy - factor[1] * yy ** 2
+ ret = np.histogramdd((xxn.flatten(), yy.flatten()),
+ bins=[x, y],
+ weights=image.flatten())[0]
+
+ return ret if axis == 1 else ret.T
+
+
+def get_spectrum(image, cal_factor=None, axis=0,
+ pixel_range=None, energy_range=None, ):
+ start, stop = (0, image.shape[axis - 1])
+ if pixel_range is not None:
+ start = max(pixel_range[0] or start, start)
+ stop = min(pixel_range[1] or stop, stop)
+
+ edge = image.sum(axis=axis)[start:stop]
+ bins = np.arange(start, stop + 1)
+ centers = (bins[1:] + bins[:-1]) * 0.5
+ if cal_factor is not None:
+ centers, edge = calibrate(centers, edge,
+ factor=cal_factor,
+ range_=energy_range)
+
+ return centers, edge
+
+
+# -----------------------------------------------------------------------------
+# Energy calibration
+
+
+def energy_calibration(channels, energies):
+ return np.polyfit(channels, energies, deg=1)
+
+
+def calibrate(x, y=None, factor=None, range_=None):
+ if factor is not None:
+ x = np.polyval(factor, x)
+
+ if y is not None and range_ is not None:
+ start = np.argmin(np.abs((x - range_[0])))
+ stop = np.argmin(np.abs((x - range_[1])))
+ # Calibrated energies have a different direction
+ x, y = x[stop:start], y[stop:start]
+
+ return x, y
+
+
+# -----------------------------------------------------------------------------
+# Gaussian-related functions
+
+
+FWHM_COEFF = 2 * np.sqrt(2 * np.log(2))
+
+
+def gaussian_fit(x_data, y_data, offset=0):
+ """
+ Centre-of-mass and width. Lifted from image_processing.imageCentreofMass()
+ """
+
+ x0 = np.average(x_data, weights=y_data)
+ sx = np.sqrt(np.average((x_data - x0) ** 2, weights=y_data))
+
+ # Gaussian fit
+ baseline = y_data.min()
+ p_0 = (y_data.max(), x0 + offset, sx, baseline)
+ try:
+ p_f, _ = curve_fit(gauss1d, x_data, y_data, p_0, maxfev=10000)
+ return p_f
+ except (RuntimeError, TypeError) as e:
+ print(e)
+ return None
+
+
+def gauss1d(x, height, x0, sigma, offset):
+ return height * np.exp(-0.5 * ((x - x0) / sigma) ** 2) + offset
+
+
+def to_fwhm(sigma):
+ return abs(sigma * FWHM_COEFF)
--
GitLab