hRIXS functions

src/toolbox_scs/detectors/hrixs.py

@@ -156,7 +156,7 @@ CURVE_A = 2.19042931e-02  # curvature parameters as determined elsewhere
 CURVE_B = -3.02191568e-07
 
 
-def centroid(image, threshold=THRESHOLD, curvature=(CURVE_A, CURVE_B)):
+def _esrf_centroid(image, threshold=THRESHOLD, curvature=(CURVE_A, CURVE_B)):
     gs = 2
     base = image.mean()
     cp = np.argwhere(image[gs // 2: -gs // 2, gs // 2: -gs // 2] > threshold) + np.array([gs // 2, gs // 2])
@@ -175,6 +175,43 @@ def centroid(image, threshold=THRESHOLD, curvature=(CURVE_A, CURVE_B)):
     return res
 
 
+def _new_centroid(image, threshold=THRESHOLD, curvature=(CURVE_A, CURVE_B)):
+    """find the position of photons with sub-pixel precision
+
+    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.
+
+    Return the list of x,y coordinate pairs, corrected by the curvature.
+    """
+    base = image.mean()
+    corners = image[1:, 1:] + image[:-1, 1:] + image[1:, :-1] + image[:-1, :-1]
+    threshold = corners.mean() + 3.5 * corners.std()
+    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 = []
+    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
+        res.append((my, mx))
+    return res
+
+
+centroid = _new_centroid
+
+
 def decentroid(res):
     res = np.array(res)
     ret = np.zeros(shape=(res.max(axis=0) + 1).astype(int))
@@ -270,12 +307,10 @@ class hRIXS:
         if proposal is None:
             proposal = cls.PROPOSAL
 
-        run, data = tb.load(proposal,
-                            runNB=runNB,
-                            fields=['hRIXS_det', 'SCS_slowTrain'])
+        run, data = tb.load(proposal, runNB=runNB, fields=['hRIXS_det'])
 
         # Get slow train data
-        mnemo = tb.mnemonics['SCS_slowTrain'][0]
+        mnemo = tb.mnemonics_for_run(run)['SCS_slowTrain']
         slow_train = run[mnemo['source'], mnemo['key']].ndarray().sum()
 
         return cls(images=data['hRIXS_det'][1 if first_wrong else 0:].data,