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> ()    Dataset        <no unit>  [<scipp.Dataset>
 Dimensions: Sizes[frame:6, ]
 Data:
   cutout_angles_center      float64            [rad]  (frame)  [1.60454, 2.48242, ..., 4.80143, 5.48021]
   cutout_angles_width       float64            [rad]  (frame)  [0.100208, 0.156644, ..., 0.305723, 0.349358]
   frequency                 float64             [Hz]  ()  [56]
   kind                       string        <no unit>  ()  ["wfm"]
   phase                     float64            [deg]  ()  [0]
   position                  vector3              [m]  ()  [(0, 0, 6.775)]

 ],
 'chopper_wfm_2': <scipp.Variable> ()    Dataset        <no unit>  [<scipp.Dataset>
 Dimensions: Sizes[frame:6, ]
 Data:
   cutout_angles_center      float64            [rad]  (frame)  [1.70474, 2.63907, ..., 5.10716, 5.82957]
   cutout_angles_width       float64            [rad]  (frame)  [0.100208, 0.156644, ..., 0.305723, 0.349358]
   frequency                 float64             [Hz]  ()  [56]
   kind                       string        <no unit>  ()  ["wfm"]
   phase                     float64            [deg]  ()  [0]
   position                  vector3              [m]  ()  [(0, 0, 7.225)]

 ],
 '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 useable 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 0x7f5b7856da60>

../../_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"].data) - intercepts) / slopes
    times_at_wfm2 = (sc.norm(far_wfm_chopper["position"].data) - 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]).data &
                 (item["wavelengths"] <= frames["wavelength_max"]["frame", i]).data)
    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]:

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]:

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_4099/3750900606.py:3: RuntimeWarning: invalid value encountered in true_divide
  values=events["sample"]["wavelength_counts"] /

[21]:

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.

Working in event mode#

It is also possible to work with WFM data in event mode. The stitch utility will accept both histogrammed and binned (event) data.

We first create a new dataset, with the same events as in the first example, but this time we bin the data with sc.bin instead of using sc.histogram, so we can retain the raw events.

[29]:

for item in events.values():
    item["valid_times"] = item["arrival_times"].values[item["valid_indices"]]

tmin = min([item["valid_times"].min() for item in events.values()])
tmax = max([item["valid_times"].max() for item in events.values()])

dt = 0.1 * (tmax - tmin)
time_coord = sc.linspace(dim='time',
                         start=tmin - dt,
                         stop=tmax + dt,
                         num=257,
                         unit=events["sample"]["arrival_times"].unit)

ds_event = sc.Dataset(coords=coords)

# Bin the data
for key, item in events.items():
    da = sc.DataArray(
        data=sc.ones(dims=['event'], shape=[len(item["valid_times"])], unit=sc.units.counts,
                     with_variances=True),
        coords={
            'time': sc.array(dims=['event'], values=item["valid_times"], unit=sc.units.us)})
    ds_event[key] = da.bin(time=time_coord)

ds_event

The underlying events can be inspected by using the .bins.constituents['data'] property of our objects:

[30]:

ds_event["sample"].bins.constituents['data']

[30]:

scipp.DataArray (4.49 MB)

Dimensions:
- event: 196024

Coordinates: (1)

time

(event)

float64

µs

1.812e+04, 1.822e+04, ..., 1.333e+05, 1.334e+05

Values:
array([ 18119.40258962,  18224.68653735,  18104.79393712, ...,
       133196.6964019 , 133288.52178163, 133352.36035354])

Data: (1)

(event)

float64

counts

1.0, 1.0, ..., 1.0, 1.0

σ = 1.0, 1.0, ..., 1.0, 1.0

Values:
array([1., 1., 1., ..., 1., 1., 1.])

Variances (σ²):
array([1., 1., 1., ..., 1., 1., 1.])

We can visualize this to make sure it looks the same as the histogrammed case above:

[31]:

sc.plot(ds_event.hist())

[31]:

As explained above, the stitch routine accepts both histogrammed and binned (event) data. So stitching the binned data works in the exact same way as above, namely

[32]:

stitched_event = wfm.stitch(frames=frames,
                            data=ds_event,
                            dim='time')
stitched_event

The stitch function will return a data structure with a single bin in the 'tof' dimension. Visualizing this data is therefore slightly more tricky, because the data needs to be histogrammed using a finer binning before a useful plot can be made.

[33]:

sc.plot({key: item.hist(tof=256) for key, item in stitched_event.items()})

[33]:

At this point, it may be useful to compare the results of the two different stitching operations.

[34]:

rebinned = stitched_event["sample"].bin(tof=stitched["sample"].coords['tof'])
sc.plot({"events": rebinned.hist(), "histogram": stitched["sample"]}, errorbars=False)

[34]:

We note that histogramming the data early introduces some smoothing in the data.

We can of course continue in event mode and perform the unit conversion and normalization to the Vanadium.

[35]:

converted_event = stitched_event.transform_coords("wavelength", graph=graph)

[36]:

# Normalizing binned data is done using the sc.lookup helper
hist = converted_event["vanadium"].hist(wavelength=converted["sample"].coords['wavelength'])
normalized_event = converted_event["sample"].bins / sc.lookup(func=hist, dim='wavelength')

Finally, we compare the end result with the original wavelengths, and see that the agreement is once again good.

[37]:

to_plot = normalized_event.bin(wavelength=converted["sample"].coords['wavelength']).hist()
sc.plot({"stitched_event": to_plot, "original": original})

[37]:

We can also compare directly to the histogrammed version, to see that both methods remain in agreement to a high degree.

[38]:

sc.plot({"stitched": normalized['wavelength', w_min:w_max],
         "original": original['wavelength', w_min:w_max],
         "stitched_event": to_plot['wavelength', w_min:w_max]})

[38]:

Reducing WFM data

Contents

Reducing WFM data#

Create beamline components#

Generate some fake data#

Visualize the beamline’s chopper cascade#

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

Create a realistic Dataset#

Stitch the frames#

Convert to wavelength#

Normalization#

Comparing to the raw wavelengths#

Without WFM choppers#

Working in event mode#