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

[nvaytet]

Computing the wavelength LUT in 'analytical' mode had a flaw where in the case of very slow neutrons making it through some much later chopper openings, the choppers were not peforming enough rotations (only 2 pulse periods) and they were blocking the slow neutrons.

The tof simulation saw these slow neutrons 'polluting' subsequent pulses

Before: the frame sequence only had a single subframe at the detector

After: two subframes at the detector

[nvaytet]

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

[nvaytet]

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

[SimonHeybrock]

This review was created with the help of AI. If anything below strikes you as unpleasantly formulated (tone, verbosity, ...), as too nit-picky, or as otherwise improvable, please tell me -- I am actively trying to improve this.

The new rotation count breaks every chopper whose frequency is not a multiple of 14 Hz / nperiods. That is why CI is red: all 16 test_pulse_skipping_*[analytical] tests in unwrap_test.py fail with "The chopper is out of phase with the source". Details and a possible fix are in the inline comments.

This raises a policy question we should settle explicitly: what should happen if a chopper runs at a frequency that does not match the source, e.g. 5 Hz? esslivedata builds the table from live setpoints, so this can happen in production. The chopper pattern then differs from frame to frame, and no static lookup table is correct. A table that lets nothing through would produce empty wavelength spectra (per decision form this morning), or a table that lets sth. arbitrary through might produce silently wrong spectra.

Should we raise, with a message naming the chopper and the condition? esslivedata reports job errors, so the problem would be visible to users. The condition should be the physical one: |f| * pulse_stride * pulse_period is an integer. The pulse_frequency trick in time_offset_open gives only an accidental check, which the PR currently changes as a side effect. Alternatively, produce a table that lets nothing pass, so detector views do not show new neutrons?

[SimonHeybrock]

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.

[SimonHeybrock]

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.

[nvaytet]

Why do we need to find the max_chopper_distance? Can't we just use the distance to the present chopper?
Each chopper will compute a different travel_time, but I think that's fine?
It would perform less rotations for choppers closer to the source, thus reducing the total number of rotations even further?

The +1 is required

I did not get that part. Which +1 are we talking about here?

[SimonHeybrock]

Good point, the per-chopper distance is enough, and it is what 428e6cf does now. Frame.chop propagates the frame to chopper.distance and intersects it with that chopper's openings only (scippneutron tof/chopper_cascade.py:216-222). A chopper therefore only needs openings up to the time the slowest neutron of the last pulse reaches it. The distances of other choppers do not matter. I rechecked on the current head with the choppers from test_choppers_rotate_enough_times_with_slow_neutrons_to_pollute_lut, LtotalRange 60-300 m and 15 Å: the per-chopper version, the max-chopper-distance version, and the old 500 m version all give bit-identical tables. The saving is small here: 3 instead of 4 rotations for the chopper at 28.4 m.

The +1 is the one in nrot = int(np.ceil((travel_time * freq).value)) + 1. DiskChopper._apply_angle_repetitions starts the repetitions at -1 (sc.arange(..., -1, n_repetitions), to catch a rotation that is still finishing when the pulse starts). So nrot repetitions give openings only up to rotation nrot - 1. Without the +1, the last rotation needed can be missing. With the same setup and ceil(...) alone, 30618 table cells become NaN, all in the 12.6-14.1 Å range. Maybe a short comment on that line would help, since the reason sits in scippneutron.

[nvaytet]

The saving is small here: 3 instead of 4 rotations for the chopper at 28.4 m.

Also that we don't need to figure out which one is the last chopper, just take the current chopper. So it's simpler.

[SimonHeybrock]

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?

[SimonHeybrock]

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.

[nvaytet]

I decided to raise instead. Not sure a whole other PR is warranted.

[SimonHeybrock]

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.

[SimonHeybrock]

Thoughts (partially beyond the scope of this PR):

Context: scipp/esslivedata#1314 builds the table from live chopper setpoints and, by team decision, substitutes a chopper that is out of phase with the source by one with no slits. The table then blanks downstream of it and the consumers keep publishing empty results instead of reducing with a stale table. It uses a per-chopper condition against the source frequency, because the provider producing DiskChoppers cannot see PulseStride (which is derived from them), and it relies on scippneutron's check as the fallback for anything it misses. The per-chopper rotation count suggested above removes that fallback, which makes the following relevant for how we would use this PR:

