Merge remote-tracking branch 'origin/main' into beer-powdiff

Author: jokasimr

docs/user-guide/dream/index.md

@@ -54,5 +54,5 @@ dream-instrument-view
 workflow-widget-dream
 dream-detector-diagnostics
 dream-visualize-absorption
-dream-make-tof-lookup-table
+dream-make-wavelength-lookup-table
 ```

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_wavelength_in_each_cluster(
     da: StreakClusteredData[RunType],
+    chopper_delay: WavelengthDefinitionChopperDelay,
     mod_period: ModulationPeriod,
     graph: GeometryCoordTransformGraph,
 ) -> WavelengthDetector[RunType]:
@@ -47,7 +48,7 @@ def compute_wavelength_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)
@@ -59,12 +60,33 @@ def compute_wavelength_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)))
-    da = da.transform_coords(('wavelength',), graph=graph)
+    # 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]:
@@ -124,14 +146,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)
 
 
@@ -141,13 +163,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'
@@ -156,12 +178,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
@@ -205,7 +227,7 @@ def tof_from_known_dhkl_graph(
     da: RawDetector[RunType],
     mod_period: ModulationPeriod,
     pulse_length: PulseLength,
-    time0: WavelengthDefinitionChopperDelay,
+    chopper_delay: WavelengthDefinitionChopperDelay,
     dhkl_list: DHKLList,
     gg: GeometryCoordTransformGraph,
 ) -> ElasticCoordTransformGraph[RunType]:
@@ -236,7 +258,7 @@ def _compute_coarse_dspacing(
         **gg,
         '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

@@ -257,8 +257,8 @@ def _load_beer_mcstas(f, north_or_south=None, number=None):
     # 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.
@@ -302,7 +302,7 @@ def _load_beer_mcstas(f, north_or_south=None, number=None):
     # 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,10 +13,15 @@
     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.types import DetectorBank, DHKLList
+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
 from ess.reduce.unwrap.types import WavelengthDetector
 
@@ -41,11 +47,15 @@ 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
-    da = wf.compute(WavelengthDetector[SampleRun])
+    wf.insert(mcstas_chopper_delay_from_mode_new_simulations)
+    da = wf.compute(TofDetector[SampleRun])
     da = da.transform_coords(
         ('dspacing',),
         graph=scn.conversion.graph.tof.elastic('tof'),
@@ -54,7 +64,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