Ticket42 sum detect monitor spectra

doc/devices/dae.md

@@ -205,6 +205,58 @@ Published signals:
 - `reducer.mon_counts_stddev` - uncertainty (standard deviation) of the summed monitor counts
 - `reducer.intensity_stddev` - uncertainty (standard deviation) of the normalised intensity
 
+### Time of Flight and Wavelength Bounding Spectra
+
+Scalar Normalizers (such as PeriodGoodFramesNormalizer, GoodFramesNormalizer) can be passed a
+summing function which can optionally sum spectra between provided time of flight or wavelength bounds.
+
+{py:obj}`ibex_bluesky_core.devices.simpledae.reducers.PeriodGoodFramesNormalizer`
+
+
+Here is an example showing creating a scalar normalizer with time of flight bounds from 15000 to 25000 μs, and summing 2 detectors:
+```
+import scipp
+
+bounds=scipp.array(dims=["tof"], values=[15000.0, 25000.0], unit=scipp.units.us)
+
+reducer = PeriodGoodFramesNormalizer(
+    prefix=get_pv_prefix(),
+    detector_spectra=[1, 2],
+    summer=tof_bounded_spectra(bounds)
+)
+```
+
+{py:obj}`ibex_bluesky_core.devices.simpledae.reducers.tof_bounded_spectra` 
+
+
+Monitor Normalizers, which have both a monitor as well as detector, can be passed a summing function for each of these components independently, e.g. the detector can use time of flight while the monitor uses wavelength. tof_bounded_spectra assumes that all pixels being summed share the same flight-path length. Where two separate instances of tof_bounded_spectra are used, such as in DetectorMonitorNormalizer, these may have different flight path lengths from each other.
+
+Here is an example with wavelength bounding used to sum the monitor component, and time of flight bounding for the detector summing spectra: 
+
+```
+import scipp
+
+wavelength_bounds = scipp.array(dims=["tof"], values=[0.0, 5.1], unit=scipp.units.angstrom, dtype="float64")
+total_flight_path_length = scipp.scalar(value=85.0, unit=sc.units.m),
+tof_bounds = scipp.array(dims=["tof"], values=[15000, 25000], unit=scipp.units.us)
+
+reducer = MonitorNormalizer(
+    prefix=get_pv_prefix(),
+    detector_spectra=[1],
+    monitor_spectra=[2],
+    detector_summer=wavelength_bounded_spectra(wavelength_bounds, total_flight_path_length),
+    monitor_summer=tof_bounded_spectra(tof_bounds)
+)
+```
+{py:obj}`ibex_bluesky_core.devices.simpledae.reducers.wavelength_bounded_spectra`
+
+
+- In either case, the bounds are passed as a scipp array, which needs a `dims` attribute, `values` passed
+as a list, and `units` (μs/microseconds for time of flight bounding, and angstrom for wavelength bounding)
+
+- If you don't specify either of these options, they will default to an summing over the entire spectrum.
+
+
 ## Waiters
 
 A `waiter` defines an arbitrary strategy for how long to count at each point.
@@ -362,3 +414,17 @@ A `DaeSpectrum` object provides 3 arrays:
 
 The `Dae` base class does not provide any spectra by default. User-level classes should specify 
 the set of spectra which they are interested in.
+
+
+
+
+Spectra can be summed between two bounds based on time of flight bounds, or wavelength bounds, for both detector and monitor normalizers.
+
+Both Scalar Normalizers (PeriodGoodFramesNormalizer, GoodFramesNormalizer) and MonitorNormalizers
+accept the following arguments: 
+- `detector_summer`: sums counts using pre-existing bounds, or sums using time of flight bounds,
+or wavelength bounds.
+- `monitor_summer` (MonitorNormalizer only): sums counts using pre-existing bounds,
+or sums using time of flight bounds, or wavelength bounds.
+
+For both options, the default, if none is specified, is to use pre-existing bounds.

pyproject.toml

@@ -49,6 +49,7 @@ dependencies = [
   "numpy",
   "scipp",
   "tzdata",
+  "scippneutron",
 ]
 
 [project.optional-dependencies]

src/ibex_bluesky_core/devices/dae/dae_spectra.py

@@ -125,6 +125,9 @@ async def read_spectrum_dataarray(self) -> sc.DataArray:
                 values=counts,
                 variances=counts + VARIANCE_ADDITION,
                 unit=sc.units.counts,
+                dtype="float64",
             ),
-            coords={"tof": sc.array(dims=["tof"], values=tof_edges, unit=sc.Unit(unit))},
+            coords={
+                "tof": sc.array(dims=["tof"], values=tof_edges, unit=sc.Unit(unit), dtype="float64")
+            },
         )

src/ibex_bluesky_core/devices/simpledae/reducers.py

