Reducing WFM data#

This notebook aims to illustrate how to work with the wavelength frame multiplication submodule wfm.

We will create a beamline that resembles the ODIN instrument beamline, generate some fake neutron data, and then show how to convert the neutron arrival times at the detector to neutron time-of-flight, from which a wavelength can then be computed (or process also commonly known as ‘stitching’).

[1]:

import numpy as np
import matplotlib.pyplot as plt
import scipp as sc
from scipp import constants
import scippneutron as scn
import ess.wfm as wfm
import ess.choppers as ch

np.random.seed(1)  # Fixed for reproducibility

Create beamline components#

We first create all the components necessary to a beamline to run in WFM mode (see Introduction to WFM for the meanings of the different symbols). The beamline will contain

a neutron source, located at the origin (\(x = y = z = 0\))
a pulse with a defined length (\(2860 ~\mu s\)) and \(t_0\) (\(130 ~\mu s\))
a single pixel detector, located at \(z = 60\) m
two WFM choppers, located at \(z = 6.775\) m and \(z = 7.225\) m, each with 6 frame windows/openings

The wfm module provides a helper function to quickly create such a beamline. It returns a dict of coordinates, that can then be subsequently added to a data container.

[2]:

coords = wfm.make_fake_beamline(nframes=6)
coords

[2]:

{'chopper_wfm_1': <scipp.Variable> ()   PyObject        <no unit>  DataGroup(sizes={'frame': 6}, keys=['frequency', 'position', 'phase', 'cutout_angles_center', 'cutout_angles_width', 'kind']),
 'chopper_wfm_2': <scipp.Variable> ()   PyObject        <no unit>  DataGroup(sizes={'frame': 6}, keys=['frequency', 'position', 'phase', 'cutout_angles_center', 'cutout_angles_width', 'kind']),
 'position': <scipp.Variable> ()    vector3              [m]  (0, 0, 60),
 'source_pulse_length': <scipp.Variable> ()    float64            [µs]  2860,
 'source_pulse_t_0': <scipp.Variable> ()    float64            [µs]  130,
 'source_position': <scipp.Variable> ()    vector3              [m]  (0, 0, 0)}

Generate some fake data#

Next, we will generate some fake imaging data (no scattering will be considered), that is supposed to mimic a spectrum with a Bragg edge located at \(4\unicode{x212B}\). We start with describing a function which will act as our underlying distribution

[3]:

x = np.linspace(1, 10.0, 100000)
a = 20.0
b = 4.0
y1 = 0.7 / (np.exp(-a * (x - b)) + 1.0)
y2 = 1.4 - 0.2 * x
y = y1 + y2
fig1, ax1 = plt.subplots()
ax1.plot(x, y)
ax1.set_xlabel("Wavelength [angstroms]")

[3]:

Text(0.5, 0, 'Wavelength [angstroms]')

../../_images/techniques_wfm_reducing-wfm-data_5_1.png

We then proceed to generate two sets of 1,000,000 events: - one for the sample using the distribution defined above - and one for the vanadium which will be just a flat random distribution

For the events in both sample and vanadium, we define a wavelength for the neutrons as well as a birth time, which will be a random time between the pulse \(t_0\) and the end of the usable pulse \(t_0\) + pulse_length.

[4]:

nevents = 1_000_000
events = {
    "sample": {
        "wavelengths": sc.array(
            dims=["event"],
            values=np.random.choice(x, size=nevents, p=y / np.sum(y)),
            unit="angstrom",
        ),
        "birth_times": sc.array(
            dims=["event"],
            values=np.random.random(nevents) * coords["source_pulse_length"].value,
            unit="us",
        )
        + coords["source_pulse_t_0"],
    },
    "vanadium": {
        "wavelengths": sc.array(
            dims=["event"],
            values=np.random.random(nevents) * 9.0 + 1.0,
            unit="angstrom",
        ),
        "birth_times": sc.array(
            dims=["event"],
            values=np.random.random(nevents) * coords["source_pulse_length"].value,
            unit="us",
        )
        + coords["source_pulse_t_0"],
    },
}

We can then take a quick look at our fake data by histogramming the events

[5]:

# Histogram and plot the event data
bins = np.linspace(1.0, 10.0, 129)
fig2, ax2 = plt.subplots()
for key in events:
    h = ax2.hist(events[key]["wavelengths"].values, bins=128, alpha=0.5, label=key)
ax2.set_xlabel("Wavelength [angstroms]")
ax2.set_ylabel("Counts")
ax2.legend()

[5]:

<matplotlib.legend.Legend at 0x7f54de2ebb20>

../../_images/techniques_wfm_reducing-wfm-data_9_1.png

We can also verify that the birth times fall within the expected range:

[6]:

for key, item in events.items():
    print(key)
    print(sc.min(item["birth_times"]))
    print(sc.max(item["birth_times"]))

