Implement AngleMappingReducer to support Psi scanning on reflectometers with linear detectors (OFFSPEC, POLREF)

doc/plans/reflectometry/psi_scan_linear_detector.md

@@ -0,0 +1,41 @@
+# Psi scanning
+
+On a reflectometer, the Psi parameter is a rotation of the sample, side-to-side as viewed from a neutron's perspective.
+The Psi parameter causes the reflected beam profile, which can be approximated as an elongated ellipse, to be rotated
+when viewed at the detector. Reflectometers want to align the beam such that the beam profile is parallel to detector
+pixels (or, in other words, a horizontal line as seen on the detector).
+
+Psi does not get _redefined_ to zero when the beam is level, but reflectometers do need to know what value of Psi
+_causes_ the beam to be level. The value of Psi which causes the beam to be level will be different for each sample,
+because of differences in how each sample is mounted.
+
+## Linear detectors
+
+This section applies to instruments like OFFSPEC and POLREF, which have a 1-D bank of detectors. These pixels are much
+wider than they are tall; representative pixel sizes are `0.5 x 50mm` on OFFSPEC or `0.5 x 30mm` on POLREF.
+
+A linear pixel configuration means that as Psi is varied, the beam profile reflected off the sample(dark blue) hits a
+different number of detector pixels (light orange):
+
+![Psi not equal to zero - broad peak](psi_scan_linear/psi_neq_0.png)
+
+![Psi equal to zero - narrow peak](psi_scan_linear/psi_eq_0.png)
+
+For a linear detector, the optimum value of Psi occurs when the beam profile is focused onto the smallest number of
+pixels (when the beam profile is parallel to the pixels).
+
+A Psi scan on a linear detector is therefore implemented using the following steps:
+- Physically scan over Psi
+- At each scan point (implemented by {py:obj}`~ibex_bluesky_core.devices.reflectometry.AngleMappingReducer`:
+  - Count & acquire a spectrum-data map describing the counts in all detector pixels
+  - Filter to a range of interesting pixels (for example using {py:obj}`~ibex_bluesky_core.utils.centred_pixel`)
+  - Map those pixels to relative angle, to account for pixel spacing which may not be exactly even between all pixels
+  - Optionally divide counts in each pixel by a flood map, to account for different pixel efficiencies
+  - Fit angle against counts using a {py:obj}`~ibex_bluesky_core.fitting.Gaussian` model
+  - Return the fit parameters, including {py:obj}`~~ibex_bluesky_core.devices.reflectometry.AngleMappingReducer.sigma`
+(width) as the data from each scan point
+- Plotting the returned sigma (width) parameter at each scan point against the scanned variable, Psi, the optimum value
+of Psi occurs when the width is minimised
+  - This is performed using an _outer_ {py:obj}`~ibex_bluesky_core.fitting.Gaussian` fit, using standard
+[callbacks and fitting](/callbacks/fitting/fitting) infrastructure. This is expected to give a negative
+{py:obj}`~ibex_bluesky_core.fitting.Gaussian` curve, where the `x0` parameter describes the optimum value of Psi

src/ibex_bluesky_core/devices/reflectometry.py

@@ -3,23 +3,31 @@
 import asyncio
 import logging
 
+import numpy as np
+import numpy.typing as npt
+import scipp as sc
 from bluesky.protocols import NamedMovable
 from ophyd_async.core import (
     AsyncStatus,
+    Device,
     SignalR,
     SignalW,
     StandardReadable,
     StandardReadableFormat,
     observe_value,
+    soft_signal_r_and_setter,
 )
 from ophyd_async.epics.core import epics_signal_r, epics_signal_rw, epics_signal_w
 
 from ibex_bluesky_core.devices import NoYesChoice
+from ibex_bluesky_core.devices.dae import Dae
+from ibex_bluesky_core.devices.simpledae import Reducer
+from ibex_bluesky_core.fitting import Gaussian
 from ibex_bluesky_core.utils import get_pv_prefix
 
 logger = logging.getLogger(__name__)
 
-__all__ = ["ReflParameter", "ReflParameterRedefine", "refl_parameter"]
+__all__ = ["AngleMappingReducer", "ReflParameter", "ReflParameterRedefine", "refl_parameter"]
 
 
 class ReflParameter(StandardReadable, NamedMovable[float]):
@@ -125,3 +133,145 @@ def refl_parameter(
     return ReflParameter(
         prefix=prefix, name=name, changing_timeout_s=changing_timeout_s, has_redefine=has_redefine
     )
+
+
+class AngleMappingReducer(Reducer, StandardReadable):
+    """Reflectometry angle-mapping reducer."""
+
+    def __init__(
+        self,
+        *,
+        detectors: npt.NDArray[np.int32],
+        angle_map: npt.NDArray[np.float64],
+        flood: sc.Variable | None = None,
+    ) -> None:
+        """Angle-mapping reducer describing parameters of beam on a 1-D angular detector.
+
+        This :py:obj:`~ibex_bluesky_core.devices.simpledae.Reducer` fits the
+        counts on each pixel of a detector, against the relative angular positions
+        of those pixels. It then exposes fitted quantities, and their standard deviations,
+        as signals from this reducer.
+
+        The fitting model used is a :py:obj:`~ibex_bluesky_core.fitting.Gaussian`.
+        Uncertainties from the counts data are taken into account on these fits -
+        the variances are set to `counts + 0.5`
+        (see :doc:`ADR5 </architectural_decisions/005-variance-addition>`
+        for justification).
+
+        Optionally, a flood map can be provided to normalise for pixel efficiencies
+        before fitting. This flood map is provided as a :py:obj:`scipp.Variable`,
+        which means that it may contain variances.
+
+        .. warning::
+
+            The uncertainties (signals ending in ``_err``) from this reducer are
+            obtained from ``lmfit``. The exact method is described in
+            lmfit's :py:obj:`~lmfit.minimizer.MinimizerResult` documentation.
+            For a perfect fit, which might result from fitting a limited number
+            of points or flat data, uncertainties may be zero.
+
+            Because the uncertainties may be zero, using them in an 'outer' fit is
+            discouraged; an uncertainty of zero implies infinite weighting,
+            which will likely cause the outer fit to fail to converge.
+
+        Args:
+            detectors: numpy array of detector spectra to include.
+            angle_map: numpy array of relative pixel angles for each
+                selected detector
+            flood: Optional :py:obj:`scipp.Variable` describing a flood-correction.
+                This array should be aligned along a "spectrum" dimension; counts are
+                divided by this array before being used in fits. This is used to
+                normalise the intensities detected by each detector pixel.
+
+        """
+        self.amp, self._amp_setter = soft_signal_r_and_setter(float, 0.0)
+        """Amplitude of fitted Gaussian"""
+        self.amp_err, self._amp_err_setter = soft_signal_r_and_setter(float, 0.0)
+        """Amplitude standard deviation of fitted Gaussian.
+
+        This is the standard error reported by :py:obj:`lmfit.model.Model.fit`.
+        """
+        self.sigma, self._sigma_setter = soft_signal_r_and_setter(float, 0.0)
+        """Width (sigma) of fitted Gaussian"""
+        self.sigma_err, self._sigma_err_setter = soft_signal_r_and_setter(float, 0.0)
+        """Width (sigma) standard deviation of fitted Gaussian.
+
+        This is the standard error reported by :py:obj:`lmfit.model.Model.fit`.
+        """
+        self.x0, self._x0_setter = soft_signal_r_and_setter(float, 0.0)
+        """Centre (x0) of fitted Gaussian"""
+        self.x0_err, self._x0_err_setter = soft_signal_r_and_setter(float, 0.0)
+        """Centre (x0) standard deviation of fitted Gaussian.
+
+        This is the standard error reported by :py:obj:`lmfit.model.Model.fit`.
+        """
+        self.background, self._background_setter = soft_signal_r_and_setter(float, 0.0)
+        """Background of fitted Gaussian"""
+        self.background_err, self._background_err_setter = soft_signal_r_and_setter(float, 0.0)
+        """Background standard deviation of fitted Gaussian.
+
+        This is the standard error reported by :py:obj:`lmfit.model.Model.fit`.
+        """
+
+        self.r_squared, self._r_squared_setter = soft_signal_r_and_setter(float, 0.0)
+        """R-squared (goodness of fit) parameter reported by :py:obj:`lmfit.model.Model.fit`."""
+
+        super().__init__()
+        self._detectors = detectors
+        self._angle_map = angle_map
+        self._flood = flood if flood is not None else sc.scalar(value=1.0, dtype="float64")
+        self._fit_method = Gaussian()
+
+    def additional_readable_signals(self, dae: Dae) -> list[Device]:
+        """Expose fit parameters as readable signals.
+
+        :meta private:
+        """
+        return [
+            self.amp,
+            self.amp_err,
+            self.sigma,
+            self.sigma_err,
+            self.x0,
+            self.x0_err,
+            self.background,
+            self.background_err,
+        ]
+
+    async def reduce_data(self, dae: Dae) -> None:
+        """Perform the 'reduction'.
+
+        :meta private:
+        """
+        data = await dae.trigger_and_get_specdata()
+
+        # Filter to relevant detectors
+        data = data[self._detectors]
+
+        # Sum in ToF
+        data = data.sum(axis=1)
+
+        data = sc.array(dims=["spectrum"], values=data, variances=data + 0.5)
+        data /= self._flood
+
+        fit_method = Gaussian()
+
+        # Generate initial guesses and fit
+        guess = fit_method.guess()(self._angle_map, data.values)
+        result = fit_method.model().fit(
+            data.values,
+            x=self._angle_map,
+            **guess,  # pyright: ignore
+            weights=1 / (data.variances**0.5),
+        )
+
+        self._amp_setter(result.params["amp"].value)
+        self._amp_err_setter(result.params["amp"].stderr)
+        self._x0_setter(result.params["x0"].value)
+        self._x0_err_setter(result.params["x0"].stderr)
+        self._sigma_setter(result.params["sigma"].value)
+        self._sigma_err_setter(result.params["sigma"].stderr)
+        self._background_setter(result.params["background"].value)
+        self._background_err_setter(result.params["background"].stderr)
+
+        self._r_squared_setter(result.rsquared)

tests/devices/test_reflectometry.py

@@ -1,14 +1,17 @@
 # pyright: reportMissingParameterType=false
 import asyncio
-from unittest.mock import patch
+from unittest.mock import AsyncMock, MagicMock, patch
 
 import bluesky.plan_stubs as bps
+import numpy as np
 import pytest
 from ophyd_async.plan_stubs import ensure_connected
 from ophyd_async.testing import callback_on_mock_put, get_mock_put, set_mock_value
 
 from ibex_bluesky_core.devices import NoYesChoice
+from ibex_bluesky_core.devices.dae import Dae
 from ibex_bluesky_core.devices.reflectometry import (
+    AngleMappingReducer,
     ReflParameter,
     ReflParameterRedefine,
     refl_parameter,
@@ -59,3 +62,46 @@ async def test_fails_to_redefine_and_raises_if_not_in_manager_mode(RE):
         r" as not in manager mode.",
     ):
         await param.set(new_value)
+
+
+async def test_angle_mapping_reducer(dae):
+    reducer = AngleMappingReducer(
+        detectors=np.array([0, 1, 2, 3, 4], dtype=np.int32),
+        angle_map=np.array([10, 11, 12, 13, 14], dtype=np.float64),
+    )
+
+    await reducer.connect(mock=True)
+
+    fakedae = Dae(prefix="unittest:")
+    await fakedae.connect(mock=True)
+    fakedae.trigger_and_get_specdata = AsyncMock(
+        return_value=np.array([[0, 0], [0, 1], [0, 10], [0, 1], [0, 0]], dtype=np.float64)
+    )
+
+    await reducer.reduce_data(fakedae)
+
+    # Approximate fit parameters for above data
+    assert await reducer.x0.get_value() == pytest.approx(12.0)
+    assert await reducer.amp.get_value() == pytest.approx(10.0, abs=0.01)
+    assert await reducer.background.get_value() == pytest.approx(0.0, abs=0.01)
+
+    # Should get a 'good' fit to the above data
+    assert await reducer.r_squared.get_value() == pytest.approx(1.0, abs=0.01)
+
+
+def test_angle_mapping_reducer_interesting_signals():
+    reducer = AngleMappingReducer(
+        detectors=np.array([], dtype=np.int32),
+        angle_map=np.array([], dtype=np.float64),
+    )
+
+    assert reducer.additional_readable_signals(MagicMock()) == [
+        reducer.amp,
+        reducer.amp_err,
+        reducer.sigma,
+        reducer.sigma_err,
+        reducer.x0,
+        reducer.x0_err,
+        reducer.background,
+        reducer.background_err,
+    ]