Wavelength frame multiplication

Wavelength frame multiplication (WFM) is a technique commonly used at long-pulse facilities to improve the resolution of the results measured at the neutron detectors. See for example the article by Schmakat et al. (2020) for a description of how WFM works.

In this notebook, we show how to use essreduce’s time_of_flight module to compute an accurate a time-of-flight coordinate, from which a wavelength can be computed.

import numpy as np
import plopp as pp
import scipp as sc
import scippnexus as snx
from scippneutron.chopper import DiskChopper
from ess.reduce.nexus.types import AnyRun, RawDetector, SampleRun, NeXusDetectorName
from ess.reduce.time_of_flight import *

Setting up the beamline

Creating the beamline choppers

We begin by defining the chopper settings for our beamline. In principle, the chopper setting could simply be read from a NeXus file.

For this example, we create choppers modeled on the V20 ESS beamline at HZB. It consists of 5 choppers:

• 2 WFM choppers

• 2 frame-overlap choppers

• 1 pulse-overlap chopper

The first 4 choppers have 6 openings (also known as cutouts), while the last one only has a single opening.

wfm1 = DiskChopper(
    frequency=sc.scalar(-70.0, unit="Hz"),
    beam_position=sc.scalar(0.0, unit="deg"),
    phase=sc.scalar(-47.10, unit="deg"),
    axle_position=sc.vector(value=[0, 0, 6.6], unit="m"),
    slit_begin=sc.array(
        dims=["cutout"],
        values=np.array([83.71, 140.49, 193.26, 242.32, 287.91, 330.3]) + 15.0,
        unit="deg",
    ),
    slit_end=sc.array(
        dims=["cutout"],
        values=np.array([94.7, 155.79, 212.56, 265.33, 314.37, 360.0]) + 15.0,
        unit="deg",
    ),
    slit_height=sc.scalar(10.0, unit="cm"),
    radius=sc.scalar(30.0, unit="cm"),
)

wfm2 = DiskChopper(
    frequency=sc.scalar(-70.0, unit="Hz"),
    beam_position=sc.scalar(0.0, unit="deg"),
    phase=sc.scalar(-76.76, unit="deg"),
    axle_position=sc.vector(value=[0, 0, 7.1], unit="m"),
    slit_begin=sc.array(
        dims=["cutout"],
        values=np.array([65.04, 126.1, 182.88, 235.67, 284.73, 330.32]) + 15.0,
        unit="deg",
    ),
    slit_end=sc.array(
        dims=["cutout"],
        values=np.array([76.03, 141.4, 202.18, 254.97, 307.74, 360.0]) + 15.0,
        unit="deg",
    ),
    slit_height=sc.scalar(10.0, unit="cm"),
    radius=sc.scalar(30.0, unit="cm"),
)

foc1 = DiskChopper(
    frequency=sc.scalar(-56.0, unit="Hz"),
    beam_position=sc.scalar(0.0, unit="deg"),
    phase=sc.scalar(-62.40, unit="deg"),
    axle_position=sc.vector(value=[0, 0, 8.8], unit="m"),
    slit_begin=sc.array(
        dims=["cutout"],
        values=np.array([74.6, 139.6, 194.3, 245.3, 294.8, 347.2]),
        unit="deg",
    ),
    slit_end=sc.array(
        dims=["cutout"],
        values=np.array([95.2, 162.8, 216.1, 263.1, 310.5, 371.6]),
        unit="deg",
    ),
    slit_height=sc.scalar(10.0, unit="cm"),
    radius=sc.scalar(30.0, unit="cm"),
)

foc2 = DiskChopper(
    frequency=sc.scalar(-28.0, unit="Hz"),
    beam_position=sc.scalar(0.0, unit="deg"),
    phase=sc.scalar(-12.27, unit="deg"),
    axle_position=sc.vector(value=[0, 0, 15.9], unit="m"),
    slit_begin=sc.array(
        dims=["cutout"],
        values=np.array([98.0, 154.0, 206.8, 255.0, 299.0, 344.65]),
        unit="deg",
    ),
    slit_end=sc.array(
        dims=["cutout"],
        values=np.array([134.6, 190.06, 237.01, 280.88, 323.56, 373.76]),
        unit="deg",
    ),
    slit_height=sc.scalar(10.0, unit="cm"),
    radius=sc.scalar(30.0, unit="cm"),
)

