From 9d3af2199d9f7c6fd0cb8435108d3b6d9d5b1648 Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Tue, 3 Mar 2026 11:16:13 +0100 Subject: [PATCH 01/12] fix: make sure coord remains after binnning --- src/ess/beer/clustering.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/ess/beer/clustering.py b/src/ess/beer/clustering.py index e37395dc..405ef4e8 100644 --- a/src/ess/beer/clustering.py +++ b/src/ess/beer/clustering.py @@ -22,6 +22,9 @@ def cluster_events_by_streak( # adding them to the binned data coords achieves this. for coord in ('two_theta', 'Ltotal'): da.bins.coords[coord] = sc.bins_like(da, da.coords[coord]) + # Making them scalar also keeps them after binning, + # but without having to store them for each event. + da.coords['tc'] = da.coords['tc'].mean() h = da.bins.concat().hist(coarse_d=1000) i_peaks, _ = find_peaks( From e525a3ba86ce135c76881416b2ec1f6a549f1adf Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Tue, 3 Mar 2026 13:30:12 +0100 Subject: [PATCH 02/12] test: skip if cif-db not reachable --- tests/diffraction/test_peaks.py | 21 +++++++++++++++++++++ 1 file changed, 21 insertions(+) diff --git a/tests/diffraction/test_peaks.py b/tests/diffraction/test_peaks.py index c77b47a3..b58e9026 100644 --- a/tests/diffraction/test_peaks.py +++ b/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 From 561b787dd8d43d9aa64d6edbdf18a2e0b165c595 Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Tue, 3 Mar 2026 11:16:13 +0100 Subject: [PATCH 03/12] fix: make sure coord remains after binnning --- src/ess/beer/clustering.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/ess/beer/clustering.py b/src/ess/beer/clustering.py index e37395dc..405ef4e8 100644 --- a/src/ess/beer/clustering.py +++ b/src/ess/beer/clustering.py @@ -22,6 +22,9 @@ def cluster_events_by_streak( # adding them to the binned data coords achieves this. for coord in ('two_theta', 'Ltotal'): da.bins.coords[coord] = sc.bins_like(da, da.coords[coord]) + # Making them scalar also keeps them after binning, + # but without having to store them for each event. + da.coords['tc'] = da.coords['tc'].mean() h = da.bins.concat().hist(coarse_d=1000) i_peaks, _ = find_peaks( From 28545d62c8bb70341a96d20fdae81c89dbf117ca Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Tue, 3 Mar 2026 13:30:12 +0100 Subject: [PATCH 04/12] test: skip if cif-db not reachable --- tests/diffraction/test_peaks.py | 21 +++++++++++++++++++++ 1 file changed, 21 insertions(+) diff --git a/tests/diffraction/test_peaks.py b/tests/diffraction/test_peaks.py index c77b47a3..b58e9026 100644 --- a/tests/diffraction/test_peaks.py +++ b/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 From d9890bc7ec251ac939a8205490505011ebf9aef1 Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Tue, 24 Mar 2026 14:10:21 +0100 Subject: [PATCH 05/12] feat: round fitted t0 to nearest chopper opening --- src/ess/beer/clustering.py | 3 +-- src/ess/beer/conversions.py | 23 ++++++++++++++++++++--- 2 files changed, 21 insertions(+), 5 deletions(-) diff --git a/src/ess/beer/clustering.py b/src/ess/beer/clustering.py index 405ef4e8..c996c132 100644 --- a/src/ess/beer/clustering.py +++ b/src/ess/beer/clustering.py @@ -20,11 +20,10 @@ 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', 'tc'): da.bins.coords[coord] = sc.bins_like(da, da.coords[coord]) # Making them scalar also keeps them after binning, # but without having to store them for each event. - da.coords['tc'] = da.coords['tc'].mean() h = da.bins.concat().hist(coarse_d=1000) i_peaks, _ = find_peaks( diff --git a/src/ess/beer/conversions.py b/src/ess/beer/conversions.py index 152f559b..0820f58c 100644 --- a/src/ess/beer/conversions.py +++ b/src/ess/beer/conversions.py @@ -18,6 +18,7 @@ def compute_tof_in_each_cluster( da: StreakClusteredData[RunType], + time0: WavelengthDefinitionChopperDelay, mod_period: ModulationPeriod, ) -> TofDetector[RunType]: """Fits a line through each cluster, the intercept of the line is t0. @@ -42,7 +43,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['tc'], ) for _ in range(15): s, t0 = _linear_regression_by_bin(sin_theta_L, t, da.data) @@ -54,11 +55,27 @@ 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))) + t0 = _round_t0_to_nearest_chopper_opening(sc.values(t0), mod_period, time0) + 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, + time0: sc.Variable, +) -> sc.Variable: + out = t0.copy() + out -= time0 + out /= mod_period + out += 0.5 + sc.floor(out, out=out) + out *= mod_period + out += time0 + return out + + def _linear_regression_by_bin( x: sc.Variable, y: sc.Variable, w: sc.Variable ) -> tuple[sc.Variable, sc.Variable]: From cbe5572c3687103365c8c14e0dfb4ffbcd8b978a Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Tue, 24 Mar 2026 14:13:00 +0100 Subject: [PATCH 06/12] remove comment --- src/ess/beer/clustering.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/ess/beer/clustering.py b/src/ess/beer/clustering.py index c996c132..64abb6ab 100644 --- a/src/ess/beer/clustering.py +++ b/src/ess/beer/clustering.py @@ -22,8 +22,6 @@ def cluster_events_by_streak( # adding them to the binned data coords achieves this. for coord in ('two_theta', 'Ltotal', 'tc'): da.bins.coords[coord] = sc.bins_like(da, da.coords[coord]) - # Making them scalar also keeps them after binning, - # but without having to store them for each event. h = da.bins.concat().hist(coarse_d=1000) i_peaks, _ = find_peaks( From a0a4715c53f7afdafcf99fa763e6c9e2b3c1fc71 Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Tue, 24 Mar 2026 14:18:42 +0100 Subject: [PATCH 07/12] refactor: better name for frame cutoff time coordinate --- src/ess/beer/clustering.py | 2 +- src/ess/beer/conversions.py | 8 ++++---- src/ess/beer/io.py | 6 +++--- 3 files changed, 8 insertions(+), 8 deletions(-) diff --git a/src/ess/beer/clustering.py b/src/ess/beer/clustering.py index 64abb6ab..a6c931c0 100644 --- a/src/ess/beer/clustering.py +++ b/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', 'tc'): + 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) diff --git a/src/ess/beer/conversions.py b/src/ess/beer/conversions.py index 0820f58c..b2b4fdba 100644 --- a/src/ess/beer/conversions.py +++ b/src/ess/beer/conversions.py @@ -43,7 +43,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.bins.coords['tc'], + da.bins.coords['frame_cutoff_time'], ) for _ in range(15): s, t0 = _linear_regression_by_bin(sin_theta_L, t, da.data) @@ -135,14 +135,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) diff --git a/src/ess/beer/io.py b/src/ess/beer/io.py index 91282dbe..b333a8f7 100644 --- a/src/ess/beer/io.py +++ b/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'] From 01672d20994f89569f82cc2c077f494e922149fd Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Tue, 24 Mar 2026 15:20:21 +0100 Subject: [PATCH 08/12] test: add test reducing 3D dataset using automatic peak finder workflow --- src/ess/beer/data.py | 5 +++++ tests/beer/mcstas_reduction_test.py | 22 ++++++++++++++++++---- 2 files changed, 23 insertions(+), 4 deletions(-) diff --git a/src/ess/beer/data.py b/src/ess/beer/data.py index 21cd6dc5..114bfec6 100644 --- a/src/ess/beer/data.py +++ b/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') diff --git a/tests/beer/mcstas_reduction_test.py b/tests/beer/mcstas_reduction_test.py index 40caec52..93dfd5f3 100644 --- a/tests/beer/mcstas_reduction_test.py +++ b/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 from ess.reduce.nexus.types import Filename, SampleRun from ess.reduce.time_of_flight.types import TofDetector @@ -41,10 +47,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',), @@ -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'), ) From 715c38983ab2c78a597df47621f21351e87b26e5 Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Tue, 24 Mar 2026 15:25:14 +0100 Subject: [PATCH 09/12] refactor: better name for chopper delay --- src/ess/beer/conversions.py | 25 ++++++++++++------------- 1 file changed, 12 insertions(+), 13 deletions(-) diff --git a/src/ess/beer/conversions.py b/src/ess/beer/conversions.py index b2b4fdba..4afbc248 100644 --- a/src/ess/beer/conversions.py +++ b/src/ess/beer/conversions.py @@ -18,7 +18,7 @@ def compute_tof_in_each_cluster( da: StreakClusteredData[RunType], - time0: WavelengthDefinitionChopperDelay, + chopper_delay: WavelengthDefinitionChopperDelay, mod_period: ModulationPeriod, ) -> TofDetector[RunType]: """Fits a line through each cluster, the intercept of the line is t0. @@ -55,7 +55,7 @@ def compute_tof_in_each_cluster( too_far_from_center=(distance_to_self > max_distance_from_streak_line), ) - t0 = _round_t0_to_nearest_chopper_opening(sc.values(t0), mod_period, time0) + 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 @@ -64,15 +64,14 @@ def compute_tof_in_each_cluster( def _round_t0_to_nearest_chopper_opening( t0: sc.Variable, mod_period: sc.Variable, - time0: sc.Variable, + chopper_delay: sc.Variable, ) -> sc.Variable: - out = t0.copy() - out -= time0 + out = t0 - chopper_delay out /= mod_period out += 0.5 sc.floor(out, out=out) out *= mod_period - out += time0 + out += chopper_delay return out @@ -152,13 +151,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' @@ -167,12 +166,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 @@ -185,7 +184,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: @@ -217,7 +216,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, From 80863ab0f5e224e0ce543e08e4f899b6762ce182 Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Wed, 25 Mar 2026 11:50:22 +0100 Subject: [PATCH 10/12] docs: explain rounding step --- src/ess/beer/conversions.py | 1 + 1 file changed, 1 insertion(+) diff --git a/src/ess/beer/conversions.py b/src/ess/beer/conversions.py index 4afbc248..6901eccb 100644 --- a/src/ess/beer/conversions.py +++ b/src/ess/beer/conversions.py @@ -33,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( From 2cfcc836b2b324670b539ac0354571246c2db079 Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Thu, 26 Mar 2026 09:11:30 +0100 Subject: [PATCH 11/12] docs: motivation for clamping origin times --- src/ess/beer/conversions.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/ess/beer/conversions.py b/src/ess/beer/conversions.py index 6901eccb..0ad7f306 100644 --- a/src/ess/beer/conversions.py +++ b/src/ess/beer/conversions.py @@ -56,6 +56,12 @@ def compute_tof_in_each_cluster( too_far_from_center=(distance_to_self > max_distance_from_streak_line), ) + # The t0 estimate from fitting is influenced by peak overlap, background, + # and other factors that can make the estimate noisy and biased. + # 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)) From 8250cc0cd50d2c36773ca734f83e91c2791d8a49 Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Thu, 26 Mar 2026 09:20:29 +0100 Subject: [PATCH 12/12] docs: reformulate --- src/ess/beer/conversions.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/ess/beer/conversions.py b/src/ess/beer/conversions.py index 0ad7f306..654c2d4a 100644 --- a/src/ess/beer/conversions.py +++ b/src/ess/beer/conversions.py @@ -57,7 +57,8 @@ def compute_tof_in_each_cluster( ) # The t0 estimate from fitting is influenced by peak overlap, background, - # and other factors that can make the estimate noisy and biased. + # 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