diff --git a/doc/plans/reflectometry/psi_scan_linear/psi_eq_0.png b/doc/plans/reflectometry/psi_scan_linear/psi_eq_0.png new file mode 100644 index 00000000..f7acc5fa Binary files /dev/null and b/doc/plans/reflectometry/psi_scan_linear/psi_eq_0.png differ diff --git a/doc/plans/reflectometry/psi_scan_linear/psi_neq_0.png b/doc/plans/reflectometry/psi_scan_linear/psi_neq_0.png new file mode 100644 index 00000000..489e6d1a Binary files /dev/null and b/doc/plans/reflectometry/psi_scan_linear/psi_neq_0.png differ diff --git a/doc/plans/reflectometry/psi_scan_linear_detector.md b/doc/plans/reflectometry/psi_scan_linear_detector.md new file mode 100644 index 00000000..c67d496d --- /dev/null +++ b/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 diff --git a/src/ibex_bluesky_core/devices/reflectometry.py b/src/ibex_bluesky_core/devices/reflectometry.py index 2cb43349..c22d68f2 100644 --- a/src/ibex_bluesky_core/devices/reflectometry.py +++ b/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 ` + 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) diff --git a/tests/devices/test_reflectometry.py b/tests/devices/test_reflectometry.py index 515725db..6204d0e7 100644 --- a/tests/devices/test_reflectometry.py +++ b/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, + ]