sample
<scipp.Variable> ()    float64            [µs]  130.008
<scipp.Variable> ()    float64            [µs]  2990
vanadium
<scipp.Variable> ()    float64            [µs]  130.004
<scipp.Variable> ()    float64            [µs]  2990

We can then compute the arrival times of the events at the detector pixel

[7]:

# The ratio of neutron mass to the Planck constant
alpha = sc.to_unit(constants.m_n / constants.h, 'us/m/angstrom')
# The distance between the source and the detector
dz = sc.norm(coords['position'] - coords['source_position'])
for key, item in events.items():
    item["arrival_times"] = alpha * dz * item["wavelengths"] + item["birth_times"]
events["sample"]["arrival_times"]

[7]:

scipp.Variable (7.63 MB)

(event: 1000000)

float64

µs

6.326e+04, 9.431e+04, ..., 2.788e+04, 5.809e+04

Values:
array([63257.61773626, 94307.95311489, 16899.66119625, ...,
       70065.41396431, 27877.08441761, 58091.769017  ])

Discard neutrons that do not make it through the chopper windows#

Once we have the parameters of the 6 wavelength frames, we need to run through all our generated neutrons and filter out all the neutrons with invalid flight paths, i.e. the ones that do not make it through both chopper openings in a given frame.

[11]:

for item in events.values():
    item["valid_indices"] = []
near_wfm_chopper = ds.coords["chopper_wfm_1"].value
far_wfm_chopper = ds.coords["chopper_wfm_2"].value
near_time_open = ch.time_open(near_wfm_chopper)
near_time_close = ch.time_closed(near_wfm_chopper)
far_time_open = ch.time_open(far_wfm_chopper)
far_time_close = ch.time_closed(far_wfm_chopper)

for item in events.values():
    # Compute event arrival times at wfm choppers 1 and 2
    slopes = 1.0 / (alpha * item["wavelengths"])
    intercepts = -slopes * item["birth_times"]
    times_at_wfm1 = (sc.norm(near_wfm_chopper["position"]) - intercepts) / slopes
    times_at_wfm2 = (sc.norm(far_wfm_chopper["position"]) - intercepts) / slopes
    # Create a mask to see if neutrons go through one of the openings
    mask = sc.zeros(dims=times_at_wfm1.dims, shape=times_at_wfm1.shape, dtype=bool)
    for i in range(len(frames["time_min"])):
        mask |= (
            (times_at_wfm1 >= near_time_open["frame", i])
            & (times_at_wfm1 <= near_time_close["frame", i])
            & (item["wavelengths"] >= frames["wavelength_min"]["frame", i])
            & (item["wavelengths"] <= frames["wavelength_max"]["frame", i])
        )
    item["valid_indices"] = np.ravel(np.where(mask.values))

Convert to wavelength#

Now that the data coordinate is time-of-flight (tof), we can use scippneutron to perform the unit conversion from tof to wavelength.

[18]:

from scippneutron.conversion.graph import beamline, tof

graph = {**beamline.beamline(scatter=False), **tof.elastic("tof")}
converted = stitched.transform_coords("wavelength", graph=graph)
converted.plot()

[18]:

../../_images/techniques_wfm_reducing-wfm-data_32_0.svg

Normalization#

Normalization is performed simply by dividing the counts of the sample run by the counts of the vanadium run.

[19]:

normalized = converted['sample'] / converted['vanadium']
normalized.plot()

[19]:

../../_images/techniques_wfm_reducing-wfm-data_34_0.svg

Comparing to the raw wavelengths#

The final step is a sanity check to verify that the wavelength-dependent data obtained from the stitching process agrees (to within the beamline resolution) with the original wavelength distribution that was generated at the start of the workflow.

For this, we simply histogram the raw neutron events using the same bins as the normalized data, filtering out the neutrons with invalid flight paths.

[20]:

for item in events.values():
    item["wavelength_counts"], _ = np.histogram(
        item["wavelengths"].values[item["valid_indices"]],
        bins=normalized.coords['wavelength'].values,
    )

We then normalize the sample by the vanadium run, and plot the resulting spectrum alongside the one obtained from the stitching.

[21]:

original = sc.DataArray(
    data=sc.array(
        dims=['wavelength'],
        values=events["sample"]["wavelength_counts"]
        / events["vanadium"]["wavelength_counts"],
    ),
    coords={"wavelength": normalized.coords['wavelength']},
)

sc.plot({"stitched": normalized, "original": original})

/tmp/ipykernel_4380/3402972659.py:4: RuntimeWarning: invalid value encountered in divide
  values=events["sample"]["wavelength_counts"]

[21]:

../../_images/techniques_wfm_reducing-wfm-data_38_1.svg

We can see that the counts in the stitched data agree very well with the original data. There is some smoothing of the data seen in the stitched result, and this is expected because of the resolution limitations of the beamline due to its long source pulse. This smoothing (or smearing) would, however, be much stronger if WFM choppers were not used.