pol = DiskChopper(
    frequency=sc.scalar(-14.0, unit="Hz"),
    beam_position=sc.scalar(0.0, unit="deg"),
    phase=sc.scalar(0.0, unit="deg"),
    axle_position=sc.vector(value=[0, 0, 17.0], unit="m"),
    slit_begin=sc.array(
        dims=["cutout"],
        values=np.array([40.0]),
        unit="deg",
    ),
    slit_end=sc.array(
        dims=["cutout"],
        values=np.array([240.0]),
        unit="deg",
    ),
    slit_height=sc.scalar(10.0, unit="cm"),
    radius=sc.scalar(30.0, unit="cm"),
)

disk_choppers = {"wfm1": wfm1, "wfm2": wfm2, "foc1": foc1, "foc2": foc2, "pol": pol}

It is possible to visualize the properties of the choppers by inspecting their repr:

wfm1

DiskChopper

(cutout: 6)

• axle_position

  scipp

  Variable

  ()

  vector3

  m

  [0. 0. 6.6]
• frequency

  scipp

  Variable

  ()

  float64

  Hz

  -70.0
• beam_position

  scipp

  Variable

  ()

  float64

  deg

  0.0
• phase

  scipp

  Variable

  ()

  float64

  deg

  -47.1
• slit_begin

  scipp

  Variable

  (cutout: 6)

  float64

  deg

  98.71, 155.490, ..., 302.910, 345.3
• slit_end

  scipp

  Variable

  (cutout: 6)

  float64

  deg

  109.7, 170.790, ..., 329.370, 375.0
• slit_height

  scipp

  Variable

  (cutout: 6)

  float64

  cm

  10.0, 10.0, ..., 10.0, 10.0
• radius

  scipp

  Variable

  ()

  float64

  cm

  30.0

Define the source position which is required to compute the distance that neutrons travelled. In this example, chopper positions are given relative to the source, so we set the source position to the origin.

source_position = sc.vector([0, 0, 0], unit="m")

Adding a detector

We also have a detector, which we place 26 meters away from the source.

Ltotal = sc.scalar(26.0, unit="m")

Creating some neutron events

We create a semi-realistic set of neutron events based on the ESS pulse.

from ess.reduce.time_of_flight.fakes import FakeBeamline

ess_beamline = FakeBeamline(
    choppers=disk_choppers,
    source_position=source_position,
    monitors={"detector": Ltotal},
    run_length=sc.scalar(1 / 14, unit="s") * 14,
    events_per_pulse=200_000,
)

The initial birth times and wavelengths of the generated neutrons can be visualized (for a single pulse):

one_pulse = ess_beamline.source.data["pulse", 0]
one_pulse.hist(birth_time=300).plot() + one_pulse.hist(wavelength=300).plot()

From this fake beamline, we extract the raw neutron signal at our detector:

raw_data = ess_beamline.get_monitor("detector")[0]

# Visualize
raw_data.hist(event_time_offset=300).squeeze().plot()

The total number of neutrons in our sample data that make it through the to detector is:

raw_data.sum().value

np.float64(204163.0)

Computing time-of-flight

Next, we use a workflow that provides an estimate of the real time-of-flight as a function of neutron time-of-arrival.

Setting up the workflow

wf = GenericTofWorkflow(run_types=[SampleRun], monitor_types=[])

wf[RawDetector[SampleRun]] = raw_data
wf[DetectorLtotal[SampleRun]] = Ltotal
wf[NeXusDetectorName] = "detector"

wf.visualize(TofDetector[SampleRun])

