diff --git a/src/toolbox_scs/detectors/hrixs.py b/src/toolbox_scs/detectors/hrixs.py
index 31f6984f3b653e470d53e4130c9c19a682d192c4..aa40480e53977097d437fe68464fe9daa5512955 100644
--- a/src/toolbox_scs/detectors/hrixs.py
+++ b/src/toolbox_scs/detectors/hrixs.py
@@ -196,6 +196,50 @@ def _esrf_centroid(image, threshold=THRESHOLD, curvature=(CURVE_A, CURVE_B)):
res.append((my, mx))
return res
+def _esrf_centroid_double(image, energy=700, avoid_dbl=False):
+
+ # Multiplication factor * ADU/photon
+ SpotLOW = 0.2 * energy/3.6/1.06
+ SpotHIGH = 1.5 * energy/3.6/1.06
+ low_th_px = 0.2 * energy/3.6/1.06
+ high_th_px = 1 * energy/3.6/1.06
+
+ if avoid_dbl:
+ SpotHIGH = 100000
+ print("Double hit events are not taken into account\n")
+
+ img = image-image.mean()
+
+ gs = 2
+
+ # Find potential hits on a clipped image where 2 rows/columns are excluded
+ # from the edges because centroiding can't be done at the edge
+ cp = np.argwhere((img[gs//2 : -gs//2, gs//2 : -gs//2] > low_th_px)*
+ (img[gs//2 : -gs//2, gs//2 : -gs//2] < high_th_px))+np.array([gs//2, gs//2])
+
+ if len(cp) > 100000:
+ raise RuntimeError('Threshold to low or acquisition time to long')
+
+ res = []
+ double = []
+ for cy, cx in cp:
+ spot = img[cy - gs//2 : cy + gs//2 + 1, cx - gs//2 : cx + gs//2 + 1]
+ #spot = img[cy - gs//2 : cy + gs//2 + 1, cx - gs//2 : cx + gs//2 + 1]
+ #if spot.sum() < 2 * (THRESHOLD - base):
+ # continue
+ spot[spot < 0] = 0
+ if (spot > img[cy, cx]).sum() == 0:
+ mx = np.average(np.arange(cx - gs//2, cx + gs//2 + 1), weights=spot.sum(axis=0))
+ my = np.average(np.arange(cy - gs//2, cy + gs//2 + 1), weights=spot.sum(axis=1))
+ my -= (CURVE_A + CURVE_B * mx) * mx
+ if (spot.sum()>SpotLOW) and (spot.sum()<=SpotHIGH):
+ res.append((my, mx))
+ elif (spot.sum()>SpotHIGH):
+ res.append((my, mx))
+ res.append((my, mx))
+ double.append((my, mx))
+ return res, double
+
def _new_centroid(image, threshold=THRESHOLD, std_threshold=3.5, curvature=(CURVE_A, CURVE_B)):
"""find the position of photons with sub-pixel precision
@@ -223,6 +267,7 @@ def _new_centroid(image, threshold=THRESHOLD, std_threshold=3.5, curvature=(CURV
"too many peaks, threshold too low or acquisition time too high")
res = []
+ double = []
for cy, cx in cp:
spot = image[cy: cy + 4, cx: cx + 4] - base
mx = np.average(np.arange(cx, cx + 4), weights=spot.sum(axis=0))
@@ -231,9 +276,58 @@ def _new_centroid(image, threshold=THRESHOLD, std_threshold=3.5, curvature=(CURV
res.append((my, mx))
return res
+def _new_centroid_double(image, threshold=THRESHOLD, std_threshold=3.5,
+ dbl_threshold=1.5, curvature=(CURVE_A, CURVE_B)):
+ """
+ A photon is supposed to have hit the detector if the intensity within a
+ 2-by-2 square exceeds a threshold. In this case the position of the photon
+ is calculated as the center-of-mass in a 4-by-4 square.
-centroid = _new_centroid
+ Return the list of x,y coordinate pairs, corrected by the curvature.
+ This implementation includes standard energy dependent threshholds and
+ checks for double hits. If a double hit is found two instances of the
+ photon are returned and the double hits are also returned
+ """
+
+ # Implement energy dependent threshhold later
+ #SpotLOW = 0.2 * energy/3.6/1.06
+ #SpotHIGH = 1.5 * energy/3.6/1.06
+ #low_th_px = 0.2 * energy/3.6/1.06
+ #high_th_px = 1 * energy/3.6/1.06
+
+
+ base = image.mean()
+ corners = image[1:, 1:] + image[:-1, 1:] + image[1:, :-1] + image[:-1, :-1]
+ if threshold is None:
+ threshold = corners.mean() + std_threshold * corners.std()
+ SpotHIGH=1.5*threshold
+ middle = corners[1:-1, 1:-1]
+ candidates = (
+ (middle > threshold)
+ * (middle >= corners[:-2, 1:-1]) * (middle > corners[2:, 1:-1])
+ * (middle >= corners[1:-1, :-2]) * (middle > corners[1:-1, 2:])
+ * (middle >= corners[:-2, :-2]) * (middle > corners[2:, :-2])
+ * (middle >= corners[:-2, 2:]) * (middle > corners[2:, 2:]))
+ cp = np.argwhere(candidates)
+ if len(cp) > 10000:
+ raise RuntimeError(
+ "too many peaks, threshold too low or acquisition time too high")
+
+ res = []
+ dres = []
+ for cy, cx in cp:
+ spot = image[cy: cy + 4, cx: cx + 4] - base
+ mx = np.average(np.arange(cx, cx + 4), weights=spot.sum(axis=0))
+ my = np.average(np.arange(cy, cy + 4), weights=spot.sum(axis=1))
+ my -= (curvature[0] + curvature[1] * mx) * mx
+ if spot.sum()<=SpotHIGH:
+ res.append((my, mx))
+ else:
+ res.append((my, mx))
+ res.append((my, mx))
+ dres.append((my,mx))
+ return res, dres
def decentroid(res):
res = np.array(res)
@@ -286,6 +380,7 @@ class hRIXS:
# centroid
THRESHOLD = None # pixel counts above which a hit candidate is assumed
STD_THRESHOLD = 3.5 # same as THRESHOLD, in standard deviations
+ DBL_THRESHOLD = 1.5 # threshhold for doubles if not given
CURVE_A = CURVE_A # curvature parameters as determined elsewhere
CURVE_B = CURVE_B
@@ -313,7 +408,9 @@ class hRIXS:
ENERGY_INTERCEPT = 0
ENERGY_SLOPE = 1
- FIELDS = ['hRIXS_det', 'hRIXS_index', 'hRIXS_delay', 'hRIXS_norm']
+
+ # Load hrixs detector details + mono energy for automatic threshold value
+ FIELDS = ['hRIXS_det', 'hRIXS_index', 'hRIXS_delay', 'hRIXS_norm', 'nrj']
def set_params(self, **params):
for key, value in params.items():
@@ -418,7 +515,35 @@ class hRIXS:
def centroid(self, data, bins=None, return_hits=False, fit_elastic=False,
- hv=None, fit_delta=None):
+ hv=None, fit_delta=None, algorithm='martin2'):
+ """
+ Carry out hit finding and bin them on grid
+
+ Parameters
+ ----------
+ data: xarray.DataArray
+ container containing run data and images to process
+ bins: int
+ number of bins to use in histogram of hits
+ return_hits: bool
+ option to return list of hits or only return the binned spectrum
+ fit_elastic: bool
+ option to fit the elastic line in each image and stores the energy
+ and position of the elastic line in the data container
+ hv: float
+ Incident energy
+ fit_delta: int
+ Width of the fit window in bins
+ algorithm: str
+ one of 'esrf', 'esrf2', 'martin', and 'martin2' where martin2 and
+ esrf2 account for double hits
+
+ Returns
+ -------
+ data: xarray.DataArray
+ container containing data and centroided spectra
+
+ """
#*************************************************************
# Carry out hit finding on data and bin them in grid
# Allows for
@@ -432,6 +557,7 @@ class hRIXS:
#*************************************************************
# If new bins are given, use them
#*************************************************************
+ print('message')
if bins is None:
bins = self.BINS
#*************************************************************
@@ -440,7 +566,14 @@ class hRIXS:
hit_x = []
hit_y = []
hits = []
+ if algorithm == 'martin2' or algorithm=='esrf2':
+ dbl_x = []
+ dbl_y = []
+ dbl_hits = []
+ do_doubles = True
+
ret = np.zeros((len(data["hRIXS_det"]), bins))
+ retd = np.zeros((len(data["hRIXS_det"]), bins))
#*************************************************************
# If the 'actual' photon energy is desired, we'll need to know
# where to start looking for the elastic line
@@ -451,7 +584,7 @@ class hRIXS:
if fit_delta == None:
fit_delta = self.HV_FIT_DELTA
if hv == None:
- # ?? We'll make a guess later
+ #ave_nrj = data.nrj.mean().mean().data
hv = -1 * np.ones(len(data["hRIXS_det"]))
else:
# hv is given as one value or a train-vector
@@ -464,7 +597,7 @@ class hRIXS:
#*************************************************************
# Handle each Aq image separately
#*************************************************************
- for i, (image, r) in enumerate(zip(data["hRIXS_det"], ret)):
+ for i, (image, r, rr) in enumerate(zip(data["hRIXS_det"], ret, retd)):
use_image = image.to_numpy()
#*************************************************************
# Treat background by optionally
@@ -494,14 +627,34 @@ class hRIXS:
#*************************************************************
# Run centroiding on the preprocessed image
#*************************************************************
- c = centroid(
- use_image[self.X_RANGE, self.Y_RANGE].T,
- threshold=self.THRESHOLD,
- std_threshold=self.STD_THRESHOLD,
- curvature=(self.CURVE_A, self.CURVE_B))
+ if algorithm == 'martin':
+ c = _new_centroid(
+ use_image[self.X_RANGE, self.Y_RANGE].T,
+ threshold=self.THRESHOLD,
+ std_threshold=self.STD_THRESHOLD,
+ curvature=(self.CURVE_A, self.CURVE_B))
+ elif algorithm == 'martin2':
+ c,d = _new_centroid_double(
+ use_image[self.X_RANGE, self.Y_RANGE].T,
+ threshold=self.THRESHOLD,
+ std_threshold=self.STD_THRESHOLD,
+ dbl_threshold=self.DBL_THRESHOLD,
+ curvature=(self.CURVE_A, self.CURVE_B))
+ elif algorithm=='esrf':
+ c = _esrf_centroid(
+ use_image[self.X_RANGE, self.Y_RANGE].T,
+ threshold=self.THRESHOLD,
+ curvature=(self.CURVE_A, self.CURVE_B))
+ elif algorithm=='esrf2':
+ ave_nrj = data.nrj.mean().mean().data
+ c,d = _esrf_centroid_double(
+ use_image[self.X_RANGE, self.Y_RANGE].T,
+ energy=ave_nrj,avoid_dbl=False)
+
if not len(c):
continue
rc = np.array(c)
+
#*************************************************************
# If hits have been requested, append the hit data of the
# image to the lists of hit lists
@@ -510,14 +663,37 @@ class hRIXS:
hit_x.append(rc[:, 0])
hit_y.append(rc[:, 1])
hits.append(rc)
- #*************************************************************
- # Assign the spectrum to the spectrum matrix ret. Iteration
- # variable r points to the proper location of ret.
- #*************************************************************
+
+ # Store histogram of hits in ret
hy, hx = np.histogram(
rc[:, 0], bins=bins,
range=(0, self.Y_RANGE.stop - self.Y_RANGE.start))
r[:] = hy
+
+ if do_doubles:
+ # Check to make sure there are actually double hits.
+ # Make empty arrays if not, so hit lits still work
+ if not len(d):
+ print("No Double hits found in image number ",str(i))
+ doubles_found = False
+ if return_hits:
+ dbl_x.append([])
+ dbl_y.append([])
+ dbl_hits.append([])
+ else:
+ doubles_found = True
+ rd = np.array(d)
+ if return_hits:
+ dbl_x.append(rd[:,0])
+ dbl_y.append(rd[:,1])
+ dbl_hits.append(rd)
+ # Make histogram of doubles hits store in retd
+ dy, dx = np.histogram(
+ rd[:,0], bins=bins,
+ range=(0, self.Y_RANGE.stop - self.Y_RANGE.start))
+ rr[:]= dy
+
+
#*************************************************************
# Return the found elastic peak energy, if requested
#*************************************************************
@@ -550,6 +726,11 @@ class hRIXS:
data = data.assign(hits=(("trainId"), hits),
xhits=(("trainId"), hit_x),
yhits=(("trainId"), hit_y))
+ if do_doubles:
+ data = data.assign(dbl_hits=(("trainId"), dbl_hits),
+ dblx=(("trainId"), dbl_x),
+ dbly=(("trainId"), dbl_y))
+
#**********************************************
# If hv was fitted, return it
#**********************************************
@@ -560,6 +741,8 @@ class hRIXS:
# Always assign the spectrum to data
#**********************************************
data = data.assign(spectrum=(("trainId", "energy"), ret))
+ if do_doubles:
+ data = data.assign(dbl_spectrum=(("trainId", "energy"), retd))
return data