Skip to content

[essreduce] Account for slow neutrons in LUT chopper frame sequence #751

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
nvaytet wants to merge 3 commits into
base: main
Choose a base branch
from
Open

[essreduce] Account for slow neutrons in LUT chopper frame sequence #751

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
0829440
rotate the choppers for long enough to cover the entire range of poss…
nvaytet Sep 15, 2026
858e545
add test
nvaytet Sep 15, 2026
731b8d4
fix test
nvaytet Sep 15, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
60 changes: 43 additions & 17 deletions packages/essreduce/src/ess/reduce/unwrap/lut.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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')
Comment on lines +29 to +35

@SimonHeybrock SimonHeybrock Sep 15, 2026 •
edited
Loading

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Could we avoid the magic 500 m? The chopper opening times only matter where the frame is chopped, i.e. at the chopper positions. Beyond the last chopper the frame is only propagated, and wrapping at the detector is handled by the period copies in _estimate_wavelength_by_polygon_centers. So the time range to cover ends when the slowest neutron of the last source pulse reaches the farthest chopper:

travel_time = (
    source_bounds.time[1]
    + (pulse_stride - 1) * pulse_period
    + max_chopper_distance / _wavelength_to_speed(source_bounds.wavelength[1])
)

The (pulse_stride - 1) term is needed because from_source_pulse(npulses=pulse_stride) creates later pulses at i * pulse_period. With 500 m this term is hidden by the margin. This also reduces the rotation count a lot, e.g. about 3 pulse periods instead of 27 for the NMX-like setup in the new test.

Edit: The subframes do keep spreading out after the last chopper, but that needs more polygon copies (line 649), not more chopper rotations. propagate_to only shifts the polygon vertices by d / v, and the chopper opening times are not used after the last chopper. At 300 m the surviving subframe of the setup in the new test spans pulse periods 1.65 to 3.86.

I checked this numerically with the choppers from the new test and LtotalRange 60-300 m. The table built with the farthest chopper as the distance (plus the +1 rotation from the other comment) is bit-identical to the table built with 500 m.

The +1 is required: DiskChopper starts its repetitions at rotation -1, so n repetitions only give openings up to about (n - 1) / f. With the current formula and 52 m instead of 500 m (3 periods), the 12-15 Å band is missing from the table. With 70 m (4 periods) the table is identical to the one built with 500 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:
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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

@SimonHeybrock SimonHeybrock Sep 15, 2026

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

nperiods is computed from the absolute max_time, but the copies are shifted by noffset + i. So the first noffset extra copies end up at negative times and only contribute NaNs. This is correct, but for long flight paths it adds work in the per-distance loop. int(max_time / frame_period) - noffset + 1 would be sufficient.

Also, no test covers this change: with range(nperiods) reverted to (0, 1) the new test still passes. Could the new test also compute the LookupTable and check that the 12-15 Å band shows up at the detector?


polygons = [
np.stack(
[
Expand All @@ -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)
Expand Down Expand Up @@ -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)
Comment on lines +701 to +705

@SimonHeybrock SimonHeybrock Sep 15, 2026

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

time_offset_open requires |chopper.frequency| / pulse_frequency to be an integer or the inverse of an integer. With the defaults, nperiods = ceil((5 ms + 500 m / v(15 Å)) * 14 Hz) = 27. A 7 Hz pulse-skipping chopper then gives a ratio of 27/2 and raises. The comment removed here explained exactly this, and why the old code divided by an even number.

Suggestion: pick the rotation count per chopper, from that chopper's own frequency. Then the ratio is an integer by construction:

freq = abs(ch.frequency).to(unit='Hz')
nrot = int(np.ceil((travel_time * freq).value)) + 1
time_open = ch.time_offset_open(pulse_frequency=freq / nrot)

I tried this locally: all tests in tests/unwrap pass, including your new test and the unmodified version of test_lut_does_not_raise_if_no_neutrons_make_it_through.

Note that this no longer rejects any chopper frequency. Frequencies that do not match the source would then need an explicit check, see the review body.


chops = {
key: chopper_cascade.Chopper(
Expand Down Expand Up @@ -714,7 +735,7 @@ def make_wavelength_lut_from_polygons(
time_resolution: TimeResolution,
pulse_period: PulsePeriod,
pulse_stride: PulseStride[RunType],
frames: ChopperFrameSequence,
frames: ChopperFrameSequence[RunType],

@nvaytet nvaytet Sep 15, 2026

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This was a typo; I don't know how it still worked...

) -> LookupTable[RunType, Component]:
"""
Compute a lookup table for wavelength as a function of distance and
Expand Down Expand Up @@ -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

@nvaytet nvaytet Sep 15, 2026

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is semi-unrelated, but I discovered it by messing around on the workflow while looking at the choppers.

# 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)

@SimonHeybrock SimonHeybrock Sep 15, 2026

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

LtotalRange is documented as (min, max). If someone sets it the wrong way round, that is probably a mistake in their setup. I would rather raise a ValueError than silently swap the values. Either way, this could be a separate PR.

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.
Expand Down
60 changes: 59 additions & 1 deletion packages/essreduce/tests/unwrap/lut_test.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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'),
Expand All @@ -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
Comment on lines +411 to +412

@SimonHeybrock SimonHeybrock Sep 15, 2026

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This change to the test is a symptom of the phase issue above. A 0.1 Hz chopper with a 14 Hz source works on main, and changing the source period to 10 s changes what the test covers.

If we adopt the explicit frequency check from the review body, a 0.1 Hz chopper would (correctly) raise. The test could then block the beam with a 14 Hz chopper whose opening is out of phase with the pulse.

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.
Expand Down Expand Up @@ -487,3 +490,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')
Loading