---------------------------------------------------------------------------
FileNotFoundError                         Traceback (most recent call last)
File ~/work/ess/ess/.pixi/envs/docs-essreduce/lib/python3.11/site-packages/graphviz/backend/execute.py:76, in run_check(cmd, input_lines, encoding, quiet, **kwargs)
     75         kwargs['stdout'] = kwargs['stderr'] = subprocess.PIPE
---> 76     proc = _run_input_lines(cmd, input_lines, kwargs=kwargs)
     77 else:

File ~/work/ess/ess/.pixi/envs/docs-essreduce/lib/python3.11/site-packages/graphviz/backend/execute.py:96, in _run_input_lines(cmd, input_lines, kwargs)
     95 def _run_input_lines(cmd, input_lines, *, kwargs):
---> 96     popen = subprocess.Popen(cmd, stdin=subprocess.PIPE, **kwargs)
     98     stdin_write = popen.stdin.write

File ~/work/ess/ess/.pixi/envs/docs-essreduce/lib/python3.11/subprocess.py:1026, in Popen.__init__(self, args, bufsize, executable, stdin, stdout, stderr, preexec_fn, close_fds, shell, cwd, env, universal_newlines, startupinfo, creationflags, restore_signals, start_new_session, pass_fds, user, group, extra_groups, encoding, errors, text, umask, pipesize, process_group)
   1023             self.stderr = io.TextIOWrapper(self.stderr,
   1024                     encoding=encoding, errors=errors)
-> 1026     self._execute_child(args, executable, preexec_fn, close_fds,
   1027                         pass_fds, cwd, env,
   1028                         startupinfo, creationflags, shell,
   1029                         p2cread, p2cwrite,
   1030                         c2pread, c2pwrite,
   1031                         errread, errwrite,
   1032                         restore_signals,
   1033                         gid, gids, uid, umask,
   1034                         start_new_session, process_group)
   1035 except:
   1036     # Cleanup if the child failed starting.

File ~/work/ess/ess/.pixi/envs/docs-essreduce/lib/python3.11/subprocess.py:1955, in Popen._execute_child(self, args, executable, preexec_fn, close_fds, pass_fds, cwd, env, startupinfo, creationflags, shell, p2cread, p2cwrite, c2pread, c2pwrite, errread, errwrite, restore_signals, gid, gids, uid, umask, start_new_session, process_group)
   1954 if err_filename is not None:
-> 1955     raise child_exception_type(errno_num, err_msg, err_filename)
   1956 else:

FileNotFoundError: [Errno 2] No such file or directory: PosixPath('dot')

The above exception was the direct cause of the following exception:

ExecutableNotFound                        Traceback (most recent call last)
File ~/work/ess/ess/.pixi/envs/docs-essreduce/lib/python3.11/site-packages/IPython/core/formatters.py:1036, in MimeBundleFormatter.__call__(self, obj, include, exclude)
   1033     method = get_real_method(obj, self.print_method)
   1035     if method is not None:
-> 1036         return method(include=include, exclude=exclude)
   1037     return None
   1038 else:

File ~/work/ess/ess/.pixi/envs/docs-essreduce/lib/python3.11/site-packages/graphviz/jupyter_integration.py:98, in JupyterIntegration._repr_mimebundle_(self, include, exclude, **_)
     96 include = set(include) if include is not None else {self._jupyter_mimetype}
     97 include -= set(exclude or [])
---> 98 return {mimetype: getattr(self, method_name)()
     99         for mimetype, method_name in MIME_TYPES.items()
    100         if mimetype in include}

File ~/work/ess/ess/.pixi/envs/docs-essreduce/lib/python3.11/site-packages/graphviz/jupyter_integration.py:98, in <dictcomp>(.0)
     96 include = set(include) if include is not None else {self._jupyter_mimetype}
     97 include -= set(exclude or [])
---> 98 return {mimetype: getattr(self, method_name)()
     99         for mimetype, method_name in MIME_TYPES.items()
    100         if mimetype in include}

