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 sciline as sl
from scippneutron.chopper import DiskChopper
from ess.reduce import time_of_flight

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

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,
    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(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(204430.0)

Computing time-of-flight

The chopper information is next used to construct a lookup table that provides an estimate of the real time-of-flight as a function of neutron time-of-arrival.

We use the tof module 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.

Setting up the workflow

workflow = sl.Pipeline(
    time_of_flight.providers(), params=time_of_flight.default_parameters()
)
workflow[time_of_flight.RawData] = raw_data
workflow[time_of_flight.LtotalRange] = Ltotal, Ltotal
workflow[time_of_flight.LookupTableRelativeErrorThreshold] = 0.1

workflow.visualize(time_of_flight.TofData)

We can see from the workflow diagram that we are still missing the simulated neutrons that are used to build the lookup table.

Those are obtained by running a quick tof simulation of the beamline:

workflow[time_of_flight.SimulationResults] = time_of_flight.simulate_beamline(
    choppers=disk_choppers,
    neutrons=3_000_000
)

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 = workflow.compute(time_of_flight.SimulationResults)
# Compute time-of-arrival at the detector
tarrival = sim.time_of_arrival + ((Ltotal - sim.distance) / sim.speed).to(unit="us")
# Compute time-of-flight at the detector
tflight = (Ltotal / sim.speed).to(unit="us")

events = sc.DataArray(
    data=sim.weight,
    coords={"wavelength": sim.wavelength, "toa": tarrival, "tof": tflight},
)
fig1 = events.hist(wavelength=300, toa=300).plot(norm="log")
fig2 = events.hist(tof=300, toa=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 = workflow.compute(time_of_flight.TimeOfFlightLookupTable).squeeze()

# Overlay mean on the figure above
table["distance", 1].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.

Computing a time-of-flight coordinate

We will now use our workflow to obtain our event data with a time-of-flight coordinate:

tofs = workflow.compute(time_of_flight.TofData)
tofs

scipp.DataArray (7.80 MB)

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

    Values:
    array([26.])

• Data: (1)
  • (detector_number)
    DataArrayView
    binned data [len=204430]

    dim='event',
    content=DataArray(
              dims=(event: 204430),
              data=float64[counts],
              coords={'id':int64, 'event_time_zero':datetime64[µs], 'event_time_offset':float64[µ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),
    }
)