Merge pull request #264 from ISISComputingGroup/polref

Reflectometry: add flood-map correction & change `Guassian` fit to use centre-of-mass as an initial `x0` guess

Author: rerpha

doc/plans/reflectometry/detector_mapping_alignment.md

@@ -9,6 +9,10 @@ The plans in this module expect:
 {py:obj}`ibex_bluesky_core.devices.simpledae.PeriodSpecIntegralsReducer`.
 - An angle map, as a `numpy` array of type `float64`, which has the same dimensionality as the set of selected detectors. This
 maps each configured detector pixel to its angular position.
+- An optional flood map, as a {external+scipp:py:obj}`scipp.Variable`. This should have a dimension label of "spectrum"
+and have the same dimensionality as the set of selected detectors. This array may have variances. This is used to
+normalise pixel efficiencies: raw counts are divided by the flood to get scaled counts. If no flood map is provided, no
+normalisation will be performed.
 
 ## Angle scan
 

src/ibex_bluesky_core/callbacks/__init__.py

@@ -229,7 +229,11 @@ def start(self, doc: RunStart) -> None:
                     # where a fit result can be returned before
                     # the QtAwareCallback has had a chance to process it.
                     self._subs.append(self._live_fit)
-                self._subs.append(LiveFitPlot(livefit=self._live_fit, ax=ax))
+
+                # Sample 5000 points as this strikes a reasonable balance between displaying
+                # 'enough' points for almost any scan (even after zooming in on a peak), while
+                # not taking 'excessive' compute time to generate these samples.
+                self._subs.append(LiveFitPlot(livefit=self._live_fit, ax=ax, num_points=5000))
             else:
                 self._subs.append(self._live_fit)
 

src/ibex_bluesky_core/callbacks/_plotting.py

@@ -206,7 +206,7 @@ def __init__(
         y: str,
         ax: Axes,
         postfix: str,
-        output_dir: str | os.PathLike[str] | None,
+        output_dir: str | os.PathLike[str] | None = None,
     ) -> None:
         """Initialise the PlotPNGSaver callback.
 

src/ibex_bluesky_core/callbacks/reflectometry/_det_map.py

@@ -46,12 +46,15 @@ class DetMapHeightScanLiveDispatcher(LiveDispatcher):
     monitor spectrum used for normalization).
     """
 
-    def __init__(self, *, mon_name: str, det_name: str, out_name: str) -> None:
+    def __init__(
+        self, *, mon_name: str, det_name: str, out_name: str, flood: sc.Variable | None = None
+    ) -> None:
         """Init."""
         super().__init__()
         self._mon_name = mon_name
         self._det_name = det_name
         self._out_name = out_name
+        self._flood = flood if flood is not None else sc.scalar(value=1, dtype="float64")
 
     def event(self, doc: Event, **kwargs: dict[str, Any]) -> Event:
         """Process an event."""
@@ -62,6 +65,8 @@ def event(self, doc: Event, **kwargs: dict[str, Any]) -> Event:
         det = sc.Variable(dims=["spectrum"], values=det_data, variances=det_data, dtype="float64")
         mon = sc.Variable(dims=["spectrum"], values=mon_data, variances=mon_data, dtype="float64")
 
+        det /= self._flood
+
         det_sum = det.sum()
         mon_sum = mon.sum()
 