@@ -4,8 +4,8 @@
 import logging
 import math
 from abc import ABC, abstractmethod
-from collections.abc import Collection, Sequence
-from typing import TYPE_CHECKING
+from collections.abc import Awaitable, Collection, Sequence
+from typing import TYPE_CHECKING, Callable
 
 import scipp as sc
 from ophyd_async.core import (
@@ -15,6 +15,7 @@
     StandardReadable,
     soft_signal_r_and_setter,
 )
+from scippneutron import conversion
 
 from ibex_bluesky_core.devices.dae.dae_spectra import DaeSpectra
 from ibex_bluesky_core.devices.simpledae.strategies import Reducer
@@ -33,6 +34,8 @@ async def sum_spectra(spectra: Collection[DaeSpectra]) -> sc.Variable | sc.DataA
 
     Returns a scipp scalar, which has .value and .variance properties for accessing the sum
     and variance respectively of the summed counts.
+
+    More info on scipp scalars can be found here: https://scipp.github.io/generated/functions/scipp.scalar.html
     """
     logger.info("Summing %d spectra using scipp", len(spectra))
     summed_counts = sc.scalar(value=0, unit=sc.units.counts, dtype="float64")
@@ -42,15 +45,99 @@ async def sum_spectra(spectra: Collection[DaeSpectra]) -> sc.Variable | sc.DataA
     return summed_counts
 
 
+def tof_bounded_spectra(
+    bounds: sc.Variable,
+) -> Callable[[Collection[DaeSpectra]], Awaitable[sc.Variable | sc.DataArray]]:
+    """Sum a set of neutron spectra between the specified time of flight bounds.
+
+    Args:
+        bounds: A scipp array of size 2, no variances, unit of us,
+            where the second element must be larger than the first.
+
+    Returns a scipp scalar, which has .value and .variance properties for accessing the sum
+    and variance respectively of the summed counts.
+
+    More info on scipp arrays and scalars can be found here: https://scipp.github.io/generated/functions/scipp.scalar.html
+
+    """
+    bounds_value = 2
+    if "tof" not in bounds.dims:
+        raise ValueError("Should contain tof dims")
+    if bounds.sizes["tof"] != bounds_value:
+        raise ValueError("Should contain lower and upper bound")
+
+    async def sum_spectra_with_tof(spectra: Collection[DaeSpectra]) -> sc.Variable | sc.DataArray:
+        """Sum spectra bounded by a time of flight upper and lower bound."""
+        summed_counts = sc.scalar(value=0, unit=sc.units.counts, dtype="float64")
+        for spec in asyncio.as_completed([s.read_spectrum_dataarray() for s in spectra]):
+            tof_bound_spectra = await spec
+            summed_counts += tof_bound_spectra.rebin({"tof": bounds}).sum()
+        return summed_counts
+
+    return sum_spectra_with_tof
+
+
+def wavelength_bounded_spectra(
+    bounds: sc.Variable, total_flight_path_length: sc.Variable
+) -> Callable[[Collection[DaeSpectra]], Awaitable[sc.Variable | sc.DataArray]]:
+    """Sum a set of neutron spectra between the specified wavelength bounds.
+
+    Args:
+        bounds: A scipp array of size 2 of wavelength bounds, in units of angstrom,
+            where the second element must be larger than the first.
+        total_flight_path_length: A scipp scalar of Ltotal (total flight path length), the path
+            length from neutron source to detector or monitor, in units of meters.
+
+    Time of flight is converted to wavelength using scipp neutron's library function
+        `wavelength_from_tof`, more info on which can be found here:
+        https://scipp.github.io/scippneutron/generated/modules/scippneutron.conversion.tof.wavelength_from_tof.html
+
+    Returns a scipp scalar, which has .value and .variance properties for accessing the sum
+    and variance respectively of the summed counts.
+
+    """
+    bounds_value = 2
+
+    if "tof" not in bounds.dims:
+        raise ValueError("Should contain tof dims")
+    if bounds.sizes["tof"] != bounds_value:
+        raise ValueError("Should contain lower and upper bound")
+
+    async def sum_spectra_with_wavelength(
+        spectra: Collection[DaeSpectra],
+    ) -> sc.Variable | sc.DataArray:
+        """Sum a set of spectra between the specified wavelength bounds."""
+        summed_counts = sc.scalar(value=0, unit=sc.units.counts, dtype="float64")
+        for spec in asyncio.as_completed([s.read_spectrum_dataarray() for s in spectra]):
+            wavelength_bounded_spectra = await spec
+            wavelength_coord = conversion.tof.wavelength_from_tof(
+                tof=wavelength_bounded_spectra.coords["tof"], Ltotal=total_flight_path_length
+            )
+            wavelength_bounded_spectra.coords["tof"] = wavelength_coord
+            summed_counts += wavelength_bounded_spectra.rebin({"tof": bounds}).sum()
+        return summed_counts
+
+    return sum_spectra_with_wavelength
+
+
 class ScalarNormalizer(Reducer, StandardReadable, ABC):
     """Sum a set of user-specified spectra, then normalize by a scalar signal."""
 
-    def __init__(self, prefix: str, detector_spectra: Sequence[int]) -> None:
+    def __init__(
+        self,
+        prefix: str,
+        detector_spectra: Sequence[int],
+        sum_detector: Callable[
+            [Collection[DaeSpectra]], Awaitable[sc.Variable | sc.DataArray]
+        ] = sum_spectra,
+    ) -> None:
         """Init.
 
         Args:
             prefix: the PV prefix of the instrument to get spectra from (e.g. IN:DEMO:)
             detector_spectra: a sequence of spectra numbers (detectors) to sum.