• An explicit phase condition in compute_frame_sequence: |f| * pulse_stride * pulse_period is a whole number. That is the only function where choppers, period, and stride are all in scope, so a check anywhere earlier runs into the cycle above.
• The stride guess must ignore choppers that cannot be locked. guess_pulse_stride_from_choppers rounds F / f for every turning chopper, so a 5 Hz chopper gives stride 3, and with the condition above DREAM's 7 Hz pulse-shaping choppers then fail it as well (this is DREAM's configuration as observed on 2026-09-10). Guessing only from choppers with a whole F / f gives 2, and only the 5 Hz chopper fails.
• Behaviour on failure as a workflow parameter: raise or block, default raise. Block means the chopper contributes no openings (empty open and close arrays). Batch reduction from a file wants the raise, live reduction from setpoints wants block. Publish a blocking lookup table when a chopper is out of phase esslivedata#1314 shows that a chopper with no openings already goes all the way to a table.
• A public function returning the names of the out-of-phase choppers for given choppers, period, and stride, so a consumer can log which chopper emptied the table. In esslivedata that would be a provider consuming the derived stride, with no cycle.
• 0 Hz: the stride guess treats a stopped chopper as inactive, while the frame sequence has it block the beam (time_offset_open gives [-inf, inf], time_offset_close gives [nan, inf]). One of the two should change. Could be a separate issue.

With those, scipp/esslivedata#1314 reduces to setting the parameter and logging.

[SimonHeybrock]

This fails for any chopper that declares slit_height, which includes DREAM's: slit_height keeps its one entry while slit_begin and slit_end become empty, and DiskChopper raises DimensionError: Expected (cutout: 0) to include (slit: 1). The tests here only use slit_height=None. The fixed cutout dim also differs from the chopper's own slit dim (slit or dim_0 in the files we load).

Something like this works in esslivedata's tests:

dim = ch.slit_begin.dim
empty = sc.array(dims=[dim], values=[], unit='deg')
height = (
    None
    if ch.slit_height is None
    else sc.array(dims=[dim], values=[], unit=ch.slit_height.unit)
)
out[key] = replace(
    ch,
    frequency=pulse_frequency,
    slit_begin=empty,
    slit_end=empty,
    slit_height=height,
)

A test with a chopper that has slit heights would catch this.

[SimonHeybrock]

compute_frame_sequence now receives ProcessedDiskChoppers, so the 0 Hz skip and the phase check here are already done by process_disk_choppers. Can this branch go, so the condition lives in one place?

[SimonHeybrock]

This is the check against the source frequency for each chopper on its own, not the check for the whole cascade. With rotations now counted per chopper, scippneutron no longer catches the rest either. Example: choppers at 14/3 Hz and 14/4 Hz each pass here, the guessed stride is 4, and the table is built without error (34951 of 120120 entries finite), although the 14/3 Hz chopper turns 4/3 times per frame.

The full condition is that |f| * pulse_stride * pulse_period is a whole number for every chopper. compute_frame_sequence is the only place where the choppers and the stride are both available, so it would have to go there. It is an unlikely configuration, so a separate issue is fine.

[nvaytet]

Not sure I followed.
My understanding was that say we have 3 choppers in the cascade, and the last one is at a weird frequency. That one becomes a blocking chopper but not the others.
We should still able to compute a LUT from the source to the blocking chopper.
If a monitor is placed before the blocking chopper, we can still compute wavelengths for that monitor.

So I'm not sure why we have to consider all choppers at once...

[SimonHeybrock]

The point is that those frequencies above are not weird -- running at a fraction of source frequency is valid, isn't it? They make sense for pulse-skipping modes, but if you have two PSCs that are set to different "skipping modes" we do not detect that.

[nvaytet]

So which way forward do you suggest here? Not sure I understood what you wanted changed.

[nvaytet]

After in-person discussion, we decided to not worry about this for now.
--> add a comment in the code about this possible edge case

[nvaytet]

