21 generate dae uncertainties

doc/architectural_decisions/002-use-ophyd-async.md

@@ -34,4 +34,4 @@ We will use `ophyd-async`.
 - `ophyd-async` will allow us to do fly scanning, if required in future, more easily than `ophyd`
 - Developers will need to understand more details about `asyncio`
 - Developers will have to write somewhat less boilerplate code in `ophyd-async` as compared to `ophyd`
-- Our devices should be more comparable to other facilities on site (e.g. DLS) who are using `ophyd-async`.
+- Our devices should be more comparable to other facilities on site (e.g. DLS) who are using `ophyd-async`.

doc/architectural_decisions/004-use-scipp.md

@@ -0,0 +1,23 @@
+# Use `scipp`
+
+## Status
+
+Current
+
+## Context
+
+A decision needs to be made about whether to use scipp, numpy, uncertainties or develop our own library for the purpose of providing support for generating uncertanties on our counts data.
+
+# Decision
+
+We will be using scipp.
+
+# Justification & Consequences
+
+- `scipp` is being developed at ESS with past input from STFC, so is well suited for neutron counts data.
+- `scipp` has a `numpy`-like interface but handles units and uncertainties by default under-the-hood.
+- Neither `numpy` or `uncertanties` have exactly the functionality we would need, so the solution using them would be a mix of the libraries and our own code, there would be more places to go wrong. Maintainability.
+- `uncertainties` package tracks correlations so may have bad scaling on "large" arrays like counts data from the DAE.
+- Developing our own uncertainties library will take time to understand and then implement. All of the functionality that we need has been done beforehand, so better to not waste time & effort.
+- Less expertise with this library on site (mitigation: don't do too much which is very complicated with it)
+- Potentially duplicates some of `mantid`'s functionality: (mitigation: use `scipp` for "simple" things, use `mantid` in future if people want to do "full" data reduction pipelines)

doc/devices/dae.md

@@ -177,6 +177,8 @@ Published signals:
 - `simpledae.good_frames` - the number of good frames reported by the DAE
 - `reducer.det_counts` - summed detector counts for all of the user-provided spectra
 - `reducer.intensity` - normalized intensity (`det_counts / good_frames`)
+- `reducer.det_counts_stddev` - uncertainty (standard deviation) of the summed detector counts
+- `reducer.intensity_stddev` - uncertainty (standard deviation) of the normalised intensity
 
 ### PeriodGoodFramesNormalizer
 
@@ -187,6 +189,8 @@ Published signals:
 - `simpledae.period.good_frames` - the number of good frames reported by the DAE
 - `reducer.det_counts` - summed detector counts for all of the user-provided spectra
 - `reducer.intensity` - normalized intensity (`det_counts / good_frames`)
+- `reducer.det_counts_stddev` - uncertainty (standard deviation) of the summed detector counts
+- `reducer.intensity_stddev` - uncertainty (standard deviation) of the normalised intensity
 
 ### DetectorMonitorNormalizer
 
@@ -197,6 +201,9 @@ Published signals:
 - `reducer.det_counts` - summed detector counts for the user-provided detector spectra
 - `reducer.mon_counts` - summed monitor counts for the user-provided monitor spectra
 - `reducer.intensity` - normalized intensity (`det_counts / mon_counts`)
+- `reducer.det_counts_stddev` - uncertainty (standard deviation) of the summed detector counts
+- `reducer.mon_counts_stddev` - uncertainty (standard deviation) of the summed monitor counts
+- `reducer.intensity_stddev` - uncertainty (standard deviation) of the normalised intensity
 
 ## Waiters
 

manual_system_tests/dae_scan.py

@@ -87,7 +87,9 @@ def dae_scan_plan() -> Generator[Msg, None, None]:
                     block.name,
                     controller.run_number.name,
                     reducer.intensity.name,
+                    reducer.intensity_stddev.name,
                     reducer.det_counts.name,
+                    reducer.det_counts_stddev.name,
                     dae.good_frames.name,
                 ]
             ),

src/ibex_bluesky_core/devices/simpledae/reducers.py

@@ -1,6 +1,7 @@
 """DAE data reduction strategies."""
 
 import asyncio
+import math
 from abc import ABCMeta, abstractmethod
 from typing import TYPE_CHECKING, Collection, Sequence
 
@@ -20,6 +21,9 @@
     from ibex_bluesky_core.devices.simpledae import SimpleDae
 
 
