Update of WFM tof workflow to take into account DREAM choppers #560

nvaytet · 2024-10-14T09:50:32Z

The DREAM pulse-shaping choppers (PSC) are different from the V20 or ODIN WFM choppers in the way that they have 8 cutouts that can shape the pulse in different ways.

In the example McStas simulation that we have, the detector sees 2 frames made from 8 cutouts.
The workflow that computed time-of-flight from event-time-offset that we had assumed that the number of frames was equal to number of cutouts.

This PR updates the workflow to handle both cases.
Note that we re-wrote the workflow graph to simplify it, from

to

…rames

SimonHeybrock

Big congrats on getting this to work! 👍

I have only looked at the code changes, not the notebooks. Maybe a walkthrough would be best?

SimonHeybrock · 2024-10-15T04:46:50Z

src/scippneutron/tof/unwrap.py

+    -----
+    To find the time-of-flight origin, we ray-trace the fastest and slowest neutrons of
+    the subframes back to the first chopper to determine the time-of-flight origin.
+    The assumption here is that the first chopper is one of the two wfm choppers.


Is this the case for all ESS instruments with WFM?

SimonHeybrock · 2024-10-15T04:50:04Z

src/scippneutron/tof/unwrap.py

+    # Ray-trace back to the position of the first chopper (note that frame 0 is the
+    # pulse itself)
+    distance = frames[1].distance
+    dist = ltotal - distance
+    t1 = tmin - dist * chopper_cascade.wavelength_to_inverse_velocity(wmin)
+    t2 = tmax - dist * chopper_cascade.wavelength_to_inverse_velocity(wmax)


chopper_cascade can propagate forwards, e.g., to the detector position. Can't the same code be used to go backwards?

That would be great, but I wasn't immediately able to figure that out. Would be super if you can help with that!

I think you could just use

scippneutron/src/scippneutron/tof/chopper_cascade.py

Line 28 in 17f727d

