From 0829440bfd178e5120afd3bd962052030f6fc9eb Mon Sep 17 00:00:00 2001 From: Neil Vaytet Date: Tue, 15 Sep 2026 13:15:25 +0200 Subject: [PATCH 1/3] rotate the choppers for long enough to cover the entire range of possible arrival times --- .../essreduce/src/ess/reduce/unwrap/lut.py | 58 ++++++++++++++----- 1 file changed, 42 insertions(+), 16 deletions(-) diff --git a/packages/essreduce/src/ess/reduce/unwrap/lut.py b/packages/essreduce/src/ess/reduce/unwrap/lut.py index 13d4fdd6e..bcd6aacb7 100644 --- a/packages/essreduce/src/ess/reduce/unwrap/lut.py +++ b/packages/essreduce/src/ess/reduce/unwrap/lut.py @@ -26,6 +26,26 @@ WavelengthLutMode, ) +# We define a maximum instrument length which is used to determine how many chopper +# rotations should be performed when computing the chopper frame sequence. +# We need to rotate the choppers for long enough to make sure we capture cases where +# very slow neutrons pass through chopper openings multiple pulse periods later. +# The most robust way is to define the longest possible distance that could be traveled +# and compute how long it would take the slowest neutrons to reach it. +MAXIMUM_INSTRUMENT_LENGTH = sc.scalar(500.0, unit='m') + + +def _wavelength_to_speed(wavelength: sc.Variable) -> sc.Variable: + """ + Convert wavelength to speed. + + Parameters + ---------- + wavelength: + Wavelength of the neutrons. + """ + return (sc.constants.h / sc.constants.m_n) / wavelength + @dataclass class BeamlineComponentReading: @@ -58,7 +78,7 @@ class BeamlineComponentReading: distance: sc.Variable def __post_init__(self): - self.speed = (sc.constants.h / sc.constants.m_n) / self.wavelength + self.speed = _wavelength_to_speed(self.wavelength).to(unit='m/s') @dataclass @@ -623,6 +643,11 @@ def _estimate_wavelength_by_polygon_centers( # This is because neutrons that arrive after the frame period will wrap around and # appear in the next pulse, which is equivalent to the original pulse but shifted # by the frame period. + # We determine the number of frame periods to shift by calculating how many periods + # are needed to cover the maximum arrival time in the subframes. + max_time = sc.reduce([f.time.max() for f in subframes]).max() + nperiods = int(max_time.to(unit=time_unit).value / frame_period.value) + 1 + polygons = [ np.stack( [ @@ -632,7 +657,7 @@ def _estimate_wavelength_by_polygon_centers( axis=1, ) for f in subframes - for i in (0, 1) + for i in range(nperiods) ] wavs, stddevs = _polygon_intersections(polygons, time_edges.values) @@ -670,18 +695,14 @@ def compute_frame_sequence( # The `pulse_frequency` parameter in time_offset_open and time_offset_close below # decides how many rotations the chopper will perform when computing the open and - # close times. Because we want to cover a number of pulses equal to `pulse_stride`, - # we need to set the pulse frequency to be `pulse_stride` times smaller than the - # actual pulse frequency. - # - # In addition, the time_offset_open and time_offset_close below require the - # pulse_frequency to be an integer multiple of the pulse frequency or vice versa. - # A simple trick is to make sure that the requested pulse frequency is divided by - # an even number. We need to rotate the chopper for long enough to cover wrapping - # around the frame period, so we cover two pulses strides. - frequency_for_chopper_rotation = (1.0 / pulse_period.to(unit='s')) / ( - pulse_stride * 2 - ) + # close times. + # We need to cover the entire time range from 0 to the time it takes the slowest + # neutron to travel the maximum instrument length. + travel_time = source_bounds.time[1].to(unit='s') + ( + MAXIMUM_INSTRUMENT_LENGTH / _wavelength_to_speed(source_bounds.wavelength[1]) + ).to(unit='s') + nperiods = sc.ceil(travel_time / pulse_period) + frequency_for_chopper_rotation = 1.0 / (nperiods * pulse_period) chops = { key: chopper_cascade.Chopper( @@ -744,8 +765,13 @@ def make_wavelength_lut_from_polygons( pulse_period = pulse_period.to(unit=time_unit) frame_period = pulse_period * pulse_stride - min_dist = ltotal_range[0].to(unit=distance_unit) - max_dist = ltotal_range[1].to(unit=distance_unit) + dist0 = ltotal_range[0].to(unit=distance_unit) + dist1 = ltotal_range[1].to(unit=distance_unit) + # By default, the minimum and maximum distances should be the first and second + # elements of the total range. But if the user set them manually on the workflow + # we need to make sure we pick the minimum and maximum distances. + min_dist = min(dist0, dist1) + max_dist = max(dist0, dist1) # We want to give the 2d interpolator a table that covers the requested range, # hence we need to extend the range by at least half a resolution in each direction. From 858e54591c86dd15bbbdec851b1a17457e845bec Mon Sep 17 00:00:00 2001 From: Neil Vaytet Date: Tue, 15 Sep 2026 13:47:30 +0200 Subject: [PATCH 2/3] add test --- .../essreduce/src/ess/reduce/unwrap/lut.py | 2 +- packages/essreduce/tests/unwrap/lut_test.py | 55 +++++++++++++++++++ 2 files changed, 56 insertions(+), 1 deletion(-) diff --git a/packages/essreduce/src/ess/reduce/unwrap/lut.py b/packages/essreduce/src/ess/reduce/unwrap/lut.py index bcd6aacb7..1cfb1bff2 100644 --- a/packages/essreduce/src/ess/reduce/unwrap/lut.py +++ b/packages/essreduce/src/ess/reduce/unwrap/lut.py @@ -735,7 +735,7 @@ def make_wavelength_lut_from_polygons( time_resolution: TimeResolution, pulse_period: PulsePeriod, pulse_stride: PulseStride[RunType], - frames: ChopperFrameSequence, + frames: ChopperFrameSequence[RunType], ) -> LookupTable[RunType, Component]: """ Compute a lookup table for wavelength as a function of distance and diff --git a/packages/essreduce/tests/unwrap/lut_test.py b/packages/essreduce/tests/unwrap/lut_test.py index ede8b0d28..3789a6de5 100644 --- a/packages/essreduce/tests/unwrap/lut_test.py +++ b/packages/essreduce/tests/unwrap/lut_test.py @@ -487,3 +487,58 @@ def test_polygon_intersections_handles_uncovered_columns_without_warning(): # Columns 0 and 2 miss the polygon (all-NaN); column 1 is covered. np.testing.assert_array_equal(np.isnan(center), [True, False, True]) np.testing.assert_array_equal(np.isnan(spread), [True, False, True]) + + +def test_choppers_rotate_enough_times_to_catch_slow_neutrons(): + choppers = { + "chopper1": DiskChopper( + axle_position=sc.vector([0, 0, 28.4], unit='m'), + frequency=sc.scalar(-14, unit='Hz'), + beam_position=sc.scalar(0, unit='deg'), + phase=sc.scalar(-112.3, unit='deg'), + slit_begin=sc.array(dims=["cutout"], values=[-38.5], unit='deg'), + slit_end=sc.array(dims=["cutout"], values=[38.5], unit='deg'), + slit_height=None, + radius=None, + ), + "chopper2a": DiskChopper( + axle_position=sc.vector([0, 0, 50.9774], unit='m'), + frequency=sc.scalar(-14, unit='Hz'), + beam_position=sc.scalar(0, unit='deg'), + phase=sc.scalar(194.1, unit='deg'), + slit_begin=sc.array(dims=["cutout"], values=[-70.0], unit='deg'), + slit_end=sc.array(dims=["cutout"], values=[70.0], unit='deg'), + slit_height=None, + radius=None, + ), + "chopper2b": DiskChopper( + axle_position=sc.vector([0, 0, 51.0024], unit='m'), + frequency=sc.scalar(-14, unit='Hz'), + beam_position=sc.scalar(0, unit='deg'), + phase=sc.scalar(168.0, unit='deg'), + slit_begin=sc.array(dims=["cutout"], values=[-70.0], unit='deg'), + slit_end=sc.array(dims=["cutout"], values=[70.0], unit='deg'), + slit_height=None, + radius=None, + ), + } + wf = _make_workflow("analytical") + wf[unwrap.DiskChoppers[AnyRun]] = choppers + wf[Position[snx.NXsource, AnyRun]] = sc.vector([0, 0, 0], unit='m') + + frames = wf.compute(unwrap.ChopperFrameSequence[AnyRun]) + + # In this configuration (based on the NMX instrument), the last frame should have + # two subframes: a main subframe containing short wavelengths 1-5 Å and a secondary + # subframe containing longer wavelengths 12-15 Å. + last_frame = frames[-1] + assert len(last_frame.subframes) == 2 + main_subframe = last_frame.subframes[0] + secondary_subframe = last_frame.subframes[1] + + # Check the wavelength ranges for the subframes + assert main_subframe.wavelength.min() > sc.scalar(1, unit='angstrom') + assert main_subframe.wavelength.max() < sc.scalar(5, unit='angstrom') + + assert secondary_subframe.wavelength.min() > sc.scalar(12, unit='angstrom') + assert secondary_subframe.wavelength.max() < sc.scalar(15, unit='angstrom') From 731b8d4ac52558b70db81b1eee158647452e3f81 Mon Sep 17 00:00:00 2001 From: Neil Vaytet Date: Tue, 15 Sep 2026 13:56:59 +0200 Subject: [PATCH 3/3] fix test --- packages/essreduce/tests/unwrap/lut_test.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/packages/essreduce/tests/unwrap/lut_test.py b/packages/essreduce/tests/unwrap/lut_test.py index 3789a6de5..b9f3ffaee 100644 --- a/packages/essreduce/tests/unwrap/lut_test.py +++ b/packages/essreduce/tests/unwrap/lut_test.py @@ -395,10 +395,11 @@ def test_lut_workflow_guesses_pulse_stride(): def test_lut_does_not_raise_if_no_neutrons_make_it_through(wavelength_from): wf = _make_workflow(wavelength_from) # Add a very slowly rotating chopper that will block all neutrons. + freq = sc.scalar(0.1, unit='Hz') wf[unwrap.DiskChoppers[AnyRun]] = { 'chopper1': DiskChopper( axle_position=sc.vector([0, 0, -15.0], unit='m'), - frequency=sc.scalar(0.1, unit='Hz'), + frequency=freq, beam_position=sc.scalar(0.0, unit='deg'), phase=sc.scalar(0.0, unit='rad'), slit_begin=sc.array(dims=['cutout'], values=[0.0], unit='deg'), @@ -407,6 +408,8 @@ def test_lut_does_not_raise_if_no_neutrons_make_it_through(wavelength_from): radius=sc.scalar(0.35, unit='m'), ) } + # Need to synchronize the source period with the chopper frequency. + wf[unwrap.PulsePeriod] = 1.0 / freq wf[Position[snx.NXsource, AnyRun]] = sc.vector([0, 0, -25.0], unit='m') # Need to force the pulse stride so that it doesn't get set to a large value due to # the slow chopper.