+INTENSITY_PRECISION = 6
+
+
 async def sum_spectra(spectra: Collection[DaeSpectra]) -> sc.Variable | sc.DataArray:
     """Read and sum a number of spectra from the DAE.
 
@@ -51,7 +55,16 @@ def __init__(self, prefix: str, detector_spectra: Sequence[int]) -> None:
         )
 
         self.det_counts, self._det_counts_setter = soft_signal_r_and_setter(float, 0.0)
-        self.intensity, self._intensity_setter = soft_signal_r_and_setter(float, 0.0, precision=6)
+        self.intensity, self._intensity_setter = soft_signal_r_and_setter(
+            float, 0.0, precision=INTENSITY_PRECISION
+        )
+
+        self.det_counts_stddev, self._det_counts_stddev_setter = soft_signal_r_and_setter(
+            float, 0.0
+        )
+        self.intensity_stddev, self._intensity_stddev_setter = soft_signal_r_and_setter(
+            float, 0.0, precision=INTENSITY_PRECISION
+        )
 
         super().__init__(name="")
 
@@ -66,14 +79,28 @@ async def reduce_data(self, dae: "SimpleDae") -> None:
         )
 
         self._det_counts_setter(float(summed_counts.value))
-        self._intensity_setter(float(summed_counts.value) / denominator)
+
+        if denominator == 0.0:  # To avoid zero division
+            self._intensity_setter(0.0)
+            intensity_var = 0.0
+        else:
+            intensity = summed_counts / denominator
+            self._intensity_setter(intensity.value)
+            intensity_var = intensity.variance if intensity.variance is not None else 0.0
+
+        detector_counts_var = 0.0 if summed_counts.variance is None else summed_counts.variance
+
+        self._det_counts_stddev_setter(math.sqrt(detector_counts_var))
+        self._intensity_stddev_setter(math.sqrt(intensity_var))
 
     def additional_readable_signals(self, dae: "SimpleDae") -> list[Device]:
         """Publish interesting signals derived or used by this reducer."""
         return [
             self.det_counts,
             self.intensity,
             self.denominator(dae),
+            self.det_counts_stddev,
+            self.intensity_stddev,
         ]
 
 
@@ -117,7 +144,19 @@ def __init__(
 
         self.det_counts, self._det_counts_setter = soft_signal_r_and_setter(float, 0.0)
         self.mon_counts, self._mon_counts_setter = soft_signal_r_and_setter(float, 0.0)
-        self.intensity, self._intensity_setter = soft_signal_r_and_setter(float, 0.0, precision=6)
+        self.intensity, self._intensity_setter = soft_signal_r_and_setter(
+            float, 0.0, precision=INTENSITY_PRECISION
+        )
+
+        self.det_counts_stddev, self._det_counts_stddev_setter = soft_signal_r_and_setter(
+            float, 0.0
+        )
+        self.mon_counts_stddev, self._mon_counts_stddev_setter = soft_signal_r_and_setter(
+            float, 0.0
+        )
+        self.intensity_stddev, self._intensity_stddev_setter = soft_signal_r_and_setter(
+            float, 0.0, precision=INTENSITY_PRECISION
+        )
 
         super().__init__(name="")
 
@@ -127,14 +166,33 @@ async def reduce_data(self, dae: "SimpleDae") -> None:
             sum_spectra(self.detectors.values()), sum_spectra(self.monitors.values())
         )
 
+        if monitor_counts.value == 0.0:  # To avoid zero division
+            self._intensity_setter(0.0)
+            intensity_var = 0.0
+
+        else:
+            intensity = detector_counts / monitor_counts
+            self._intensity_setter(float(intensity.value))
+            intensity_var = intensity.variance if intensity.variance is not None else 0.0
+
+        self._intensity_stddev_setter(math.sqrt(intensity_var))
+
         self._det_counts_setter(float(detector_counts.value))
         self._mon_counts_setter(float(monitor_counts.value))
-        self._intensity_setter(float((detector_counts / monitor_counts).value))
+
+        detector_counts_var = 0.0 if detector_counts.variance is None else detector_counts.variance
+        monitor_counts_var = 0.0 if monitor_counts.variance is None else monitor_counts.variance
+
+        self._det_counts_stddev_setter(math.sqrt(detector_counts_var))
+        self._mon_counts_stddev_setter(math.sqrt(monitor_counts_var))
 
     def additional_readable_signals(self, dae: "SimpleDae") -> list[Device]:
         """Publish interesting signals derived or used by this reducer."""
         return [
             self.det_counts,
             self.mon_counts,
             self.intensity,
+            self.det_counts_stddev,
+            self.mon_counts_stddev,
+            self.intensity_stddev,
         ]