feat: frame unwrapping for estia #185

jokasimr · 2025-11-27T12:34:00Z

Fixes #141

~~The approach requires five pieces of information:~~

1. The chopper position.
2. The chopper opening time.
3. The chopper frequency (only relevant for pulse skipping).
4. The average time at the source. (This might be fine to hard code).
5. The maximum time at the source. (This might be fine to hard code).

We're going to use the Tof lookup table approach from essreduce.
For that we need the chopper information required to run a simulation in tof.

jokasimr · 2025-11-27T12:36:12Z

src/ess/estia/conversions.py

+def mcstas_wavelength_coordinate_transformation_graph() -> CoordTransformationGraph:
+    return {
+        "wavelength": lambda wavelength_from_mcstas: wavelength_from_mcstas,
+        "theta": theta,
+        "divergence_angle": divergence_angle,
+        "Q": reflectometry_q,
+        "L1": lambda source_position, sample_position: sc.norm(
+            sample_position - source_position
+        ),  # + extra correction for guides?
+        "L2": lambda position, sample_position: sc.norm(position - sample_position),
+    }


This is added so that we can still use the true wavelengths from McStas in cases when we don't want the wavelength uncertainty to contribute to the resolution.

jokasimr · 2025-11-27T12:37:45Z

src/ess/estia/load.py

+
+    da.coords['chopper_distance'] = sc.scalar(10.895, unit='m')
+    velocity_max = sc.scalar(3.956034e3 / 3.75, unit='m/s')
+    da.coords['chopper_open_time'] = da.coords['chopper_distance'] / velocity_max


chopper_open_time and chopper_distance should be read from the file, but in the Estia McStas files the data is not stored in the DiskChopper component.

Is this info still required here or is it now encoded in the lookup table?

It's not required and will be removed 👍

nvaytet · 2025-11-27T14:23:23Z

src/ess/estia/conversions.py