File ~/work/ess/ess/.pixi/envs/docs-essreduce/lib/python3.11/site-packages/graphviz/jupyter_integration.py:112, in JupyterIntegration._repr_image_svg_xml(self)
    110 def _repr_image_svg_xml(self) -> str:
    111     """Return the rendered graph as SVG string."""
--> 112     return self.pipe(format='svg', encoding=SVG_ENCODING)

File ~/work/ess/ess/.pixi/envs/docs-essreduce/lib/python3.11/site-packages/graphviz/piping.py:104, in Pipe.pipe(self, format, renderer, formatter, neato_no_op, quiet, engine, encoding)
     55 def pipe(self,
     56          format: typing.Optional[str] = None,
     57          renderer: typing.Optional[str] = None,
   (...)     61          engine: typing.Optional[str] = None,
     62          encoding: typing.Optional[str] = None) -> typing.Union[bytes, str]:
     63     """Return the source piped through the Graphviz layout command.
     64
     65     Args:
   (...)    102         '<?xml version='
    103     """
--> 104     return self._pipe_legacy(format,
    105                              renderer=renderer,
    106                              formatter=formatter,
    107                              neato_no_op=neato_no_op,
    108                              quiet=quiet,
    109                              engine=engine,
    110                              encoding=encoding)

File ~/work/ess/ess/.pixi/envs/docs-essreduce/lib/python3.11/site-packages/graphviz/_tools.py:185, in deprecate_positional_args.<locals>.decorator.<locals>.wrapper(*args, **kwargs)
    177     wanted = ', '.join(f'{name}={value!r}'
    178                        for name, value in deprecated.items())
    179     warnings.warn(f'The signature of {func_name} will be reduced'
    180                   f' to {supported_number} positional arg{s_}{qualification}'
    181                   f' {list(supported)}: pass {wanted} as keyword arg{s_}',
    182                   stacklevel=stacklevel,
    183                   category=category)
--> 185 return func(*args, **kwargs)

File ~/work/ess/ess/.pixi/envs/docs-essreduce/lib/python3.11/site-packages/graphviz/piping.py:121, in Pipe._pipe_legacy(self, format, renderer, formatter, neato_no_op, quiet, engine, encoding)
    112 @_tools.deprecate_positional_args(supported_number=1, ignore_arg='self')
    113 def _pipe_legacy(self,
    114                  format: typing.Optional[str] = None,
   (...)    119                  engine: typing.Optional[str] = None,
    120                  encoding: typing.Optional[str] = None) -> typing.Union[bytes, str]:
--> 121     return self._pipe_future(format,
    122                              renderer=renderer,
    123                              formatter=formatter,
    124                              neato_no_op=neato_no_op,
    125                              quiet=quiet,
    126                              engine=engine,
    127                              encoding=encoding)

File ~/work/ess/ess/.pixi/envs/docs-essreduce/lib/python3.11/site-packages/graphviz/piping.py:149, in Pipe._pipe_future(self, format, renderer, formatter, neato_no_op, quiet, engine, encoding)
    146 if encoding is not None:
    147     if codecs.lookup(encoding) is codecs.lookup(self.encoding):
    148         # common case: both stdin and stdout need the same encoding
--> 149         return self._pipe_lines_string(*args, encoding=encoding, **kwargs)
    150     try:
    151         raw = self._pipe_lines(*args, input_encoding=self.encoding, **kwargs)

File ~/work/ess/ess/.pixi/envs/docs-essreduce/lib/python3.11/site-packages/graphviz/backend/piping.py:212, in pipe_lines_string(engine, format, input_lines, encoding, renderer, formatter, neato_no_op, quiet)
    206 cmd = dot_command.command(engine, format,
    207                           renderer=renderer,
    208                           formatter=formatter,
    209                           neato_no_op=neato_no_op)
    210 kwargs = {'input_lines': input_lines, 'encoding': encoding}
--> 212 proc = execute.run_check(cmd, capture_output=True, quiet=quiet, **kwargs)
    213 return proc.stdout