@@ -90,19 +95,31 @@ class DetMapAngleScanLiveDispatcher(LiveDispatcher):
     """
 
     def __init__(
-        self, x_data: npt.NDArray[np.float64], x_name: str, y_in_name: str, y_out_name: str
+        self,
+        x_data: npt.NDArray[np.float64],
+        x_name: str,
+        y_in_name: str,
+        y_out_name: str,
+        flood: sc.Variable | None = None,
     ) -> None:
         """Init."""
         super().__init__()
         self.x_data = x_data
         self.x_name = x_name
 
-        self.y_data = np.zeros_like(x_data)
+        self.y_data = sc.array(
+            dims=["spectrum"],
+            values=np.zeros_like(x_data),
+            variances=np.zeros_like(x_data),
+            dtype="float64",
+        )
         self.y_in_name: str = y_in_name
         self.y_out_name: str = y_out_name
 
         self._descriptor_uid: str | None = None
 
+        self._flood = flood if flood is not None else sc.scalar(value=1, dtype="float64")
+
     def descriptor(self, doc: EventDescriptor) -> None:
         """Process a descriptor."""
         self._descriptor_uid = doc["uid"]
@@ -118,7 +135,10 @@ def event(self, doc: Event, **kwargs: dict[str, Any]) -> Event:
                 f"Shape of data ({data.shape} does not match x_data.shape ({self.x_data.shape})"
             )
 
-        self.y_data += data
+        scaled_data = (
+            sc.array(dims=["spectrum"], values=data, variances=data, dtype="float64") / self._flood
+        )
+        self.y_data += scaled_data
         return doc
 
     def stop(self, doc: RunStop, _md: dict[str, Any] | None = None) -> None:
@@ -128,13 +148,13 @@ def stop(self, doc: RunStop, _md: dict[str, Any] | None = None) -> None:
             return super().stop(doc, _md)
 
         current_time = time.time()
-        for x, y in zip(self.x_data, self.y_data, strict=True):
+        for x, y in zip(self.x_data, self.y_data, strict=True):  # type: ignore (pyright doesn't understand scipp)
             logger.debug("DetMapAngleScanLiveDispatcher emitting event with x=%f, y=%f", x, y)
             event = {
                 "data": {
                     self.x_name: x,
-                    self.y_out_name: y,
-                    self.y_out_name + "_err": np.sqrt(y + 0.5),
+                    self.y_out_name: y.value,
+                    self.y_out_name + "_err": np.sqrt(y.variance + 0.5),
                 },
                 "timestamps": {
                     self.x_name: current_time,

src/ibex_bluesky_core/fitting.py

@@ -1,5 +1,6 @@
 """Fitting methods used by the LiveFit callback."""
 
+import math
 from abc import ABC, abstractmethod
 from collections.abc import Callable
 
@@ -102,6 +103,19 @@ def fit(cls, *args: int) -> FitMethod:
         return FitMethod(model=cls.model(*args), guess=cls.guess(*args))
 
 
+def _guess_cen_and_width(
+    x: npt.NDArray[np.float64], y: npt.NDArray[np.float64]
+) -> tuple[float, float]:
+    """Guess the center and width of a positive peak."""
+    com, total_area = center_of_mass_of_area_under_curve(x, y)
+    y_range = np.max(y) - np.min(y)
+    if y_range == 0.0:
+        width = (np.max(x) - np.min(x)) / 2
+    else:
+        width = total_area / y_range
+    return com, width
+
+
 class Gaussian(Fit):
     """Gaussian Fitting."""
 
@@ -130,8 +144,9 @@ def guess(
         def guess(
             x: npt.NDArray[np.float64], y: npt.NDArray[np.float64]
         ) -> dict[str, lmfit.Parameter]:
-            mean = np.sum(x * y) / np.sum(y)
-            sigma = np.sqrt(np.sum(y * (x - mean) ** 2) / np.sum(y))
+            cen, width = _guess_cen_and_width(x, y)
+            sigma = width / math.sqrt(2 * math.pi)  # From expected area under gaussian
+
             background = np.min(y)
 
             if np.max(y) > abs(np.min(y)):
@@ -142,7 +157,7 @@ def guess(
             init_guess = {
                 "amp": lmfit.Parameter("amp", amp),
                 "sigma": lmfit.Parameter("sigma", sigma, min=0),
-                "x0": lmfit.Parameter("x0", mean),
+                "x0": lmfit.Parameter("x0", cen),
                 "background": lmfit.Parameter("background", background),
             }
 
@@ -547,19 +562,6 @@ def guess(
         return guess
 
 
-def _guess_cen_and_width(
-    x: npt.NDArray[np.float64], y: npt.NDArray[np.float64]
-) -> tuple[float, float]:
-    """Guess the center and width of a positive peak."""
-    com, total_area = center_of_mass_of_area_under_curve(x, y)
-    y_range = np.max(y) - np.min(y)
-    if y_range == 0.0:
-        width = (np.max(x) - np.min(x)) / 2
-    else:
-        width = total_area / y_range
-    return com, width
-
-
 class TopHat(Fit):
     """Top Hat Fitting."""