[essreduce] Fast wavelength lookup table workflow by nvaytet · Pull Request #589 · scipp/ess

nvaytet · 2026-05-18T08:39:47Z

Alternative to #557:

Instead of replacing the existing workflow with the one that uses scippneutron.tof.chopper_cascade to compute the wavelength LUT, we simply add a new workflow.

This can then be rolled out to different packages as desired.

The next step will be to use the providers from the new workflow and inject them directly into the GenericUnwrapWorkflow to compute the LUT on-the-fly. We defer that to a later PR.

Additional info (copied from #557)

To approximate the polygons as single lines, we produce a set of vertical lines that represent the bins in event_time_offset.
For each line, we find the intersections with the polygon edges, and then compute the midpoint between the min and max y coordinate of the intersections.

The table below is computed in ~0.7s compared to ~16s before (using 5M neutrons).

To use:

import scipp as sc
from ess.reduce import unwrap
from ess.reduce.nexus.types import AnyRun
from ess.odin.beamline import choppers

source_position = sc.vector([0, 0, 0], unit='m')
disk_choppers = choppers(source_position)

wf = unwrap.FastLookupTableWorkflow()
wf[unwrap.DiskChoppers[AnyRun]] = disk_choppers
wf[unwrap.SourcePosition] = source_position
wf[unwrap.PulseStride] = 2
wf[unwrap.LtotalRange] = sc.scalar(5.0, unit="m"), sc.scalar(65.0, unit="m")
wf[unwrap.DistanceResolution] = sc.scalar(0.1, unit="m")
wf[unwrap.TimeResolution] = sc.scalar(250.0, unit='us')

table = wf.compute(unwrap.LookupTable)
table.plot(vmax=10.0)

Needs scipp/scippneutron#700

github-actions · 2026-05-18T08:40:38Z

Hi! Your PR was missing some labels 🔖 so I added them: essreduce

jokasimr · 2026-05-18T10:55:18Z

+        and event_time_offset in the lookup table.
+    """
+    distance_unit = "m"
+    time_unit = "us"


why not ns, since event_time_offset typically has unit ns?

It was because scientists usually think about time-of-flight in μs, and so when looking at lookup tables, it make more sense visually when the axis is in μs.
That said, you are right that this means we are almost always making a copy of the array of event_time_offset because we comvert it to μs...

I think I will use ns and save the compute 👍

Aha I see, in that case either unit works for me

Hmm actually it's the other way round:

we take the event_time_offset unit: https://github.com/scipp/ess/blob/fast-lut/packages/essreduce/src/ess/reduce/unwrap/to_wavelength.py#L270

and convert the axis of the LUT to that unit: https://github.com/scipp/ess/blob/fast-lut/packages/essreduce/src/ess/reduce/unwrap/to_wavelength.py#L93

jokasimr · 2026-05-18T10:56:31Z

+    min_dist = ltotal_range[0].to(unit=distance_unit)
+    max_dist = ltotal_range[1].to(unit=distance_unit)
+
+    # We need to bin the data below, to compute the weighted mean of the wavelength.


What does this refer to? I don't see any binning or weighted mean here. Was the comment duplicated from the other lookup-table-computation function?

Yes, copied over from old computation.

jokasimr · 2026-05-18T11:06:56Z

+        # Find the correct simulation reading
+        selected_frame = None
+        for frame in sorted_frames:
+            if dist.value >= frame.distance.to(unit=dist.unit).value:


I don't quite understand this, what is frame here? Why are there multiple of them? Later on in this function we compute "subframes" associated with the selected frame, what are those?

The 'frame' and 'subframe' terminology is from the chopper_cascade module: https://scipp.github.io/scippneutron/user-guide/chopper/chopper-cascade.html

A frame sequence contains frames, which is a description of how the pulse looks like at different stages along the beamline. There is a frame at the source (when the pulse is just a rectangle), one frame at each chopper, and a final frame that was propagated to the detector distance.

Then, each frame consists of one or multiple Subframe which are the polygons on the figure. Basically, one Frame is depicted by one colour on the figure. And there are one or more polygons of the same colour.

jokasimr · 2026-05-18T11:10:30Z

+        pulse_period=pulse_period,
+        pulse_stride=pulse_stride,
+        distance_resolution=table.coords["distance"][1] - table.coords["distance"][0],
+        time_resolution=table.coords["event_time_offset"][1]


Is there a reason why we are not using the time_resolution and distance_resolution arguments here?

The final time_resolution may not be exactly what the value in the argument was, because we want it to fit exactly into one or more pulse periods. We give the step size which is the closest to what the user requested but fits exactly in the range (resolution will always be at least what the user requested or finer, never coarser).

The distance_resolution is preserved and in principle we could just use the argument. In this case, it is the ltotal_range which is not exactly preserved. But in case we decided later to change that behaviour in the future, I wanted to be sure that bugs were not introduced here (like forgetting to update the resolution in the table).
So I thought the safest is just to take it from the array.

jokasimr · 2026-05-18T11:41:12Z

+    # around the frame period, so we cover two pulses strides.
+    frequency_for_chopper_rotation = (1.0 / pulse_period.to(unit='s')) / (
+        pulse_stride * 2
+    )


I'd prefer if pulse_stride was not an argument to this function, because essentially the pulse stride is a function of the chopper system, and we already pass the chopper parameters as an argument. Can we derive the pulse_stride from the chopper parameters?

The initial thought was that trying to guess from chopper settings was dangerous.
Were you thinking of guessing from the rotation frequencies?

There may be chopper configurations we did not think of?

For example, can there be a chopper that instead of having one opening and spins at 14 Hz, it has 2 symmetric openings and spins at 7 Hz?

Or maybe during hot commissioning, they may operate the source at 7Hz, meaning that Odin might just decide to park the pulse skipping chopper, meaning that we find only choppers with frequencies of 14 Hz and higher, thus guessing a pulse_stride of 1?

Yes I was thinking we could guess from the chopper frequencies.

Is it a problem if we guess a "too large" pulse stride?
My understanding is that it's only a problem if we guess a "too small" pulse stride. So in the case of a 7Hz chopper with two openings wouldn't be a problem, because we would guess pulse_stride=2, while it could have been treated as 1 as well.

I don't really understand the Odin example, if all choppers are operating at 14 Hz then the $\lambda(t_{oa})$ function will also be periodic with a frequency of 14 Hz right?

jokasimr · 2026-05-18T11:45:54Z

+    return 0.5 * (y_min + y_max), 0.5 * (y_max - y_min)
+
+
+def _compute_mean_wavelength_in_polygons(


Maybe rename to _estimate_wavelength_by_polygon_centers to be more specific.

jokasimr · 2026-05-18T11:57:59Z



+@dataclass
+class SourcePulse:


Maybe rename to UniformSourcePulse to make it distinct from the non-uniform source pulse that is used in tof.

How about SourceExtent or SourceBounds?
I agree SourcePulse isn't great, but I'm not sure users will immediately understand why Uniform comes into play?

nvaytet added 3 commits May 18, 2026 10:19

add FastLookupTableWorkflow based on chopper_cascade

5377379

bump scippneutron version

3d493f0

update pixi lock file

4246775

github-actions Bot added the essreduce Issues for essreduce. label May 18, 2026

nvaytet added 2 commits May 18, 2026 11:21

add missing types in init

94cbe24

add notebook that uses new workflow

7a17b4b

nvaytet requested review from SimonHeybrock and jokasimr May 18, 2026 09:24

lint

bd62ad0

nvaytet mentioned this pull request May 18, 2026

[essreduce] Compute wavelength LUT from chopper cascade #557

Closed

jokasimr reviewed May 18, 2026

View reviewed changes

nvaytet added 2 commits May 18, 2026 14:43

remove old comments and give function a better name

2261165

SourcePulse -> SourceBounds

ba01ece

		return 0.5 * (y_min + y_max), 0.5 * (y_max - y_min)


		def _compute_mean_wavelength_in_polygons(



		@dataclass
		class SourcePulse:

Conversation

nvaytet commented May 18, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Alternative to #557:

Additional info (copied from #557)

Uh oh!

github-actions Bot commented May 18, 2026

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

jokasimr May 18, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

nvaytet commented May 18, 2026 •

edited

Loading

jokasimr May 18, 2026 •

edited

Loading