File ~/work/ess/ess/.pixi/envs/docs-essreduce/lib/python3.11/site-packages/graphviz/backend/execute.py:81, in run_check(cmd, input_lines, encoding, quiet, **kwargs)
     79 except OSError as e:
     80     if e.errno == errno.ENOENT:
---> 81         raise ExecutableNotFound(cmd) from e
     82     raise
     84 if not quiet and proc.stderr:

ExecutableNotFound: failed to execute PosixPath('dot'), make sure the Graphviz executables are on your systems' PATH

<graphviz.graphs.Digraph at 0x7ff2969f2410>

By default, the workflow tries to load a TofLookupTable from a file.

In this notebook, instead of using such a pre-made file, we will build our own lookup table from the chopper information and apply it to the workflow.

Building the time-of-flight lookup table

We use the Tof package to propagate a pulse of neutrons through the chopper system to the detectors, and predict the most likely neutron wavelength for a given time-of-arrival and distance from source.

From this, we build a lookup table on which bilinear interpolation is used to compute a wavelength (and its corresponding time-of-flight) for every neutron event.

lut_wf = TofLookupTableWorkflow()
lut_wf[DiskChoppers[AnyRun]] = disk_choppers
lut_wf[SourcePosition] = source_position
lut_wf[LtotalRange] = sc.scalar(5, unit='m'), sc.scalar(35, unit='m')
lut_wf.visualize(TofLookupTable)

---------------------------------------------------------------------------
FileNotFoundError                         Traceback (most recent call last)
File ~/work/ess/ess/.pixi/envs/docs-essreduce/lib/python3.11/site-packages/graphviz/backend/execute.py:76, in run_check(cmd, input_lines, encoding, quiet, **kwargs)
     75         kwargs['stdout'] = kwargs['stderr'] = subprocess.PIPE
---> 76     proc = _run_input_lines(cmd, input_lines, kwargs=kwargs)
     77 else:

File ~/work/ess/ess/.pixi/envs/docs-essreduce/lib/python3.11/site-packages/graphviz/backend/execute.py:96, in _run_input_lines(cmd, input_lines, kwargs)
     95 def _run_input_lines(cmd, input_lines, *, kwargs):
---> 96     popen = subprocess.Popen(cmd, stdin=subprocess.PIPE, **kwargs)
     98     stdin_write = popen.stdin.write

File ~/work/ess/ess/.pixi/envs/docs-essreduce/lib/python3.11/subprocess.py:1026, in Popen.__init__(self, args, bufsize, executable, stdin, stdout, stderr, preexec_fn, close_fds, shell, cwd, env, universal_newlines, startupinfo, creationflags, restore_signals, start_new_session, pass_fds, user, group, extra_groups, encoding, errors, text, umask, pipesize, process_group)
   1023             self.stderr = io.TextIOWrapper(self.stderr,
   1024                     encoding=encoding, errors=errors)
-> 1026     self._execute_child(args, executable, preexec_fn, close_fds,
   1027                         pass_fds, cwd, env,
   1028                         startupinfo, creationflags, shell,
   1029                         p2cread, p2cwrite,
   1030                         c2pread, c2pwrite,
   1031                         errread, errwrite,
   1032                         restore_signals,
   1033                         gid, gids, uid, umask,
   1034                         start_new_session, process_group)
   1035 except:
   1036     # Cleanup if the child failed starting.

File ~/work/ess/ess/.pixi/envs/docs-essreduce/lib/python3.11/subprocess.py:1955, in Popen._execute_child(self, args, executable, preexec_fn, close_fds, pass_fds, cwd, env, startupinfo, creationflags, shell, p2cread, p2cwrite, c2pread, c2pwrite, errread, errwrite, restore_signals, gid, gids, uid, umask, start_new_session, process_group)
   1954 if err_filename is not None:
-> 1955     raise child_exception_type(errno_num, err_msg, err_filename)
   1956 else:

FileNotFoundError: [Errno 2] No such file or directory: PosixPath('dot')

