QENS data reduction

QENS data reduction#

This notebook will guide you through the data reduction for the QENS experiment that you simulated with McStas yesterday.

The following is a basic outline of what this notebook will cover:

Loading the NeXus files that contain the data
Inspect/visualize the data contents
How to convert the raw time-of-flight coordinate to something more useful (\(\Delta E\))
Write the results to file

import matplotlib.pyplot as plt
import numpy as np
import scipp as sc
import plopp as pp
import scippneutron as scn
import qens_utils as utils

Process the elastic sample data#

Visualize the raw data#

Here is a histogram of the raw data in tof:

events.hist(tof=100).plot()

../_images/dd8a3998e86b2efc5763988e7aea9bc334d03542f22d6e5258be5eb082ad2f2d.svg

It is more enlightening to look at a 2D histogram in tof and y as this shows a defect of the detector for a certain y-range:

events.hist(tof=100, y=100).plot()

../_images/89b8586c2b9418dbc65516ac817d1e55f2c786c818e07d8c710b07ff6dd26adb.svg

Transform to energy transfer and momentum transfer#

The next step is to compute the measured energy transfer and momentum transfer in the sample from time-of-flight and the other coordinates of our data. We use scipp.transform_coords for this purpose.

ScippNeutron provides some pre-defined functions for coordinate transformations, see the documentation for a list. In particular, it provides scippneutron.conversion.tof.energy_transfer_indirect_from_tof which computes the energy transfer for an indirect-geometry spectrometer:

\[\Delta E = \frac{m_n L_1^2}{2 {(t-t_0)}^2} - \mathsf{final\_energy}\,, \qquad t_0 = \sqrt{\frac{L_2^2 m_n}{2 \mathsf{final\_energy}}}\,,\]

where \(m_n\) is the neutron mass, \(L_1\) is the primary flight path length (from source to sample), and \(L_2\) is the secondary flight path length (from sample to detector). The intermediate variable \(t_0\) is the the time of flight for the secondary path such that \((t - t_0)\) is the arrival time at the sample.

While ScippNeutron provides most of what we need, we have to define some custom components for our specific instrument. We require functions to compute

L2, the path length of the scattered neutron. The default function in ScippNeutron assumes a straight path but here, we need to take the reflection from the analyzer crystal into account: \(L_2 = |\overline{\mathsf{sample},\mathsf{analyzer}}| + |\overline{\mathsf{analyser},\mathsf{detector}}|\)
final_wavelength, the neutron wavelength that the analyzer selects for. We compute it from the known \(d\)-spacing of the analyzer \(d_a\) and the scattering angle off the analyzer \(\theta_a\): \(\mathsf{final\_wavelength} = 2 d_a \sin{\theta_a}\)
Q, the momentum transfer: \(Q = |k_i - k_f|\), where \(k_i\) and \(k_f\) are the incident and final wavevectors, respectively. Here, the neutrons scatter almost elasitcally, such that we use \(|k_i| = |k_f| = mathsf{final\_wavelength}^{-1}\). We define a function to compute \(Q\) already now but will only use it later.
final_energy, the energy of the neutrons that hit the detector (more below): \(\mathsf{final\_energy} = \displaystyle\frac{m_n}{2} v_f^2\)

The first three are implemented below. They use the positions and analyzer parameters in the input data.

def backscattered_l2(position, sample_position, analyzer_position):
    """
    Compute the length of the secondary flight path for backscattering off an analyzer.
    """
    return sc.norm(position - analyzer_position) + sc.norm(
        analyzer_position - sample_position
    )


def wavelength_from_analyzer(analyzer_dspacing, analyzer_angle):
    """
    Compute the neutron wavelength after scattering from the analyzer's d-spacing.

    Assuming Bragg scattering in the analyzer, the wavelength is
        wavelength = 2 * d * sin(theta)

    Where
        d is the analyzer's d-spacing,
        theta is the scattering angle or equivalently, the tilt of the analyzer
              w.r.t. to the sample-analyzer axis.
    """
    # 2*theta is the angle between transmitted and scattered beam.
    return (
        2
        * analyzer_dspacing
        * sc.sin(sc.scalar(np.pi / 2, unit="rad") - analyzer_angle.to(unit="rad"))
    )


def momentum_transfer(
    incident_beam, sample_position, analyzer_position, final_wavelength
):
    """
    Compute the absolute value of the momentum transfer.
    """
    scattered_beam = analyzer_position - sample_position
    kf = scattered_beam / sc.norm(scattered_beam) / final_wavelength

    # The data are close to the elastic line
    # => use incident_wavelength = final_wavelength
    ki = incident_beam / sc.norm(incident_beam) / final_wavelength

    return sc.norm(ki - kf)

We can start assembling the graph to compute the energy transfer:

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

