Alternative to time-origin tof computation #577

nvaytet · 2024-11-25T13:28:43Z

This is an alternative implementation of the computation of the tof coordinate in the unwrap submodule.

The motivations are the following:

Simplify the code in tof/unwrap.py (finding the tof origin with WFM was complex with a fake chopper and propagating frames backwards along the beamline)
Generalize the code so that we only have to use a single set of providers instead of having providers for standard unwrap, pulse-skipping and WFM (only the uncommon case where there are no choppers requires different providers).

Implementation

The implementation is basically described in https://tof.readthedocs.io/en/stable/wfm.html.

We use the chopper_cascade module to find the ranges of time and wavelength (polygons) that are reachable at the detector pixels.
We approximate the polygons with straight lines
We use these to define lookup tables which give us directly the tof as a function of arrival time at detector

…providers and types to unwrap

…ests in unwrap

jokasimr

Overall this looks more or less complete to me. Great job. (-300 SLOC is nice)
Left a couple of comments of things I think we should change before merging.

jokasimr · 2024-11-25T14:27:15Z

src/scippneutron/tof/unwrap.py

+        a /= weight.sum(dims)
+        # Compute the intercept as the mean of the intercepts of the curves that pass
+        # through the vertices
+        b = (y - a * x).mean('vertex')


I think we should use another method to determine the line that best approximates the polygon.
While this method works well in the cases we've seen so far it is not obvious that it will always produce something reasonable.

For example, if two vertices have almost the same x-coord the slope will grow a lot and the average slope might be larger than expected.

Yes I forgot to look at the notebook you sent.
Thanks for reminding me.

The suggested method indeed yields better results 👍

Old:

New:

jokasimr · 2024-11-25T14:29:12Z

src/scippneutron/tof/unwrap.py

-        # begin of the first bin, there may be a significant gap (which we could fill
-        # with NaNs or zeros), but there could also be a small overlap.
-    return UnwrappedData(da)
+    return TofCoord(toa)


Should we not assume the neutron left in the middle of the pulse? (Just because that makes a slightly better estimate.)
That is, shouldn't this be something like toa - sc.scalar(3, unit='ms')/2?

Yeah, good point, we could do that 👍

jokasimr · 2024-11-25T14:46:13Z

src/scippneutron/tof/unwrap.py

+    unsorted bin edges (due to either wrapping of `time_of_flight` or wavelength
+    overlap between subframes).
+    This function re-histograms the data to ensure that the bin edges are sorted.
+    It generates a number of events in each bin with a normal distribution around


It seems to me this function is essentially making a change of variables in a distribution function (or histogram).

That is, starting with a distribution over toa and a function tof(toa), we want to compute the distribution over tof.

To do this we should use the formula here: https://en.wikipedia.org/wiki/Probability_density_function#Function_of_random_variables_and_change_of_variables_in_the_probability_density_function (the not monotonic case)

If we used that I'm sure we'd end up with a simpler solution here.

If you could give me some pointers on how you would implement this, I'd be happy to try. I don't think I fully understood what's in the link.

Note also that this function is just to return something which is like the original data, but can be used further for normalization etc.
The real data is actually the TofData, but that has unsorted bin edges.
So maybe the details of this function are not super important because we can swap it out for a different/better provider later?
It was mostly designed so I could write some tests where I compare the wavelengths...

I was looking at this now and it became more complicated than I thought it would be 🤔
What you've already written might be the better solution.
Like you said, if we come up with something later it can be revisited.

… in tof

jokasimr · 2024-12-02T14:00:07Z

I'll take another look tomorrow!

…dling of variances in rehistogrammed data

jokasimr · 2024-12-05T11:28:58Z

I'm still working on the review, it'll be done later today

jokasimr

I think there might be a wrong assumption here that the minimum event_time_zero in the file is also the start of a stride in pulse skipping mode.

But what might happen is that the user starts writing the file in the middle of the stride when pulse skipping is used.

Is that taken into consideration here?

Note that this is probably not an issue for simulated data.

jokasimr · 2024-12-05T07:59:01Z

src/scippneutron/tof/fakes.py

+            # the boundaries have crossed, which means that the chopper is rotating
+            # must rotate at a lower frequency.


Suggested change

# the boundaries have crossed, which means that the chopper is rotating

# must rotate at a lower frequency.

# the boundaries have crossed, which means that the chopper

# is rotating at a lower frequency.