def propagate_times(

Or use

scippneutron/src/scippneutron/tof/chopper_cascade.py

Line 87 in 17f727d

def propagate_by(self, distance: sc.Variable) -> Subframe:

if we can have a proper frame setup?

SimonHeybrock · 2024-10-15T04:52:12Z

tests/tof/unwrap_test.py

-    # We currently have no way of detecting which cutouts of the "source" chopper
-    # are used (not blocked by other choppers), for this test we remove those we
-    # know are not used.


Maybe replace by comment that explains (and assertion) that there are extra openings but the code knows how to deal with them automatically?

SimonHeybrock · 2024-10-16T04:54:23Z

docs/user-guide/wfm/dream-wfm.ipynb

+   "source": [
+    "# DREAM in WFM mode\n",
+    "\n",
+    "In this notebook, we show how to use Scippneutron's `tof` module to find the boundaries of the WFM frames for the DREAM instrument.\n",


Suggested change

"In this notebook, we show how to use Scippneutron's `tof` module to find the boundaries of the WFM frames for the DREAM instrument.\n",

"In this notebook, we show how to use ScippNeutron's `tof` module to find the boundaries of the WFM frames for the DREAM instrument.\n",

SimonHeybrock · 2024-10-16T04:59:59Z

docs/user-guide/wfm/dream-wfm.ipynb

+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "psc1 = DiskChopper(\n",


Is this setup for the case of overlapping subframes, as in the McStas data? I have not seen how that is being handled?

SimonHeybrock · 2024-10-16T05:02:03Z

src/scippneutron/chopper/disk_chopper.py

+        start_rotation:
+            This can be used to include rotations before the first pulse by giving a
+            negative integer. This is useful for choppers that are not in phase with
+            the source. The default is 0.


Is there ever a case where we don't want to include such rotations? That is, is this argument even needed?

Changed this to now always have the -1 rotation.

SimonHeybrock · 2024-10-16T05:03:18Z

src/scippneutron/chopper/disk_chopper.py

-        n_repetitions:
-            Return this many times for each angle corresponding to multiple rotations
-            of the chopper.
+        start_rotation:
+            This can be used to include rotations before the first pulse by giving a
+            negative integer. This is useful for choppers that are not in phase with
+            the source. The default is 0.
+        end_rotation:
+            How many chopper rotations to perform.


Same question as above: Maybe the existing argument was fine and we should just ensure that we always include the start?

SimonHeybrock · 2024-10-16T05:07:38Z

src/scippneutron/tof/chopper_cascade.py

+                # If end is after current end, there is only partial overlap, so we clip
+                if bound.end > current.end:
+                    old_end = current.end
+                    old_wav_end = current.wav_end
+                    current = Bound(
+                        current.start,
+                        bound.start,
+                        current.wav_start,
+                        bound.wav_start,
+                    )
+                    merged_bounds.append(current)
+                    current = Bound(
+                        old_end,
+                        bound.end,
+                        old_wav_end,
+                        bound.wav_end,
+                    )
+                # If end is before current end, overlap is total and we merge


I don't get this. What is being clipped? What is meant be "overlap is partial" and why don't we just want to merge?

nvaytet · 2024-11-04T18:30:59Z

Note that CI will fail until scipp/scipp#3582 is merged and released.

nvaytet · 2024-11-04T18:31:32Z

@SimonHeybrock I believe this is ready for another look.

SimonHeybrock · 2024-11-05T07:59:45Z

src/scippneutron/tof/chopper_cascade.py

+        # Start at -1 to ensure a rotation that is finishing when the pulse begins is
+        # also included.


Not sure what this comment relates to?

Removed. Was left over from previous changes.

SimonHeybrock · 2024-11-05T08:01:50Z

src/scippneutron/tof/unwrap.py

+    sample. The frame bounds are then computed from this.
+    """
+    return FrameAtDetector(frames[-1].propagate_to(ltotal))
+    # return FrameAtDetector(frames[ltotal])


Commented code?

Yes, I actually left it in there intentionally because I wasn't sure if I should use the getitem syntax.
I believe I needed to make the change because inside the getitem it assumes that the distance passed is a scalar.

Should I fix in the getitem or should I keep the code I have now?

To fix in the getitem, we would have to handle cases where you propagate to a range of distances, some might be before the last frame, some might be after... It gets messy.
I could also just raise a helpful error message if the distance passed in the getitem is not a scalar?

I don't have a strong opinion, do what you think works best.

src/scippneutron/tof/unwrap.py

SimonHeybrock · 2024-11-05T08:05:46Z

src/scippneutron/tof/unwrap.py

+    In the case of multiple detector pixels, we find the pixel closest to the source
+    and use that as the reference for the clipping.


Hmm, shouldn't you get more overlap for pixels with the greatest distance?

I think it doesn't matter because it will just propagate to the other distances and there are no more choppers in the way.

My point was that if there are two subframes that to not overlap at L2=5 meter, they may overload at any greater L2. Or do I misunderstand the comment?

Yes, you're right, the frames may not overlap on the closest pixel.
However, if we take the overlapping times at the farthest pixel, and use those to make a fake chopper, those times would not match the overlap seen at the closest pixel, they would be shifted to the right.
So the closest frame would have a small region where there is still overlap.

Basically, I think I need to find the overlap at the farthest pixel, and propagate those backwards to the closest pixel distance, and then make my fake chopper...

As discussed, this fake chopper implementation is not great.
@jokasimr and I do have a prototype of an alternative method which is working.
Question is: should we spend time on fixing this, or should we just try to get that into a production shape instead of trying to fix this implementation?

We could leave this one as is and assume this edge case won't be an issue in the next couple of weeks.
We then can roll out the alternative version a short time after?

Why is it an edge case? Doesn't it affect all (WFM) beamlines that do not have pixels on a sphere? Or do you mean the need for clipping is an edge case (unusual DREAM settings?).

I mean that both the DREAM choppers are a little unusual, and that the case where some pixels have overlap and others don't is hopefully an edge case as the range of pixel distances is not too large for a given instrument?
(I guess Loki is the odd one out, but they don't have WFM choppers)

Then we should just use the furthest right (and live with throwing away a tiny bit too many events for the closer pixels)? Otherwise you always get overlap for all but the closest pixel.

Let's find a time to chat in a call.

Conclusion from the call: because we are making a gap when adding the fake chopper rather than clipping with an infinitely thin vertical line, the distance it would take for the gap to completely fill in again with overlap is hopefully larger than the range of pixel distances.
So the overlap should be removed for all pixels (in most cases).

SimonHeybrock · 2024-11-05T08:06:41Z

src/scippneutron/tof/unwrap.py

-    of the WFM choppers.
+    Check for time overlap between subframes.
+    If overlap is found, we clip away the regions where there is overlap.
+    This is done by adding a fake chopper which is closed during the overlapping times.


Not clear to me why we still need choppers at this point. Instead of adding a fake chopper, can't we just clip the (sub)frame bounds?

SimonHeybrock · 2024-11-05T08:07:41Z

src/scippneutron/tof/unwrap.py

+    # We cannot chop frames at multiple distances at once, so instead of chopping the
+    # frame at the detector, we chop the last frame in the chopper cascade before the
+    # detector.
+    last_frame = chopper_cascade_frames[-1]
+
+    # Chop the subframes one by one
+    # TODO: We currently need to chop them one by one as the `chop` method seems
+    # to not work correctly with overlapping subframes.
+    subframes = []
+    for subframe in last_frame.subframes:
+        f = chopper_cascade.Frame(distance=last_frame.distance, subframes=[subframe])
+        chopped = f.chop(fake_chopper)
+        subframes.extend(
+            [
+                sf
+                for sf in chopped.subframes
+                if not sc.allclose(sf.start_time, sf.end_time)
+            ]
+        )


All this (and above) looks really complicated (read: very difficult to verify). I wonder if it is due to the fake-chopper. Why can't we just compute bounds, clip, and continue?

I am using chop because I need to have a well-constructed frame, that I can then propagate backwards to find the time origin.
I could just clip the time range, but I also want to clip the wavelengths in the correct way, so the backward propagation is correct. I thought the best way was to use chop.

Do you have a another suggestion?

The frame-propagation is already a free function. Wouldn't it still be easier to compute time and wavelength bounds directly?

Finding the time bounds is easy enough (it's just about sorting the time bounds along the subframe dimension).
Do you have a vectorized way of clipping the subframe to find the wavelengths?

The code in the _chop() function does not handle this, and I can't see an easy way to adapt it...

I honestly do not remember, it has been too long since I wrote this. We can have a look together one day if you like?

Hmm I think I may have a way to do this, by taking the mean of the intersection from all vertices of the subframe.

I've looked into simplifying this without using the fake chopper, for many pixels in a vectorized way. It's proving difficult.

My suggestion would be to merge this implementation so that it doesn't block all the other workflows.

I will revisit right after with a simpler approach but it should only be internal refactoring and shouldn't change how it is used by workflows.

SimonHeybrock · 2024-11-05T08:10:06Z

src/scippneutron/tof/unwrap.py

+    # TODO: what dimension order is the best here? Where should `detector_number` be?
+    origin_time = sc.DataArray(shift.transpose(times.dims), coords={'subframe': times})


If this is meant for processing the raw data, detector_number should be first (slowest dim), since that is how data is sorted.

SimonHeybrock · 2024-11-05T08:10:34Z

src/scippneutron/tof/unwrap.py

+        # TODO: Can we do without making a copy of delta?
+        delta = sc.lookup(delta.copy(deep=True), dim='subframe')[time_offset]


Shouldn't/isn't sc.lookup doing this?

Yes, I think you're right.

Co-authored-by: Simon Heybrock <[email protected]>

nvaytet · 2024-11-12T13:23:44Z

@SimonHeybrock any final comments?

nvaytet added 16 commits October 9, 2024 22:19

fix integer phase, add a -1 rotation, make gaps between overlapping f…

560d68d

…rames

always return more than one rotation flattened in the correct order

1131aad

add step to compute all frames in workflow

c271902

start getting time of flight origin from subframes

5867ed9

use frames to get time origin instead of reading from choppers again

3757b1f

add dream WFM example notebook

ffefdb8

fix bug in time of flight origin frame selection

83e332e

add dream notebook to index

d0f2bbd

remove unused types and providers

1b32219

add some notes

fb672de

update notebooks

7f9abc1

cleanup

57e5e15

cleanup disk_chopper

7bdf2af

update disk chopper interface to handle previous rotations

0b1db14

static analysis

7705e23

update unwrap wfm test

b8781fe

nvaytet marked this pull request as ready for review October 14, 2024 15:37

update tof/wfm test

7f7d195

SimonHeybrock reviewed Oct 15, 2024

View reviewed changes

SimonHeybrock reviewed Oct 16, 2024

View reviewed changes

nvaytet mentioned this pull request Oct 17, 2024

Fix propagation of chopper cascade frames to an array of distances #561

Merged

nvaytet added 9 commits October 17, 2024 15:57

always perform one chopper rotation before start of pulse

adfe837

fix some chopper tests

d2c1161

fix remaining disk chopper tests

9c32eb3

cleanup

1aadda1

start adding method to clip frames

625cc9b

filter out empty subframes

19861ae

add clipping of overlapping frames

23d5d50

partially fix broadcasting to multiple pixels in unwrap

ce0ca39

Merge branch 'main' into dream-wfm

48f257d

nvaytet added 2 commits November 4, 2024 18:01

Merge branch 'main' into dream-wfm

e901323

add cases with multiple pixels and overlap in dream notebook

b1117eb

nvaytet marked this pull request as ready for review November 4, 2024 18:31

nvaytet requested a review from SimonHeybrock November 4, 2024 18:31

SimonHeybrock reviewed Nov 5, 2024

View reviewed changes

nvaytet and others added 8 commits November 5, 2024 11:23

Update src/scippneutron/tof/unwrap.py

8fd24fa

Co-authored-by: Simon Heybrock <[email protected]>

remove outdated comment, and resolve TODOs

d705305

fix event_time_offset in ess fake beamline

4fe63f1

cleanup commented code

a6ca1ab

Merge branch 'main' into dream-wfm

7794863

Apply automatic formatting

e1e26b6

remove unused import

447cd88

update scipp to 24.11

9971aac

nvaytet mentioned this pull request Nov 11, 2024

Bump scipp from 24.9.1 to 24.11.0 in /requirements #570

Closed

nvaytet added 9 commits November 11, 2024 10:49

fix remaining tests

2baf93f

beam_angle -> beam_position in dream notebook

53e8676

sync buildconfig to requirements files

5722c5f

make solvable env

d6ea1cd

update deps with scipp==24.11.1

af6ae21

sync conda deps versions with pip requirements

ba555e8

fix last remaining beam_angle -> beam_position

eadfe59

fix subframe bounds

2800222

put diagram in code block

7685e7a

nvaytet and others added 2 commits November 12, 2024 14:36

add comment about gap between subframes in docstring

96f9a60

Merge branch 'main' into dream-wfm

e341c80

nvaytet merged commit 7263415 into main Nov 13, 2024
5 checks passed

nvaytet deleted the dream-wfm branch November 13, 2024 13:02

	"In this notebook, we show how to use Scippneutron's `tof` module to find the boundaries of the WFM frames for the DREAM instrument.\n",
	"In this notebook, we show how to use ScippNeutron's `tof` module to find the boundaries of the WFM frames for the DREAM instrument.\n",

		# Start at -1 to ensure a rotation that is finishing when the pulse begins is
		# also included.

		In the case of multiple detector pixels, we find the pixel closest to the source
		and use that as the reference for the clipping.

		# TODO: what dimension order is the best here? Where should `detector_number` be?
		origin_time = sc.DataArray(shift.transpose(times.dims), coords={'subframe': times})

		# TODO: Can we do without making a copy of delta?
		delta = sc.lookup(delta.copy(deep=True), dim='subframe')[time_offset]

Update of WFM tof workflow to take into account DREAM choppers #560

Update of WFM tof workflow to take into account DREAM choppers #560

Conversation

nvaytet commented Oct 14, 2024

SimonHeybrock left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

SimonHeybrock Oct 15, 2024 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

nvaytet commented Nov 4, 2024

nvaytet commented Nov 4, 2024

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

nvaytet Nov 12, 2024 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

nvaytet commented Nov 12, 2024 • edited Loading

SimonHeybrock Oct 15, 2024 •

edited

Loading

nvaytet Nov 12, 2024 •

edited

Loading

nvaytet commented Nov 12, 2024 •

edited

Loading