@SimonHeybrock the essdiffraction tests are failing because beer has a mode with a chopper rotating at 63 Hz, which is 4.5 x 14.
Does this mean that this mode should actually be allowed as long as pulse_stride is set to 2?

[nvaytet]

Feel free to push fixes and merge while I'm away

[nvaytet]

I fixed by using the pulse stride to check if chopper is in sync or not.

It means I now need to guess the pulse stride from the raw choppers instead of the processed choppers.
Because I didn't want to drop the 0Hz choppers twice, I added an extra step in the workflow which removes the 0Hz choppers and yields ActiveDiskChoppers, which is then used by the pulse stride guess and the provider which computes FrameCompatibleDiskChoppers (which is the new name for ProcessedDiskChoppers).

[nvaytet]

For the record, this was the Beer case:

After the fix, the lookup table has the correct prediction at the detector:

[SimonHeybrock]

Checking against the frame frequency instead of the source frequency fixes BEER and also covers the 14/3 Hz + 14/4 Hz case we discussed. CI is red only because two tests still import process_disk_choppers, see inline.

[SimonHeybrock]

process_disk_choppers was split into get_active_choppers and close_non_synced_disk_choppers, so this test and the next one fail with an ImportError. These are the two failures in the essreduce CI job. The next test also has to pass a pulse_stride now.

[SimonHeybrock]

This comment is no longer correct now that the check uses the frame frequency. With 14/3 Hz and 14/4 Hz the guessed stride is 4, the frame frequency is 3.5 Hz, and the 14/3 Hz chopper has a ratio of 4/3, so it is closed. I checked this by calling guess_pulse_stride_from_choppers and close_non_synced_disk_choppers directly. The open question in the comment is fine too: 14/2 Hz with stride 4 gives a ratio of 2, so the chopper turns exactly twice per frame. I think the comment can be removed.

[SimonHeybrock]

Do we need the inverse branch? The physical condition is that each chopper turns a whole number of times per frame, i.e. |f| * frame_period is an integer. If the chopper period is instead n > 1 frame periods, consecutive frames see the chopper at different angles, so no single table is correct.

With the guessed stride this cannot happen, because stride >= round(F / f). However, esslivedata lets users set the stride by hand (auto_stride=False). With stride 1 and a 7 Hz chopper, the chopper passes this check and is not closed.

I tried an integer-only check:

if not bool(abs(sc.round(quot) - quot) < sc.scalar(1e-8)):

All tests in essreduce/tests/unwrap and essdiffraction/tests pass, apart from the two that fail with the ImportError. _is_int_or_inverse_int could then go as well, and the docstring of FrameCompatibleDiskChoppers would need updating.

[nvaytet]

I went with the suggested fix and added one more test to make sure a 7Hz chopper gets closed if pulse stride is forced to 1 by the user.

[SimonHeybrock]

Nit: this is the frame frequency, so frame_frequency would avoid confusing it with the source frequency. The docstring is also missing pulse_stride.

[SimonHeybrock]

Dropping 0 Hz choppers means a stopped chopper is assumed to be parked out of the beam. Before this PR, the frame sequence treated it as blocking. Could the docstring state this assumption?

[SimonHeybrock]

This is my review reply copied almost word for word, and it reads like one. Something shorter would do:

    # Copy i is shifted by (noffset + i) frame periods, so copies up to
    # i = int(max_time / frame_period) - noffset are needed to cover max_time.

[SimonHeybrock]

The +1 still has no explanation, and the reason is in scippneutron. Also, the integer requirement now holds by construction. Maybe:

        # time_offset_open/close require freq / pulse_frequency to be an integer,
        # which holds here by construction. DiskChopper starts its repetitions at
        # rotation -1, so nrot repetitions only reach rotation nrot - 1, hence +1.

[SimonHeybrock]

Typo:

Suggested change

	def test_lut_workflow_chopper_frquency_multiple_of_frame_period(wavelength_from):
	def test_lut_workflow_chopper_frequency_multiple_of_frame_period(wavelength_from):

[SimonHeybrock]

Thanks, all points are addressed. A few leftovers in comments and docstrings, all as suggestions that can be applied directly.