The above exception was the direct cause of the following exception:

ExecutableNotFound                        Traceback (most recent call last)
File ~/work/ess/ess/.pixi/envs/docs-essreduce/lib/python3.11/site-packages/IPython/core/formatters.py:1036, in MimeBundleFormatter.__call__(self, obj, include, exclude)
   1033     method = get_real_method(obj, self.print_method)
   1035     if method is not None:
-> 1036         return method(include=include, exclude=exclude)
   1037     return None
   1038 else:

File ~/work/ess/ess/.pixi/envs/docs-essreduce/lib/python3.11/site-packages/graphviz/jupyter_integration.py:98, in JupyterIntegration._repr_mimebundle_(self, include, exclude, **_)
     96 include = set(include) if include is not None else {self._jupyter_mimetype}
     97 include -= set(exclude or [])
---> 98 return {mimetype: getattr(self, method_name)()
     99         for mimetype, method_name in MIME_TYPES.items()
    100         if mimetype in include}

File ~/work/ess/ess/.pixi/envs/docs-essreduce/lib/python3.11/site-packages/graphviz/jupyter_integration.py:98, in <dictcomp>(.0)
     96 include = set(include) if include is not None else {self._jupyter_mimetype}
     97 include -= set(exclude or [])
---> 98 return {mimetype: getattr(self, method_name)()
     99         for mimetype, method_name in MIME_TYPES.items()
    100         if mimetype in include}

File ~/work/ess/ess/.pixi/envs/docs-essreduce/lib/python3.11/site-packages/graphviz/jupyter_integration.py:112, in JupyterIntegration._repr_image_svg_xml(self)
    110 def _repr_image_svg_xml(self) -> str:
    111     """Return the rendered graph as SVG string."""
--> 112     return self.pipe(format='svg', encoding=SVG_ENCODING)

File ~/work/ess/ess/.pixi/envs/docs-essreduce/lib/python3.11/site-packages/graphviz/piping.py:104, in Pipe.pipe(self, format, renderer, formatter, neato_no_op, quiet, engine, encoding)
     55 def pipe(self,
     56          format: typing.Optional[str] = None,
     57          renderer: typing.Optional[str] = None,
   (...)     61          engine: typing.Optional[str] = None,
     62          encoding: typing.Optional[str] = None) -> typing.Union[bytes, str]:
     63     """Return the source piped through the Graphviz layout command.
     64
     65     Args:
   (...)    102         '<?xml version='
    103     """
--> 104     return self._pipe_legacy(format,
    105                              renderer=renderer,
    106                              formatter=formatter,
    107                              neato_no_op=neato_no_op,
    108                              quiet=quiet,
    109                              engine=engine,
    110                              encoding=encoding)

File ~/work/ess/ess/.pixi/envs/docs-essreduce/lib/python3.11/site-packages/graphviz/_tools.py:185, in deprecate_positional_args.<locals>.decorator.<locals>.wrapper(*args, **kwargs)
    177     wanted = ', '.join(f'{name}={value!r}'
    178                        for name, value in deprecated.items())
    179     warnings.warn(f'The signature of {func_name} will be reduced'
    180                   f' to {supported_number} positional arg{s_}{qualification}'
    181                   f' {list(supported)}: pass {wanted} as keyword arg{s_}',
    182                   stacklevel=stacklevel,
    183                   category=category)
--> 185 return func(*args, **kwargs)

File ~/work/ess/ess/.pixi/envs/docs-essreduce/lib/python3.11/site-packages/graphviz/piping.py:121, in Pipe._pipe_legacy(self, format, renderer, formatter, neato_no_op, quiet, engine, encoding)
    112 @_tools.deprecate_positional_args(supported_number=1, ignore_arg='self')
    113 def _pipe_legacy(self,
    114                  format: typing.Optional[str] = None,
   (...)    119                  engine: typing.Optional[str] = None,
    120                  encoding: typing.Optional[str] = None) -> typing.Union[bytes, str]:
