Load time-dependent transformations from NeXus

packages/essreduce/pyproject.toml

@@ -31,7 +31,7 @@ dynamic = ["version"]
 
 dependencies = [
     "sciline>=25.11.0",
-    "scipp>=26.3.0",
+    "scipp>=26.3.1",
     "scippneutron>=25.11.1",
     "scippnexus>=25.06.0",
 ]

packages/essreduce/src/ess/reduce/nexus/__init__.py

@@ -15,6 +15,7 @@
 from . import types
 from ._nexus_loader import (
     compute_component_position,
+    compute_detector_position,
     extract_signal_data_array,
     group_event_data,
     load_all_components,
@@ -29,6 +30,7 @@
 __all__ = [
     'GenericNeXusWorkflow',
     'compute_component_position',
+    'compute_detector_position',
     'extract_signal_data_array',
     'group_event_data',
     'load_all_components',

packages/essreduce/src/ess/reduce/nexus/_nexus_loader.py

@@ -23,6 +23,8 @@
     NeXusFile,
     NeXusGroup,
     NeXusLocationSpec,
+    NeXusTransformation,
+    RunType,
 )
 
 
@@ -454,6 +456,55 @@ def _to_snx_selection(selection, *, for_events: bool) -> snx.typing.ScippIndex:
     return selection
 
 
+def compute_detector_position(
+    da: sc.DataArray,
+    *,
+    transform: NeXusTransformation[snx.NXdetector, RunType],
+    # Strictly speaking we could apply an offset by modifying the transformation chain,
+    # using a more generic implementation. However, this may in general require
+    # extending the chain and it is currently not clear if that is desirable. As far as
+    # I am aware the offset is currently mainly used for handling files from other
+    # facilities and it is not clear if it is needed for ESS data and should be kept at
+    # all.
+    offset: sc.Variable,
+) -> sc.Variable | sc.DataArray:
+    """Compute the positions of detector pixels.
+
+    Parameters
+    ----------
+    da:
+        Detector (event) data as returned by :func:`extract_signal_data_array`.
+    transform:
+        Transformation matrix for the detector.
+    offset:
+        Offset to add to the detector position.
+
+    Returns
+    -------
+    :
+        The detector position as a data array if ``transform`` is time-dependent
+        or as a variable otherwise.
+    """
+    # Note: We apply offset as early as possible, i.e., right in this function
+    # the detector array from the raw loader NeXus group, to prevent a source of bugs.
+    # If the NXdetector in the file is not 1-D, we want to match the order of dims.
+    # zip_pixel_offsets otherwise yields a vector with dimensions in the order given
+    # by the x/y/z offsets.
+    offsets = snx.zip_pixel_offsets(da.coords)
+    # Get the dims in the order of the detector data array, but filter out dims that
+    # don't exist in the offsets (e.g. the detector data may have a 'time' dimension).
+    dims = [dim for dim in da.dims if dim in offsets.dims]
+    offsets = offsets.transpose(dims).copy()
+    # We use the unit of the offsets as this is likely what the user expects.
+    if transform.value.unit is not None and transform.value.unit != '':
+        transform_value = transform.value.to(unit=offsets.unit)
+    else:
+        transform_value = transform.value
+    position = transform_value * offsets
+    position += offset.to(unit=position.unit, copy=False)
+    return position
+
+
 def load_data(
     file_path: FilePath | NeXusFile | NeXusGroup,
     selection: snx.typing.ScippIndex | slice = (),

packages/essreduce/src/ess/reduce/nexus/types.py

@@ -186,7 +186,61 @@ class NeXusData(sciline.Scope[Component, RunType, sc.DataArray], sc.DataArray):
 
 
 class Position(sciline.Scope[Component, RunType, sc.Variable], sc.Variable):
-    """Position of a component such as source, sample, monitor, or detector."""
+    """Position of a component that does not move, such as source or sample."""
+
+
+@dataclass(init=False, repr=False, slots=True)
+class DynamicPosition(Generic[Component, RunType]):
+    """Position of a potentially moving component such as an analyzer or detector.
+
+    The position can depend on time. In this case, a time coordinate is also stored.
+    Use ``position`` to get the position if it is scalar, or ``positions``
+    to get the position as a (potentially time-dependent) DataArray.
+    """
+
+    _position: sc.Variable
+    _time: sc.Variable | None
+
+    def __init__(self, pos: sc.DataArray | sc.Variable) -> None:
+        if pos.ndim == 0:
+            self._position = pos.data if isinstance(pos, sc.DataArray) else pos
+            self._time = None
+        else:
+            if not isinstance(pos, sc.DataArray):
+                raise sc.DimensionError(
+                    "Position is not a scalar, so it must be a DataArray"
+                )
+            self._position = pos.data
+            self._time = pos.coords['time']
+
+    @property
+    def is_dynamic(self) -> bool:
+        return self._time is not None
+
+    @property
+    def position(self) -> sc.Variable:
+        if self.is_dynamic:
+            raise sc.DimensionError(
+                "Position is time-dependent, use `positions` instead."
+            )
+        return self._position
+
+    @property
+    def positions(self) -> sc.DataArray:
+        da = sc.DataArray(self._position)
+        if self._time is not None:
+            da.coords['time'] = self._time
+        return da
+
+    def __str__(self) -> str:
+        if self.is_dynamic:
+            time_str = f", time={self._time}"
+        else:
+            time_str = ""
+        return f"Position(position={self._position}{time_str})"
+
+    def __repr__(self) -> str:
+        return f"Position(position={self._position}, time={self._time})"
 
 
 class DetectorPositionOffset(sciline.Scope[RunType, sc.Variable], sc.Variable):
@@ -225,6 +279,15 @@ class TimeInterval(Generic[RunType]):
     value: slice
 
 
+class TransformationTimeFilter(Generic[Component, RunType]):
+    """Filter for time-dependent transformations."""
+
+    def __new__(cls, x: Any) -> Any:
+        return x
+
+    def __call__(self, transform: sc.DataArray) -> sc.Variable | sc.DataArray: ...
+
+
 @dataclass
 class NeXusFileSpec(Generic[RunType]):
     value: FilePath | NeXusFile | NeXusGroup
@@ -280,25 +343,21 @@ class NeXusTransformationChain(
 
 @dataclass
 class NeXusTransformation(Generic[Component, RunType]):
-    value: sc.Variable
+    """A NeXus transformation computed from a transformation chain.
+
+    If the transformation is time-dependent, it is stored as a data array
+    with a 'time' coordinate.
+    Otherwise, the transformation is stored as a variable.
+    """
+
+    value: sc.Variable | sc.DataArray
 
     @staticmethod
     def from_chain(
         chain: NeXusTransformationChain[Component, RunType],
     ) -> 'NeXusTransformation[Component, RunType]':
-        """
-        Convert a transformation chain to a single transformation.
-
-        As transformation chains may be time-dependent, this method will need to select
-        a specific time point to convert to a single transformation. This may include
-        averaging as well as threshold checks. This is not implemented yet and we
-        therefore currently raise an error if the transformation chain does not compute
-        to a scalar.
-        """
-        if chain.transformations.sizes != {}:
-            raise ValueError(f"Expected scalar transformation, got {chain}")
-        transform = chain.compute()
-        return NeXusTransformation(value=transform)
+        """Convert a transformation chain to a single transformation."""
+        return NeXusTransformation(value=chain.compute())
 
 
 class RawChoppers(