Merge pull request #244 from scipp/histogram-monitors

Histogram monitors

Author: jl-wynen

src/ess/reduce/data.py

@@ -40,6 +40,15 @@ def get_path(self, name: str, unzip: bool = False) -> str:
 )
 
 
+_dream_registry = Registry(
+    instrument='dream',
+    files={
+        "TEST_977695_00068064.hdf": "md5:9e6ee9ec70d7c5e8c0c93b9e07e8949f",
+    },
+    version='2',
+)
+
+
 _loki_registry = Registry(
     instrument='loki',
     files={
@@ -94,3 +103,11 @@ def loki_tutorial_background_run_60393() -> str:
 def loki_tutorial_sample_transmission_run() -> str:
     """Sample transmission run (sample + sample holder/can + transmission monitor)."""
     return _loki_registry.get_path('60394-2022-02-28_2215.nxs')
+
+
+def dream_coda_test_file() -> str:
+    """CODA file for DREAM where most pulses have been removed.
+
+    See ``tools/shrink_nexus.py``.
+    """
+    return _dream_registry.get_path('TEST_977695_00068064.hdf')

src/ess/reduce/nexus/workflow.py

@@ -560,10 +560,10 @@ def __init__(
 class _DummyField:
     """Dummy field that can replace snx.Field in NXmonitor."""
 
-    def __init__(self):
+    def __init__(self, dim: str):
         self.attrs = {}
-        self.sizes = {'event_time_zero': 0}
-        self.dims = ('event_time_zero',)
+        self.sizes = {dim: 0}
+        self.dims = (dim,)
         self.shape = (0,)
 
     def __getitem__(self, key: Any) -> sc.Variable:
@@ -573,26 +573,39 @@ def __getitem__(self, key: Any) -> sc.Variable:
 class _StrippedMonitor(snx.NXmonitor):
     """Monitor definition without event data for ScippNexus.
 
-    Drops NXevent_data group, data is replaced by a dummy field.
+    Drops NXevent_data and NXdata groups, data is replaced by a dummy field.
     """
 
     def __init__(
         self, attrs: dict[str, Any], children: dict[str, snx.Field | snx.Group]
     ):
-        children = _drop(children, (snx.NXevent_data,))
-        children['data'] = _DummyField()
+        is_dense = snx.NXdata in (
+            getattr(child, 'nx_class', None) for child in children
+        )
+        children = _drop(children, (snx.NXevent_data, snx.NXdata))
+        children['data'] = _DummyField(dim='time' if is_dense else 'event_time_zero')
         super().__init__(attrs=attrs, children=children)
 
 
 def _add_variances(da: sc.DataArray) -> sc.DataArray:
     out = da.copy(deep=False)
     if out.bins is not None:
         content = out.bins.constituents['data']
-        if content.variances is None:
-            content.variances = content.values
+        content.data = _assign_values_as_variances(content.data)
+    elif out.variances is None:
+        out.data = _assign_values_as_variances(out.data)
     return out
 
 
+def _assign_values_as_variances(var: sc.Variable) -> sc.Variable:
+    try:
+        var.variances = var.values
+    except sc.VariancesError:
+        var = var.to(dtype=sc.DType.float64)
+        var.variances = var.values
+    return var
+
+
 def load_beamline_metadata_from_nexus(file_spec: NeXusFileSpec[SampleRun]) -> Beamline:
     """Load beamline metadata from a sample NeXus file."""
     return nexus.load_metadata(file_spec.value, Beamline)

tests/nexus/workflow_test.py

@@ -471,6 +471,21 @@ def monitor_event_data() -> workflow.NeXusData[FrameMonitor1, SampleRun]:
     )
 
 
+@pytest.fixture
+def monitor_histogram_data() -> workflow.NeXusData[FrameMonitor1, SampleRun]:
+    time = sc.epoch(unit='ns') + sc.arange('time', 1, 6, unit='s').to(unit='ns')
+    frame_time = sc.arange('frame_time', 12, unit='ms').to(unit='ns')
+    return workflow.NeXusData[FrameMonitor1, SampleRun](
+        sc.DataArray(
+            10.0
+            * sc.arange('x', 5 * 12, unit='counts').fold(
+                'x', sizes={'time': 5, 'frame_time': 12}
+            ),
+            coords={'time': time, 'frame_time': frame_time},
+        )
+    )
+
+
 def test_assemble_monitor_data_adds_events_as_values_and_coords(
     calibrated_monitor, monitor_event_data
 ) -> None:
@@ -482,7 +497,7 @@ def test_assemble_monitor_data_adds_events_as_values_and_coords(
     )
 
 
-def test_assemble_monitor_data_adds_variances_to_weights(
+def test_assemble_monitor_data_adds_variances_to_events(
     calibrated_monitor, monitor_event_data
 ) -> None:
     monitor_data = workflow.assemble_monitor_data(
@@ -494,6 +509,30 @@ def test_assemble_monitor_data_adds_variances_to_weights(
     )
 
 
+def test_assemble_monitor_data_adds_histogram_as_values_and_coords(
+    calibrated_monitor, monitor_histogram_data
+) -> None:
+    monitor_data = workflow.assemble_monitor_data(
+        calibrated_monitor, monitor_histogram_data
+    )
+    assert_identical(
+        monitor_data.drop_coords(tuple(calibrated_monitor.coords)),
+        monitor_histogram_data,
+    )
+
+
+def test_assemble_monitor_data_adds_variances_to_weights(
+    calibrated_monitor, monitor_histogram_data
+) -> None:
+    monitor_data = workflow.assemble_monitor_data(
+        calibrated_monitor, monitor_histogram_data
+    )
+    assert_identical(
+        sc.variances(monitor_data.drop_coords(tuple(calibrated_monitor.coords))),
+        monitor_histogram_data,
+    )
+
+
 def test_assemble_monitor_preserves_coords(calibrated_monitor, monitor_event_data):
     calibrated_monitor.coords['abc'] = sc.scalar(1.2)
     monitor_data = workflow.assemble_monitor_data(
@@ -510,7 +549,7 @@ def test_assemble_monitor_preserves_masks(calibrated_monitor, monitor_event_data
     assert 'mymask' in monitor_data.masks
 
 
-def test_load_monitor_workflow() -> None:
+def test_load_event_monitor_workflow() -> None:
     wf = LoadMonitorWorkflow(run_types=[SampleRun], monitor_types=[FrameMonitor1])
     wf[Filename[SampleRun]] = data.loki_tutorial_sample_run_60250()
     wf[NeXusName[FrameMonitor1]] = 'monitor_1'
@@ -519,6 +558,21 @@ def test_load_monitor_workflow() -> None:
     assert 'source_position' in da.coords
     assert da.bins is not None
     assert da.dims == ('event_time_zero',)
+    assert da.bins.constituents['data'].variances is not None
+
+
+def test_load_histogram_monitor_workflow() -> None:
+    wf = LoadMonitorWorkflow(run_types=[SampleRun], monitor_types=[FrameMonitor1])
+    wf[Filename[SampleRun]] = data.dream_coda_test_file()
+    wf[NeXusName[FrameMonitor1]] = 'monitor_bunker'
+    da = wf.compute(MonitorData[SampleRun, FrameMonitor1])
+    assert 'position' in da.coords
+    assert 'source_position' in da.coords
+    assert da.bins is None
+    assert set(da.dims) == {'time', 'frame_time'}
+    assert 'time' in da.coords.keys()
+    assert 'frame_time' in da.coords.keys()
+    assert da.variances is not None
 
 
 def test_load_detector_workflow() -> None:

tools/shrink_nexus.py

@@ -0,0 +1,83 @@
+# SPDX-License-Identifier: BSD-3-Clause
+# Copyright (c) 2025 Scipp contributors (https://github.com/scipp)
+"""Shrink a nexus file by removing all but the first N pulses."""
+
+import shutil
+import subprocess
+from pathlib import Path
+
+import h5py as h5
+import numpy as np
+import scippnexus as snx
+
+N_PULSE = 2
+
+DIR = Path("data/coda")
+IN_FNAME = "977695_00068064.hdf"
+TMP_FNAME = f"TMP_{IN_FNAME}"
+OUT_FNAME = f"TEST_{IN_FNAME}"
+
+
+def main() -> None:
+    shutil.copy(DIR / IN_FNAME, DIR / TMP_FNAME)
+    with snx.File(DIR / TMP_FNAME, 'r') as f:
+        entry = f['entry']
+        user_names = list(entry[snx.NXuser])
+        det_names = list(entry['instrument'][snx.NXdetector])
+        mon_names = list(entry['instrument'][snx.NXmonitor])
+
+    last_times = []
+
+    with h5.File(DIR / TMP_FNAME, 'a') as f:
+        entry = f['entry']
+        for name in user_names:
+            del entry[name]
+        del entry['scripts']
+
+        instr = entry['instrument']
+        for name in det_names:
+            det = instr[name]
+            del det['pixel_shape']
+
+            event_data_name = next(name for name in det if name.endswith('event_data'))
+            events = det[event_data_name]
+
+            last_times.append(
+                np.array(
+                    int(events['event_time_zero'][N_PULSE]),
+                    dtype=f'datetime64[{events["event_time_zero"].attrs["units"]}]',
+                )
+            )
+
+            n_events = events['event_index'][N_PULSE]
+            slice_dataset(events, 'event_index', N_PULSE)
+            slice_dataset(events, 'event_time_zero', N_PULSE)
+            slice_dataset(events, 'event_id', n_events)
+            slice_dataset(events, 'event_time_offset', n_events)
+
+        last_time = np.min(last_times)
+        for name in mon_names:
+            mon = instr[name]
+            data = mon[f'{name}_events']
+            time = data['time'][()].astype(f'datetime64[{data["time"].attrs["units"]}]')
+            if time.shape == (0,):
+                continue
+            n_times = np.argmax(time > last_time)
+            slice_dataset(data, 'time', n_times)
+
+            assert (  # noqa: S101
+                data.attrs['axes'][0] == 'time'
+            )  # required to slice as below:
+            slice_dataset(data, 'signal', n_times)
+
+    subprocess.check_call(['h5repack', DIR / TMP_FNAME, DIR / OUT_FNAME])  # noqa: S603, S607
+
+
+def slice_dataset(group, name, n):
+    attrs = group[name].attrs
+    group[name] = group.pop(name)[:n]
+    group[name].attrs.update(attrs)
+
+
+if __name__ == "__main__":
+    main()