--> 121     return self._pipe_future(format,
    122                              renderer=renderer,
    123                              formatter=formatter,
    124                              neato_no_op=neato_no_op,
    125                              quiet=quiet,
    126                              engine=engine,
    127                              encoding=encoding)

File ~/work/ess/ess/.pixi/envs/docs-essreduce/lib/python3.11/site-packages/graphviz/piping.py:149, in Pipe._pipe_future(self, format, renderer, formatter, neato_no_op, quiet, engine, encoding)
    146 if encoding is not None:
    147     if codecs.lookup(encoding) is codecs.lookup(self.encoding):
    148         # common case: both stdin and stdout need the same encoding
--> 149         return self._pipe_lines_string(*args, encoding=encoding, **kwargs)
    150     try:
    151         raw = self._pipe_lines(*args, input_encoding=self.encoding, **kwargs)

File ~/work/ess/ess/.pixi/envs/docs-essreduce/lib/python3.11/site-packages/graphviz/backend/piping.py:212, in pipe_lines_string(engine, format, input_lines, encoding, renderer, formatter, neato_no_op, quiet)
    206 cmd = dot_command.command(engine, format,
    207                           renderer=renderer,
    208                           formatter=formatter,
    209                           neato_no_op=neato_no_op)
    210 kwargs = {'input_lines': input_lines, 'encoding': encoding}
--> 212 proc = execute.run_check(cmd, capture_output=True, quiet=quiet, **kwargs)
    213 return proc.stdout

File ~/work/ess/ess/.pixi/envs/docs-essreduce/lib/python3.11/site-packages/graphviz/backend/execute.py:81, in run_check(cmd, input_lines, encoding, quiet, **kwargs)
     79 except OSError as e:
     80     if e.errno == errno.ENOENT:
---> 81         raise ExecutableNotFound(cmd) from e
     82     raise
     84 if not quiet and proc.stderr:

ExecutableNotFound: failed to execute PosixPath('dot'), make sure the Graphviz executables are on your systems' PATH

<graphviz.graphs.Digraph at 0x7ff296062490>

Inspecting the lookup table

The workflow first runs a simulation using the chopper parameters above, and the result is stored in SimulationResults (see graph above).

From these simulated neutrons, we create figures displaying the neutron wavelengths and time-of-flight, as a function of arrival time at the detector.

This is the basis for creating our lookup table.

sim = lut_wf.compute(SimulationResults)


def to_event_time_offset(sim):
    # Compute event_time_offset at the detector
    eto = (
        sim.time_of_arrival + ((Ltotal - sim.distance) / sim.speed).to(unit="us")
    ) % sc.scalar(1e6 / 14.0, unit="us")
    # Compute time-of-flight at the detector
    tof = (Ltotal / sim.speed).to(unit="us")
    return sc.DataArray(
        data=sim.weight,
        coords={"wavelength": sim.wavelength, "event_time_offset": eto, "tof": tof},
    )


events = to_event_time_offset(sim.readings["pol"])
fig1 = events.hist(wavelength=300, event_time_offset=300).plot(norm="log")
fig2 = events.hist(tof=300, event_time_offset=300).plot(norm="log")
fig1 + fig2

The lookup table is then obtained by computing the weighted mean of the time-of-flight inside each time-of-arrival bin.

This is illustrated by the orange line in the figure below:

table = lut_wf.compute(TofLookupTable)

# Overlay mean on the figure above
table.array["distance", 212].plot(ax=fig2.ax, color="C1", ls="-", marker=None)

# Zoom in
fig2.canvas.xrange = 40000, 50000
fig2.canvas.yrange = 35000, 50000
fig2

We can see that the orange lines follow the center of the colored areas.

We can also see that in regions where there is contamination from other chopper openings (overlapping regions in time), the error bars on the orange line get larger.

Another way of looking at this is to plot the entire table, showing the predicted time-of-flight as a function of event_time_offset and distance. We also show alongside it the standard deviation of the predicted time-of-flight.