graph = {
    # Functions for the beamline geometry.
    **beamline(scatter=True),
    # The default energy transfer function.
    "energy_transfer": energy_transfer_indirect_from_tof,
    # Replace L2 and two_theta with our own implementations.
    "L2": backscattered_l2,
    # Insert new functions for the wavelength and Q.
    "final_wavelength": wavelength_from_analyzer,
    "Q": momentum_transfer,
}
# Optional: remove unused functions in order to clean up the image below.
del graph["scattered_beam"]
del graph["two_theta"]
del graph["Ltotal"]
sc.show_graph(graph, simplified=True)

../_images/615a1425095019b71cd164daed60a7a79ee3d7dd6226b101a6d1b668d327098c.svg

Exercise 3: Compute energy and momentum transfer for all samples#

We have only looked at one of our three samples so far. Now, we repeat the same procedure for the remaining two.

Your task is to load the other samples, mask the broken detector region, and compute energy_transfer and Q. (You don’t need to repeat the calculation for the unmasked data.)

Hints:

Write a function that encapsulates the whole procedure.
Store the results for all samples in a dict. This way, you can easily plot it using sc.plot.
Use a scalar variable with a unit, e.g. sc.scalar(0.01, unit="meV"), for the histogramming to get a constant bin width.

folders = (
    "../3-mcstas/QENS_elastic",
    "../3-mcstas/QENS_known_quasi_elastic",
    "../3-mcstas/QENS_unknown_quasi_elastic",
)

Solution:

Save result to disk#

Finally, we need to save our results to disk so that the reduced data can be forwarded to the next step in the pipeline (data analysis). We will use a Scipp HDF5 file for this.

This code assumes that you stored your histograms in a dict called results. If this is not the case, you need to modify the code below.

from scippneutron.io import save_xye

for name, result in results.items():
    # Rename 'energy_transfer' to 'E' because the
    # analysis software expects this name.
    data = result.rename(energy_transfer="E")
    data.save_hdf5(f"energy_transfer_{name}.h5")

Sciline workflow#

Pause here if you have not yet been introduced to Sciline!

Below we build a re-usable Sciline workflow for reducing the QENS data. We wrap the operations that were carried out throughout this notebook into function that have the correct type-annotations, building a pipeline of connected steps.

The custom types have been created for you in the qens_utils module, and are imported from there.

Take some time to read through the code and understand the different parts. Ask questions if anything is unclear.

import sciline as sl
from qens_utils import *


def load(folder: Foldername) -> RawData:
    """Load raw data from file"""
    return RawData(utils.load_qens(folder))


def mask_bad_region(
    data: RawData, masked_range: MaskedRange, masked_bank: MaskedBank
) -> MaskedData:
    y = data.coords["y"]
    mask = (
        (y >= masked_range[0])
        & (y < masked_range[1])
        & (data.coords["bank"] == masked_bank)
    )
    return MaskedData(data.assign_masks(bad_timing=mask))


def to_Q_and_E(
    masked: MaskedData, graph: CoordTransformGraph
) -> QEData:
    return QEData(
        masked.transform_coords(["energy_transfer", "Q"], graph=graph)
    )


def to_histogram(
    data: QEData, bin_width: BinWidth
) -> QEHistogram:
    return QEHistogram(data.group("Q").hist(energy_transfer=bin_width))


wf = sl.Pipeline((load, mask_bad_region, to_Q_and_E, to_histogram))

wf.visualize()

../_images/dbe9b0e7a76d3189e4976485863fde916933caff9e24f8bf80d09d90071e9f6e.svg

Exercise 4: Set the workflow parameters#

In the visualization above, the red boxes indicate missing values; parameters required by the computation that have yet to be set on the workflow.

Using the wf[TYPE] = VALUE syntax, set the missing parameters on the workflow.

Once you have done so, use wf.compute(QEHistogram) multiple times to compute the final result for each sample, making sure to change the folder to load the data from each time.

Make one figure per sample.

Solution:

Exercise 5: Make a \((y, \Delta E)\) map#

We wish to make a 2d plot of counts as a function of \(y\) and energy transfer for one of runs (e.g. elastic). You have to:

identify one of the intermediate results that would contain the necessary information
compute that intermediate result
histogram the result (200 bins in y and 200 bins in energy_transfer will do nicely)
plot it on a figure

Hint: The solution is very simple, and should only be 2-3 lines of python at most!

Solution:

QENS data reduction

Contents

QENS data reduction#

Process the elastic sample data#

Load raw data#

Visualize the raw data#

Exercise 1: mask bad region#

Exercise 1.1: find bad region#

Exercise 1.2: create and apply mask#

Exercise 1.3: explore the masked data (optional)#

Transform to energy transfer and momentum transfer#

Exercise 2: Compute energy transfer#

Exercise 2.1: Finish the graph#

Exercise 2.2: Compute energy transfer with the masked data#

Exercise 2.3: Compute energy transfer for unmasked data#

Exercise 2.4: Compute energy and momentum transfer#

Exercise 2.5: Plot different Q’s separately (optional)#

Exercise 3: Compute energy and momentum transfer for all samples#

Save result to disk#

Sciline workflow#

Exercise 4: Set the workflow parameters#

Exercise 5: Make a \((y, \Delta E)\) map#