BIFROST reduction of simulated data#

This notebook demonstrates the basic data reduction workflow for BIFROST. It uses data that was simulated with McStas and a dedicated workflow that can process McStas data.

[1]:

import numpy as np
import scipp as sc
import sciline

from ess import bifrost
from ess.bifrost.data import (
    simulated_elastic_incoherent_with_phonon,
    lookup_table_simulation
)
from ess.spectroscopy.types import *

Construct the workflow which is a sciline.Pipeline and encodes the entire reduction procedure. We need to provide a couple of parameters so we can run the workflow:

[2]:

workflow = bifrost.BifrostSimulationWorkflow()
# Set the input file name:
workflow[Filename[SampleRun]] = simulated_elastic_incoherent_with_phonon()
# Set the lookup table for frame unwrapping:
workflow[LookupTableFilename] = lookup_table_simulation()
# Set thresholds to mask wavelength regions with large uncertainties.
# For this example, we do not mask anything:
workflow[LookupTableRelativeErrorThreshold] = {
    'detector': np.inf,
    'normalization_monitor': np.inf,
}
# We need to read many objects from the file,
# keeping it open improves performance: (optional)
workflow[PreopenNeXusFile] = PreopenNeXusFile(True)
# We drop uncertainties where they would otherwise lead to correlations:
workflow[UncertaintyBroadcastMode] = UncertaintyBroadcastMode.drop

Next, draw the workflow as a graph to inspect the steps it will take to reduce the data.

Note the groups where entries are labeled with (triplet). (These labels will be (dim_0) if you don’t have Pandas installed.) These groups indicate functions that are applied separately to each detector triplet.

[3]:

workflow.visualize(EnergyQDetector[SampleRun], compact=True, graph_attr={"rankdir": "LR"})

[3]:

../../_images/user-guide_bifrost_bifrost-reduction_5_0.svg

We are ready to compute the reduced data. We use the naive scheduler of sciline because it tends to perform better for BIFROST than the Dask scheduler. But this is optional.

[4]:

scheduler = sciline.scheduler.NaiveScheduler()
data = workflow.compute(EnergyQDetector[SampleRun], scheduler=scheduler)

The result contains coordinates for the sample table and detector rotation angles a3 and a4, respectively. It also contains event coordinates for energy_transfer and two momentum transfers, one in the lab frame and one in the sample table frame.

[5]:

data

[5]:

We can plot the counts as a function of energy transfer and \(a_3\) by removing the unused dimensions. As expected, it is independent of \(a_3\).

[6]:

(
    data
    .bins.concat()
    .hist(energy_transfer=500)
).plot()

[6]:

../../_images/user-guide_bifrost_bifrost-reduction_11_0.svg

We can also plot the counts as a function of the momentum transfer in the sample table frame \(Q\). For this, we first need to create a 2D slice in \(Q\) and \(\Delta E\). For simplicity, we use the x and z axes or \(Q\) (see https://scipp.github.io/scippneutron/user-guide/coordinate-transformations.html for definitions). But we could use any other normalized, orthogonal vectors in the dot products.

We can make a second plot using one of the projections and the energy transfer.

The plots are a bit coarse because we only used two triplets.

[7]:

d = data.bins.concat().copy()
x = sc.vector([1, 0, 0])
z = sc.vector([0, 0, 1])
d.bins.coords['Qx'] = sc.dot(x, d.bins.coords['sample_table_momentum_transfer'])
d.bins.coords['Qz'] = sc.dot(z, d.bins.coords['sample_table_momentum_transfer'])
d.hist(Qz=100, Qx=100).plot(norm='log') / d.hist(energy_transfer=100, Qz=100).plot(norm='log')

[7]:

../../_images/user-guide_bifrost_bifrost-reduction_13_0.svg