def plot_lut(table):
    fig = table.plot(title="Predicted time-of-flight") + sc.stddevs(table.array).plot(title="Standard deviation", vmax=1000)
    for f in (fig[0, 0], fig[0, 1]):
        f.ax.axhline(Ltotal.value, ls='dashed', color='k')
        f.ax.text(1e3, Ltotal.value, "detector", va='bottom', color='k')
        f.ax.text(1e3, Ltotal.value, "at 26m", va='top', color='k')
    return fig

plot_lut(table)

We can see that at low distances (< 7m), before the first choppers, the uncertainties are very large. Neutrons with different wavelengths, originating from different parts of the pulse, are mixing and make it very difficult to predict a good time-of-flight.

As we move to larger distances, uncertainties drop overall (colors drift to blue). However, we still see spikes of uncertainties inside the frames around 26 m where the detector is placed, which indicates contamination from neighbouring chopper openings.

Optionally masking out large uncertainties

It is actually possible to mask out regions of large uncertainty using the LookupTableRelativeErrorThreshold parameter.

Because we may want to use different uncertainty criteria for different experimental runs as well as different components in the beamline (monitors, multiple detector banks), the masking is not hard-coded in the table but a parameter that is applied on-the-fly in the original workflow which computes time-of-flight.

We thus first update that workflow by setting the newly computed table as the TofLookupTable parameter. Next, we apply a threshold for the detector component and inspect the masked table:

wf[TofLookupTable] = table
wf[LookupTableRelativeErrorThreshold] = {"detector": 0.01}

masked = wf.compute(ErrorLimitedTofLookupTable[snx.NXdetector])

plot_lut(masked)

As we can see in the left panel, the regions in the table with large uncertainties are now white (they have been assigned a NaN value). However, the masking is not very smooth, and it would also mean that a large number of events in these white regions would be discarded in the final reduced result.

The contamination from neighbouring chopper openings is actually relatively small, and we instead decide to live with that uncertainty and set the threshold to infinity to remove any masking in the table.

wf[LookupTableRelativeErrorThreshold] = {"detector": np.inf}

Computing a time-of-flight coordinate

We will now compute our event data with a time-of-flight coordinate:

tofs = wf.compute(TofDetector[SampleRun])
tofs

scipp.DataArray (7.79 MB)

• Dimensions:
  • detector_number: 1
• Coordinates: (1)
  • Ltotal
    ()
    float64
    m
    26.0

    Values:
    array(26.)

• Data: (1)
  • (detector_number)
    float64
    counts
    binned data [len=204163]

    dim='event',
    content=DataArray(
              dims=(event: 204163),
              data=float64[counts],
              coords={'id':int64, 'eto':float64[µs], 'event_time_zero':datetime64[µs],
                      'tof':float64[µs]})

Histogramming the data for a plot should show a profile with 6 bumps that correspond to the frames:

tofs.bins.concat().hist(tof=300).plot()

Converting to wavelength

We can now convert our new time-of-flight coordinate to a neutron wavelength, using tranform_coords:

from scippneutron.conversion.graph.beamline import beamline
from scippneutron.conversion.graph.tof import elastic

# Perform coordinate transformation
graph = {**beamline(scatter=False), **elastic("tof")}
wav_wfm = tofs.transform_coords("wavelength", graph=graph)

# Define wavelength bin edges
wavs = sc.linspace("wavelength", 2, 10, 301, unit="angstrom")

histogrammed = wav_wfm.hist(wavelength=wavs).squeeze()
histogrammed.plot()

Comparing to the ground truth

As a consistency check, because we actually know the wavelengths of the neutrons we created, we can compare the true neutron wavelengths to those we computed above.

ground_truth = ess_beamline.model_result["detector"].data.flatten(to="event")
ground_truth = ground_truth[~ground_truth.masks["blocked_by_others"]]

pp.plot(
    {
        "wfm": histogrammed,
        "ground_truth": ground_truth.hist(wavelength=wavs),
    }
)