fix: reduction workflow for 3D BEER detectors (#254)

* fix: make sure coord remains after binnning

* test: skip if cif-db not reachable

* fix: make sure coord remains after binnning

* test: skip if cif-db not reachable

* feat: round fitted t0 to nearest chopper opening

* remove comment

* refactor: better name for frame cutoff time coordinate

* test: add test reducing 3D dataset using automatic peak finder workflow

* refactor: better name for chopper delay

* docs: explain rounding step

* docs: motivation for clamping origin times

* docs: reformulate

Author: jokasimr

src/ess/beer/clustering.py

@@ -20,7 +20,7 @@ def cluster_events_by_streak(
 
     # We need to keep these coordinates after binning,
     # adding them to the binned data coords achieves this.
-    for coord in ('two_theta', 'Ltotal'):
+    for coord in ('two_theta', 'Ltotal', 'frame_cutoff_time'):
         da.bins.coords[coord] = sc.bins_like(da, da.coords[coord])
 
     h = da.bins.concat().hist(coarse_d=1000)

src/ess/beer/conversions.py

@@ -18,6 +18,7 @@
 
 def compute_tof_in_each_cluster(
     da: StreakClusteredData[RunType],
+    chopper_delay: WavelengthDefinitionChopperDelay,
     mod_period: ModulationPeriod,
 ) -> TofDetector[RunType]:
     """Fits a line through each cluster, the intercept of the line is t0.
@@ -32,6 +33,7 @@ def compute_tof_in_each_cluster(
        of the points in the cluster, and probably should belong to another cluster or
        are part of the background.
     3. Go back to 1) and iterate until convergence. A few iterations should be enough.
+    4. Finally, round the estimated t0 to the closest known chopper opening time.
     """
     if isinstance(da, sc.DataGroup):
         return sc.DataGroup(
@@ -42,7 +44,7 @@ def compute_tof_in_each_cluster(
     sin_theta_L = sc.sin(da.bins.coords['two_theta'] / 2) * da.bins.coords['Ltotal']
     t = time_of_arrival(
         da.bins.coords['event_time_offset'],
-        da.coords['tc'].to(unit=da.bins.coords['event_time_offset'].unit),
+        da.bins.coords['frame_cutoff_time'],
     )
     for _ in range(15):
         s, t0 = _linear_regression_by_bin(sin_theta_L, t, da.data)
@@ -54,11 +56,33 @@ def compute_tof_in_each_cluster(
             too_far_from_center=(distance_to_self > max_distance_from_streak_line),
         )
 
-    da = da.assign_coords(t0=sc.values(t0))
-    da = da.bins.assign_coords(tof=(t - sc.values(t0)))
+    # The t0 estimate from fitting is influenced by peak overlap, background,
+    # and other factors that can make the estimate offset from the true
+    # chopper opening time that it should match.
+    # We know the true chopper opening times, so instead of using the t0 estimte
+    # directly we can round the estimate to the closest chopper opening time.
+    # That way the t0 estimate becomes more robust and is guaranteed to correspond to
+    # a true chopper opening time.
+    t0 = _round_t0_to_nearest_chopper_opening(sc.values(t0), mod_period, chopper_delay)
+    da = da.assign_coords(t0=t0)
+    da = da.bins.assign_coords(tof=(t - t0))
     return da
 
 
+def _round_t0_to_nearest_chopper_opening(
+    t0: sc.Variable,
+    mod_period: sc.Variable,
+    chopper_delay: sc.Variable,
+) -> sc.Variable:
+    out = t0 - chopper_delay
+    out /= mod_period
+    out += 0.5
+    sc.floor(out, out=out)
+    out *= mod_period
+    out += chopper_delay
+    return out
+
+
 def _linear_regression_by_bin(
     x: sc.Variable, y: sc.Variable, w: sc.Variable
 ) -> tuple[sc.Variable, sc.Variable]:
@@ -118,14 +142,14 @@ def _compute_d_given_list_of_peaks(
 
 def time_of_arrival(
     event_time_offset: sc.Variable,
-    tc: sc.Variable,
+    frame_cutoff_time: sc.Variable,
 ):
     """Does frame unwrapping for pulse shaping chopper modes.
 
-    Events before the "cutoff time" `tc` are assumed to come from the previous pulse."""
+    Events before the "cutoff time" are assumed to come from the previous pulse."""
     _eto = event_time_offset
     T = sc.scalar(1 / 14, unit='s').to(unit=_eto.unit)
-    tc = tc.to(unit=_eto.unit)
+    tc = frame_cutoff_time.to(unit=_eto.unit)
     return sc.where(_eto >= tc, _eto, _eto + T)
 
 
@@ -135,13 +159,13 @@ def _tof_from_dhkl(
     coarse_dhkl: sc.Variable,
     Ltotal: sc.Variable,
     mod_period: sc.Variable,
-    time0: sc.Variable,
+    chopper_delay: sc.Variable,
 ) -> sc.Variable:
     """Computes tof for BEER given the dhkl peak that the event belongs to"""
     # Source: https://www.mcstas.org/download/components/current/contrib/NPI_tof_dhkl_detector.comp
     # tref = 2 * d_hkl * sin(theta) / hm * Ltotal
-    # tc = time_of_arrival - time0 - tref
-    # dt = floor(tc / mod_period + 0.5) * mod_period + time0
+    # tc = time_of_arrival - chopper_delay - tref
+    # dt = floor(tc / mod_period + 0.5) * mod_period + chopper_delay
     # tof = time_of_arrival - dt
     c = (-2 * 1.0 / (scipp.constants.h / scipp.constants.m_n)).to(
         unit=f'{time_of_arrival.unit}/m/angstrom'
@@ -150,12 +174,12 @@ def _tof_from_dhkl(
     out *= sc.sin(theta)
     out *= Ltotal
     out += time_of_arrival
-    out -= time0
+    out -= chopper_delay
     out /= mod_period
     out += 0.5
     sc.floor(out, out=out)
     out *= mod_period
-    out += time0
+    out += chopper_delay
     out *= -1
     out += time_of_arrival
     return out
@@ -168,7 +192,7 @@ def geometry_graph() -> GeometryCoordTransformGraph:
 def tof_from_known_dhkl_graph(
     mod_period: ModulationPeriod,
     pulse_length: PulseLength,
-    time0: WavelengthDefinitionChopperDelay,
+    chopper_delay: WavelengthDefinitionChopperDelay,
     dhkl_list: DHKLList,
     gg: GeometryCoordTransformGraph,
 ) -> TofCoordTransformGraph:
@@ -200,7 +224,7 @@ def _compute_coarse_dspacing(
         **graph.tof.elastic("tof"),
         'pulse_length': lambda: pulse_length,
         'mod_period': lambda: mod_period,
-        'time0': lambda: time0,
+        'chopper_delay': lambda: chopper_delay,
         'tof': _tof_from_dhkl,
         'time_of_arrival': time_of_arrival,
         'coarse_dhkl': _compute_coarse_dspacing,

src/ess/beer/data.py

@@ -44,6 +44,7 @@
         # Simulation with 3D detector model - almost no events
         # - only used to verify we can load the 3D geometry.
         "few_neutrons_3d_detector_example.h5": "md5:88cbe29cb539c8acebf9fd7cee9d3c57",
+        "silicon-mode9-3d-more-neutrons.h5": "md5:a6afdc4ee9827a57f88c2a3c6ef27383",
     },
 )
 
@@ -81,6 +82,10 @@ def mcstas_few_neutrons_3d_detector_example():
     return _registry.get_path('few_neutrons_3d_detector_example.h5')
 
 
+def mcstas_more_neutrons_3d_detector_example():
+    return _registry.get_path('silicon-mode9-3d-more-neutrons.h5')
+
+
 def duplex_peaks() -> Path:
     return _registry.get_path('duplex-dhkl.tab')
 

src/ess/beer/io.py

@@ -234,8 +234,8 @@ def _load_beer_mcstas(f, north_or_south=None, *, number):
     # Bin detector panel into rectangular "pixels"
     # similar in size to the physical detector pixels.
     da = da.bin(
-        y=sc.linspace('y', -0.5, 0.5, 501, unit='m'),
-        x=sc.linspace('x', -0.5, 0.5, 201, unit='m'),
+        y=sc.linspace('y', -0.5, 0.5, 201, unit='m'),
+        x=sc.linspace('x', -0.5, 0.5, 501, unit='m'),
     )
 
     # Compute the position of each pixel in the global coordinate system.
@@ -279,7 +279,7 @@ def _load_beer_mcstas(f, north_or_south=None, *, number):
     # Used in modulation mode automatic-peak-finding reduction to estimate d.
     # In practice this will probably be replaced by the regular tof workflow.
     # But I'm not 100% sure.
-    da.coords["tc"] = (
+    da.coords["frame_cutoff_time"] = (
         sc.constants.m_n
         / sc.constants.h
         * da.coords['wavelength_estimate']

tests/beer/mcstas_reduction_test.py

@@ -1,4 +1,5 @@
 import numpy as np
+import pytest
 import scipp as sc
 import scippneutron as scn
 from scipp.testing import assert_allclose
@@ -12,9 +13,14 @@
     duplex_peaks_array,
     mcstas_duplex,
     mcstas_few_neutrons_3d_detector_example,
+    mcstas_more_neutrons_3d_detector_example,
     mcstas_silicon_new_model,
 )
-from ess.beer.io import load_beer_mcstas, load_beer_mcstas_monitor
+from ess.beer.io import (
+    load_beer_mcstas,
+    load_beer_mcstas_monitor,
+    mcstas_chopper_delay_from_mode_new_simulations,
+)
 from ess.beer.types import DetectorBank, DHKLList, TofDetector
 from ess.reduce.nexus.types import Filename, SampleRun
 
@@ -40,10 +46,14 @@ def test_can_reduce_using_known_peaks_workflow():
     )
 
 
-def test_can_reduce_using_unknown_peaks_workflow():
+@pytest.mark.parametrize(
+    'fname', [mcstas_silicon_new_model(7), mcstas_more_neutrons_3d_detector_example()]
+)
+def test_can_reduce_using_unknown_peaks_workflow(fname):
     wf = BeerModMcStasWorkflow()
-    wf[Filename[SampleRun]] = mcstas_duplex(7)
+    wf[Filename[SampleRun]] = fname
     wf[DetectorBank] = DetectorBank.north
+    wf.insert(mcstas_chopper_delay_from_mode_new_simulations)
     da = wf.compute(TofDetector[SampleRun])
     da = da.transform_coords(
         ('dspacing',),
@@ -53,7 +63,11 @@ def test_can_reduce_using_unknown_peaks_workflow():
     max_peak_d = sc.midpoints(h['dspacing', np.argmax(h.values)].coords['dspacing'])[0]
     assert_allclose(
         max_peak_d,
-        sc.scalar(2.0407, unit='angstrom'),
+        # The two peaks around 1.6 are very similar in magnitude,
+        # so either of them can be bigger and that is fine.
+        sc.scalar(1.5677, unit='angstrom')
+        if max_peak_d < sc.scalar(1.6, unit='angstrom')
+        else sc.scalar(1.6374, unit='angstrom'),
         atol=sc.scalar(1e-2, unit='angstrom'),
     )
 

tests/diffraction/test_peaks.py

@@ -1,8 +1,23 @@
+import urllib.request
+
+import pytest
 import scipp as sc
 
 from ess.diffraction.peaks import dspacing_peaks_from_cif
 
 
+def _cifdb_is_reachable(timeout=5):
+    try:
+        req = urllib.request.Request('http://cod.ibt.lt', method="HEAD")
+        with urllib.request.urlopen(req, timeout=timeout):  # noqa: S310
+            return True
+    except Exception:
+        return False
+
+
+@pytest.mark.skipif(
+    not _cifdb_is_reachable(), reason='The CIF database can not be reached.'
+)
 def test_retreived_peak_list_has_expected_units():
     d = dspacing_peaks_from_cif(
         'codid::9011998',
@@ -12,6 +27,9 @@ def test_retreived_peak_list_has_expected_units():
     assert d.unit == 'barn'
 
 
+@pytest.mark.skipif(
+    not _cifdb_is_reachable(), reason='The CIF database can not be reached.'
+)
 def test_intensity_threshold_is_taken_into_account():
     d = dspacing_peaks_from_cif(
         'codid::9011998',
@@ -26,6 +44,9 @@ def test_intensity_threshold_is_taken_into_account():
     assert len(d) == 0
 
 
+@pytest.mark.skipif(
+    not _cifdb_is_reachable(), reason='The CIF database can not be reached.'
+)
 def test_retreived_peak_list_with_temperature_kwarg():
     d = dspacing_peaks_from_cif(
         'codid::9011998', sc.scalar(50, unit='barn'), uiso_temperature=300