jokasimr · 2024-12-05T08:33:48Z

src/scippneutron/tof/unwrap.py

    """
    if da.bins is None:
        # Canonical name in NXmonitor
-        return PulseWrappedTimeOffset(da.coords['time_of_flight'])
-    return PulseWrappedTimeOffset(da.bins.coords['event_time_offset'])
+        toa = da.coords['time_of_flight']


Is the time_of_flight field already unwrapped? That is, does it denote true time of flight, or does it just denote time of arrival but has a misleading name?

The latter. It's just the name that is given in NeXus files.

jokasimr · 2024-12-05T08:54:38Z

src/scippneutron/tof/unwrap.py

+        # Hence we use the smallest event_time_zero as the time origin.
+        time_zero = da.coords['event_time_zero'] - da.coords['event_time_zero'].min()
+        coord = da.bins.coords['event_time_offset']
+        toa = coord + time_zero.to(dtype=float, unit=elem_unit(coord), copy=False)


I think there might be something wrong here in the case when pulse skipping is used.

We know that

toa = event_time_offset + event_time_zero = event_time_offset + frame_time_zero + frame_time_offset

where frame_time_zero is the time the start of the latest "logical pulse" (the start of the most recent pulse that wasn't blocked by the pulse-skip chopper) and frame_time_offset = event_time_zero - frame_time_zero.
From the chopper cascade we know that the frame takes frame_travel_time to reach the detector.
Therefore the frame starts at start_time = frame_travel_time + frame_time_zero.

Then

toa_minus_start_time_modulo_frame_period = (toa - start_time) % frame_period = (event_time_offset + event_time_zero - start_time) % frame_period = (event_time_offset + frame_time_offset - frame_travel_time) % frame_period

so I'd suggest that we

Remove UnwrappedTimeOfArrival since (it seems) it's only used to compute UnwrappedTimeOfArrivalMinusStartTime.

Rename FrameAtDetectorStartTime to FrameTravelTime.

Add a provider that computes FrameTimeOffset - this might be complicated. Typically I guess we have to get this from the chopper configuration somehow. For simulated data this is easy because then frame_time_offset = (event_time_zero - minimum_event_time_zero) % frame_period just like in the code here.

I fixed the case where a pulse skipping chopper would have a phase offset of 180 deg by propagating multiple pulses through the chopper cascade to predict frame at detector location.

I fixed the potential issue of starting to record in the middle of a pulse-skipping frame by adding a PulseStrideOffset integer parameter, which can be set by the user if the data doesn't look right, and which is set to 0 by default.

jokasimr · 2024-12-05T12:20:35Z

src/scippneutron/tof/unwrap.py

-    time_zero = 0.5 * (source_time_open + source_time_close)
-    return TimeOfFlightOrigin(time=time_zero, distance=source_chopper.distance)
+    return TimeOfArrivalModuloPeriod(
+        toa % frame_period.to(unit=elem_unit(toa), copy=False)


Related to the above comment. What we want here is the time since the most recent frame_time_zero.

That is event_time_offset + frame_time_offset.

In a simulation we know that the first event_time_zero coincides with the first frame_time_zero and then what is implemented here is correct.

…0deg offset

…e stride we should start with when pulse skipping

… frame is missing from raw data

jokasimr · 2024-12-11T09:00:04Z

Looks good to me!

nvaytet and others added 19 commits November 9, 2024 23:50

update disk_chopper and essbeamline fake from other PR. Start adding …

aa5185b

…providers and types to unwrap

complete provider setup for alternate tof workflow

bfd785e

fix implementation to handle multi-pixels

4e682fd

start purging old workflow

23c896f

use FramePeriod and start looking at the tests

83a02be

Merge branch 'main' into new-wfm

de66260

fix wfm tests

65c6370

update some tests

d3f3d0a

add providers for the case where there are no choppers

fb203ef

fix for histogrammed input data

19097f4

order of dimension in operation

e1801c5

add re-histogramming provider

cf1741d

dimension fixes in rehistogram and fix test

22d7b4e

use ess fake in pulse skipping tests

bb14f32

fix pulse skipping tests

c518155

fix dream wfm test and remove old to_time_of_flight test covered by t…

75869db

…ests in unwrap

cleanup unwrap cose

e710488

cleanup

ecd6988

update notebooks in docs

cc1fbd3

nvaytet requested review from SimonHeybrock and jokasimr November 25, 2024 13:32

jokasimr reviewed Nov 25, 2024

View reviewed changes

nvaytet and others added 8 commits November 26, 2024 14:25

use polygon least-squares method which yields better results

ad490f0

use half pulse offset when we have no choppers

9849337

add event data generator to docs module. This should probably just go…

140a30a

… in tof

update frame unwrapping notebook

ff3b8e7

update tof version

15618f7

Apply automatic formatting

a8f9df5

formatting

094e9bc

Merge branch 'main' into new-wfm

e785665

nvaytet marked this pull request as ready for review December 2, 2024 13:52

nvaytet added 4 commits December 4, 2024 11:52

use nanmedian for new bin width to remove outliers, and add basic han…

97e0336

…dling of variances in rehistogrammed data

move polygon approximation code into a separate function

d0a76c6

add comments

6602bcf

use uniform distribution instead of normal to generate events in bins

85f9384

SimonHeybrock approved these changes Dec 5, 2024

View reviewed changes

jokasimr reviewed Dec 5, 2024

View reviewed changes

nvaytet and others added 16 commits December 9, 2024 13:53

fix frame at detector in the case where pulse skipping chopper has 18…

e288a5c

…0deg offset

add test with 180deg offset with pulse skipping

8a1c0fb

Apply automatic formatting

e4925c8

start adding PulseStrideOffset so we can set which pulse number in th…

1ca97b1

…e stride we should start with when pulse skipping

add default workflow params and start adding test where first half of…

7a801f8

… frame is missing from raw data

Apply automatic formatting

8e91f95

fix test with half frame missing

4b51d09

Merge branch 'new-wfm' of github.com:scipp/scippneutron into new-wfm

5de689b

add to_events function and use that in re-histogram function

accc51e

fix tof tests

6a48eab

update notebook on frame unwrapping

202a542

Merge branch 'main' into new-wfm

24ca971

update wfm notebooks

b9bed49

Apply automatic formatting

e29b9e7

fix dream notebook

9956fcd

Merge branch 'new-wfm' of github.com:scipp/scippneutron into new-wfm

4af881e

jokasimr approved these changes Dec 11, 2024

View reviewed changes

add support for masks in to_events

c3823a6

nvaytet merged commit 2fc22c5 into main Dec 11, 2024
5 checks passed

nvaytet deleted the new-wfm branch December 11, 2024 10:46

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Alternative to time-origin tof computation #577

Alternative to time-origin tof computation #577

nvaytet commented Nov 25, 2024

jokasimr left a comment •

edited

Loading

jokasimr Nov 25, 2024 •

edited

Loading

nvaytet Nov 25, 2024

nvaytet Nov 26, 2024

jokasimr Nov 25, 2024

nvaytet Nov 26, 2024

jokasimr Nov 25, 2024 •

edited

Loading

nvaytet Nov 25, 2024

nvaytet Nov 26, 2024 •

edited

Loading

jokasimr Nov 26, 2024

jokasimr commented Dec 2, 2024

jokasimr commented Dec 5, 2024

jokasimr left a comment

jokasimr Dec 5, 2024

jokasimr Dec 5, 2024

nvaytet Dec 11, 2024

jokasimr Dec 5, 2024 •

edited

Loading

nvaytet Dec 10, 2024

jokasimr Dec 5, 2024 •

edited

Loading

jokasimr commented Dec 11, 2024

		# the boundaries have crossed, which means that the chopper is rotating
		# must rotate at a lower frequency.

Alternative to time-origin tof computation #577

Alternative to time-origin tof computation #577

Conversation

nvaytet commented Nov 25, 2024

jokasimr left a comment • edited Loading

Choose a reason for hiding this comment

jokasimr Nov 25, 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

jokasimr Nov 25, 2024 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

nvaytet Nov 26, 2024 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

jokasimr commented Dec 2, 2024

jokasimr commented Dec 5, 2024

jokasimr 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

jokasimr Dec 5, 2024 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

jokasimr Dec 5, 2024 • edited Loading

Choose a reason for hiding this comment

jokasimr commented Dec 11, 2024

jokasimr left a comment •

edited

Loading

jokasimr Nov 25, 2024 •

edited

Loading

jokasimr Nov 25, 2024 •

edited

Loading

nvaytet Nov 26, 2024 •

edited

Loading

jokasimr Dec 5, 2024 •

edited

Loading

jokasimr Dec 5, 2024 •

edited

Loading