[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.

[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.

[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.