+            sum_detector: takes spectra objects, reads from them, and returns a scipp scalar
+                describing the detector intensity. Defaults to summing over the entire spectrum.
 
         """
         self.detectors = DeviceVector(
@@ -71,6 +158,7 @@ def __init__(self, prefix: str, detector_spectra: Sequence[int]) -> None:
         self.intensity_stddev, self._intensity_stddev_setter = soft_signal_r_and_setter(
             float, 0.0, precision=INTENSITY_PRECISION
         )
+        self.sum_detector = sum_detector
 
         super().__init__(name="")
 
@@ -82,7 +170,7 @@ async def reduce_data(self, dae: "SimpleDae") -> None:
         """Apply the normalization."""
         logger.info("starting reduction")
         summed_counts, denominator = await asyncio.gather(
-            sum_spectra(self.detectors.values()), self.denominator(dae).get_value()
+            self.sum_detector(self.detectors.values()), self.denominator(dae).get_value()
         )
 
         self._det_counts_setter(float(summed_counts.value))
@@ -132,14 +220,29 @@ class MonitorNormalizer(Reducer, StandardReadable):
     """Normalize a set of user-specified detector spectra by user-specified monitor spectra."""
 
     def __init__(
-        self, prefix: str, detector_spectra: Sequence[int], monitor_spectra: Sequence[int]
+        self,
+        prefix: str,
+        detector_spectra: Sequence[int],
+        monitor_spectra: Sequence[int],
+        sum_detector: Callable[
+            [Collection[DaeSpectra]], Awaitable[sc.Variable | sc.DataArray]
+        ] = sum_spectra,
+        sum_monitor: Callable[
+            [Collection[DaeSpectra]], Awaitable[sc.Variable | sc.DataArray]
+        ] = sum_spectra,
     ) -> None:
         """Init.
 
         Args:
             prefix: the PV prefix of the instrument to get spectra from (e.g. IN:DEMO:)
             detector_spectra: a sequence of spectra numbers (detectors) to sum.
             monitor_spectra: a sequence of spectra number (monitors) to sum and normalize by.
+            sum_detector: takes spectra objects, reads from them, and returns a scipp scalar
+                describing the detector intensity. Defaults to summing over the entire spectrum.
+            sum_monitor: takes spectra objects, reads from them, and returns a scipp scalar
+                describing the monitor intensity. Defaults to summing over the entire spectrum.
+
+        Scipp scalars are described in further detail here: https://scipp.github.io/generated/functions/scipp.scalar.html
 
         """
         dae_prefix = prefix + "DAE:"
@@ -165,14 +268,17 @@ def __init__(
         self.intensity_stddev, self._intensity_stddev_setter = soft_signal_r_and_setter(
             float, 0.0, precision=INTENSITY_PRECISION
         )
+        self.sum_detector = sum_detector
+        self.sum_monitor = sum_monitor
 
         super().__init__(name="")
 
     async def reduce_data(self, dae: "SimpleDae") -> None:
         """Apply the normalization."""
         logger.info("starting reduction")
         detector_counts, monitor_counts = await asyncio.gather(
-            sum_spectra(self.detectors.values()), sum_spectra(self.monitors.values())
+            self.sum_detector(self.detectors.values()),
+            self.sum_monitor(self.monitors.values()),
         )
 
         if monitor_counts.value == 0.0:  # To avoid zero division

src/ibex_bluesky_core/run_engine/__init__.py

@@ -72,6 +72,7 @@ def get_run_engine() -> RunEngine:
         RE.subscribe(lambda name, document: ...)
 
     For full documentation about the run engine, see:
+
     - https://nsls-ii.github.io/bluesky/tutorial.html#the-runengine
     - https://nsls-ii.github.io/bluesky/run_engine_api.html
     """