-    # Use frame unwrapping from scippneutron
-    raise NotImplementedError()
+    """
+    Converts event_time_offset to time_of_flight.


I think we should create a lookup table and use the GenericTofWorkflow instead?

We need the chopper parameters to make a tof model for estia, or alternatively if there is a good mcstas model, then use that?

I'm thinking that it might be overkill to use the lookup table approach here.

Estia has a single chopper, and as I understand it that chopper only eliminates frame overlap, the wavelength resolution is the natural resolution of the source.

This is approximately what it looks like:

You can still get a better estimate from making a lookup table. You get a much more uniform agreement to wavelengths. With a single time/distance origin, the results are good at the center of the frames, and get worse on the edges of the frames.

I don't think it's overkill and it would mean it's using a workflow which is common to other instruments, which is also good? Especially if there is only one chopper, making a table will be very easy (not so many configurations) 🙂

Okay that makes sense, I'll take a stab at it tomorrow 👍

See here for an example of how to make a table (Dream instrument):
https://scipp.github.io/essdiffraction/user-guide/dream/dream-make-tof-lookup-table.html

nvaytet

Looking good 👍

Can you also update the PR description to reflect the update?

nvaytet · 2025-12-01T16:41:44Z

docs/user-guide/estia/create-estia-tof-lookup-table.ipynb

+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "table.plot()"


The table has a big gap in the middle, probably because of overlap between frames. Is that expected or has something gone wrong with the chopper parameters? (in the tof time-distance diagram you posted in the PR, there didn't seem to be frame overlap?)

The time distance diagram I posted didn't use exactly correct chopper parameters, it was just meant illustrate the chopper system qualitatively.

Yeah, there's a big overlap in the lookup table ,and it is a bit problematic because we loose quite a lot of intensity there.
There is an overlap in the McStas data too, but is smaller than the overlap in the lookup table:

So something seems to be wrong.

Aha, found the issue. The LookupTableRelativeErrorThreshold was set too low.

If I don't misremember a resonable wavelength uncertainty for Estia is ~6%, is that right @paracini?
With the relative error threshold set to 6% the lookup table looks like below:

The time distance diagram I posted didn't use exactly correct chopper parameters, it was just meant illustrate the chopper system qualitatively.

I guess it was an ok representation because in the table above, around 35m (where the detector was placed in the other diagram), there is almost no overlap. You can guess from the time-distance diagram that going to higher distances would lead to more overlap.

yes resolution is between 8% for short wavelengths and 3.5% for the long wavelengths

nvaytet · 2025-12-01T16:47:52Z

docs/user-guide/estia/estia-advanced-mcstas-reduction.ipynb

+    "#from ess.estia.conversions import mcstas_wavelength_coordinate_transformation_graph\n",
+    "#wf.insert( mcstas_wavelength_coordinate_transformation_graph )\n",
+    "\n",
+    "wf[TimeOfFlightLookupTable] = estia_tof_lookup_table()\n",


You should set the TimeOfFlightLookupTableFilename on the workflow instead of the table itself, and let the workflow do the loading (there is some handling between old and new file formats that you would be missing here).

See scipp/essspectroscopy#80

nvaytet · 2025-12-01T16:59:03Z

docs/user-guide/estia/estia-advanced-mcstas-reduction.ipynb

    "        sc.scalar(float('nan'), unit=R.unit),\n",
-    "        R.data\n",
+    "        R.data,\n",
    "    )\n",


It's kind of weird. I was looking at how the results changed. To me, the new results look more like other reflectometry curves I've seen before (the old ones look very 'spiky').

But when I then look at the plot further down with the 'ground truth' in black, it seems the old results are much better? How was the ground truth result computed?

This is expected, the old results used the exact wavelengths from McStas, so the resolution has no contribution from wavelength uncertainty. When we do frame unwrapping then wavelength uncertainty contributes to the resolution (or lack thereof). As I understand it wavelength uncertainty is the dominant factor of the resolution of Estia.

The black curve shows the ground truth reflectivity of the sample. When the resolution is good the computed reflectivity overlaps it. But with frame unwrapping the sharp peaks are smeared out.

That's correct, resolution on Estia is dominated by wavelength, except at very small angles (around the critical edge) where the detector resolution dominates

nvaytet · 2025-12-01T17:01:23Z

docs/user-guide/estia/estia-advanced-mcstas-reduction.ipynb

+    "#from ess.estia.conversions import mcstas_wavelength_coordinate_transformation_graph\n",
+    "#wf.insert( mcstas_wavelength_coordinate_transformation_graph )\n",
+    "\n",
+    "wf[TimeOfFlightLookupTable] = estia_tof_lookup_table()\n",


Also set the filename

nvaytet · 2025-12-01T17:02:20Z

docs/user-guide/estia/estia-advanced-mcstas-reduction.ipynb

+    "wf[TimeOfFlightLookupTable] = estia_tof_lookup_table()\n",
    "\n",
    "# There is no proton current data in the McStas files, here we just add some fake proton current\n",
    "# data to make the workflow run.\n",


The new diagnostics figure looks weird, especially the top left subplot

nvaytet · 2025-12-01T17:05:03Z

docs/user-guide/estia/simulated-spin-flip-sample.ipynb

    "\n",
+    "\n",
+    "wf.insert( mcstas_wavelength_coordinate_transformation_graph )\n",
+    "wf[TimeOfFlightLookupTable] = estia_tof_lookup_table()\n",


Use TimeOfFlightLookupTableFilename

nvaytet · 2025-12-01T20:01:45Z

src/ess/estia/data.py



+def estia_tof_lookup_table(pulse_stride=1):
+    if pulse_stride == 1:


We should use the loader in essreduce instead.

About the pulse_stride=1 logic, this will hopefully just be handled by the future mechanism in the generic workflow that would check that the chopper params in the file that is being loaded are compatible with the lookup table being used. So it should not be needed here.

nvaytet · 2025-12-01T20:03:36Z

src/ess/estia/beamline.py



-def LoadNeXusWorkflow() -> sciline.Pipeline:
+def LoadNeXusWorkflow(**kwargs) -> sciline.Pipeline:


Maybe we should pick a different name for the workflow (something like EstiaWorkflow?) as it is no longer just loading from nexus?

Hm, but we do already define EstiaWorkflow in src/ess/estia/workflow.py.

nvaytet · 2025-12-01T20:04:27Z

tests/estia/mcstas_data_test.py

    wf[WavelengthBins] = sc.geomspace('wavelength', 3.5, 12, 2001, unit='angstrom')
    wf[QBins] = sc.geomspace('Q', 0.005, 0.1, 200, unit='1/angstrom')
+
+    wf[TimeOfFlightLookupTable] = estia_tof_lookup_table()


Use TimeOfFlightLookupTableFilename

nvaytet · 2025-12-01T20:05:50Z

src/ess/estia/load.py

+
+    da.coords['chopper_distance'] = sc.scalar(10.895, unit='m')
+    velocity_max = sc.scalar(3.956034e3 / 3.75, unit='m/s')
+    da.coords['chopper_open_time'] = da.coords['chopper_distance'] / velocity_max


Is this info still required here or is it now encoded in the lookup table?

jokasimr · 2025-12-02T08:38:27Z

Is this info still required here or is it now encoded in the lookup table?

It is encoded in the lookup table, so that code can be removed.

docs/user-guide/estia/simulated-spin-flip-sample.ipynb

nvaytet

I'm happy with the changes.

We need to make sure @paracini is also happy, as the results have changed significantly.

jokasimr · 2025-12-02T10:26:14Z

We need to make sure @paracini is also happy, as the results have changed significantly.

Note that you can still get back the old results by changing the coordinate transformations used to compute wavelength.
Specifically, if you add this provider to the workflow:

from ess.estia.conversions import mcstas_wavelength_coordinate_transformation_graph
wf.insert( mcstas_wavelength_coordinate_transformation_graph )

then wavelength will be "computed" from the exact wavelength coordinate that is loaded from McStas, like was done before.

feat: frame unwrapping for estia

6c10cd2

jokasimr requested review from nvaytet and paracini November 27, 2025 12:34

jokasimr commented Nov 27, 2025

View reviewed changes

nvaytet reviewed Nov 27, 2025

View reviewed changes

jokasimr added 8 commits November 28, 2025 17:25

feat: use generic tof workflow with lookuptable

8ef33cb

fix: use definion from scippneutron

4027d60

docs: add notebook for computing tof table

f161ba7

docs

0738aac

docs

5fbd792

docs

87ad30a

deps: add tof to docs deps

01da0da

docs: update spin-flip-example

dd4592a

jokasimr requested a review from nvaytet December 1, 2025 10:57

jokasimr added 4 commits December 1, 2025 13:47

docs: add lookup table notebook to index

548454a

fix

c1f9fd2

docs: add title and short description to tof table notebook

fec0d60

docs: remove duplicate link to tof notebook

dfe5bb2

nvaytet reviewed Dec 1, 2025

View reviewed changes

fix: change error threshold and update table

29056fb

jokasimr added 4 commits December 2, 2025 09:41

docs: cleanup

d44338f

fix: unnecessary when we use tof lut

c22760e

fix

579dbfb

fix: set lookup table filename

30cb482

jokasimr requested a review from nvaytet December 2, 2025 09:45

nvaytet reviewed Dec 2, 2025

View reviewed changes

docs/user-guide/estia/simulated-spin-flip-sample.ipynb Outdated Show resolved Hide resolved

nvaytet approved these changes Dec 2, 2025

View reviewed changes

docs

34e79d4

fix: change number of bins to make figure less noisy

26e06c0

jokasimr mentioned this pull request Dec 4, 2025

Mcstas to nexus #186

Merged

paracini merged commit 7277ed8 into main Dec 4, 2025
4 checks passed

paracini deleted the estia-tof branch December 4, 2025 13:28



		def estia_tof_lookup_table(pulse_stride=1):
		if pulse_stride == 1:



		def LoadNeXusWorkflow() -> sciline.Pipeline:
		def LoadNeXusWorkflow(**kwargs) -> sciline.Pipeline:

feat: frame unwrapping for estia #185

feat: frame unwrapping for estia #185

Uh oh!

Conversation

jokasimr commented Nov 27, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

jokasimr Nov 27, 2025 • 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!

jokasimr Nov 27, 2025 • 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!

nvaytet Nov 27, 2025 • 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!

nvaytet Nov 28, 2025 • 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!

nvaytet left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

jokasimr Dec 2, 2025 • 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!

nvaytet Dec 2, 2025 • 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!

jokasimr Dec 2, 2025 • 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!

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 commented Nov 27, 2025 •

edited

Loading

jokasimr Nov 27, 2025 •

edited

Loading

jokasimr Nov 27, 2025 •

edited

Loading

nvaytet Nov 27, 2025 •

edited

Loading

nvaytet Nov 28, 2025 •

edited

Loading

jokasimr Dec 2, 2025 •

edited

Loading

nvaytet Dec 2, 2025 •

edited

Loading

jokasimr Dec 2, 2025 •

edited

Loading