BEER McStas powder data reduction#

Audience: BEER users reducing McStas pulse-shaping powder data
Prerequisites: Basic knowledge of Scipp and Sciline workflows

This guide reduces a McStas silicon-in-vanadium-can sample run using matching vanadium and empty-can runs. The result is EmptyCanSubtractedIofDspacing, the normalized intensity as a function of \(d\)-spacing after subtracting the empty can.

The example files are downloaded through the ESSdiffraction data registry. To run the notebook you will have to install pooch to access them: pip install pooch.

[1]:

%matplotlib ipympl
import scipp as sc

from ess import beer, powder
import ess.beer.data  # noqa: F401
from ess.beer.types import DetectorBank
from ess.diffraction.peaks import dspacing_peaks_from_cif
from ess.powder.types import *

Configure the workflow#

Start with the BEER powder McStas workflow and set the run files. We use monitor-histogram normalization here because the McStas files include a wavelength monitor.

[2]:

workflow = beer.BeerPowderMcStasWorkflow(
    run_norm=powder.RunNormalization.monitor_histogram
)

workflow[Filename[SampleRun]] = beer.data.mcstas_powder_silicon_in_vanadium_can()
workflow[Filename[VanadiumRun]] = beer.data.mcstas_powder_vanadium()
workflow[Filename[EmptyCanRun]] = beer.data.mcstas_powder_empty_can()

Downloading file 'beer-powder-si-in-vcan-mode4.h5' from 'https://public.esss.dk/groups/scipp/ess/beer/1/beer-powder-si-in-vcan-mode4.h5' to '/home/runner/.cache/ess/beer'.
Downloading file 'beer-powder-vanadium-mode4.h5' from 'https://public.esss.dk/groups/scipp/ess/beer/1/beer-powder-vanadium-mode4.h5' to '/home/runner/.cache/ess/beer'.
Downloading file 'beer-powder-empty-can-mode4.h5' from 'https://public.esss.dk/groups/scipp/ess/beer/1/beer-powder-empty-can-mode4.h5' to '/home/runner/.cache/ess/beer'.

Next set the detector bank and the \(d\)-spacing range. The notebook only reduces the north bank and only histograms the range from 0.6 Å to 2.0 Å.

[3]:

workflow[DetectorBank] = DetectorBank.north
workflow[DspacingBins] = sc.linspace("dspacing", 0.6, 2.0, 701, unit="angstrom")

To keep the events in the output, to allow post reduction filtering or other processing, set KeepEvents to True.

[4]:

workflow[KeepEvents[SampleRun]] = KeepEvents(True)
workflow[KeepEvents[VanadiumRun]] = KeepEvents(True)
workflow[KeepEvents[EmptyCanRun]] = KeepEvents(True)

workflow[MaskedDetectorIDs] = MaskedDetectorIDs({})
workflow[UncertaintyBroadcastMode] = UncertaintyBroadcastMode.drop

A finer region of interest can be applied with the standard powder masks. Mask callables return True for values to remove. This example keeps a broad two-theta range for the selected detector data and leaves the other masks disabled.

[5]:

two_theta_min = sc.scalar(76.0, unit="deg").to(unit="rad")
two_theta_max = sc.scalar(104.0, unit="deg").to(unit="rad")

workflow[TwoThetaMask] = lambda two_theta: (two_theta < two_theta_min) | (
    two_theta > two_theta_max
)
workflow[TofMask] = None
workflow[WavelengthMask] = None

Some events can have ambiguous wavelength assignments. We first keep all detector events so the baseline spectrum shows the full result.

[6]:

workflow[LookupTableRelativeErrorThreshold] = {
    "detector": float("inf"),
    "monitor_bunker": float("inf"),
    "monitor_cave": float("inf"),
}

Reduce frame-overlap artifacts#

The detector wavelength-assignment relative-error threshold removes events in regions where the wavelength assignment is ambiguous. The thresholds are set separately for the detector and the monitors.

[9]:

thresholded_workflow = workflow.copy()
thresholded_workflow[LookupTableRelativeErrorThreshold] = {
    "detector": 0.005,
    "monitor_bunker": float("inf"),
    "monitor_cave": float("inf"),
}

thresholded_histogram = thresholded_workflow.compute(
    EmptyCanSubtractedIofDspacing
).hist()

The threshold removes events in specific regions of scattering angle and wavelength. We can see this directly by comparing CorrectedDetector before and after applying the threshold.

[10]:

unfiltered_detector = workflow.compute(CorrectedDetector[SampleRun])
filtered_detector = thresholded_workflow.compute(CorrectedDetector[SampleRun])

two_theta_bins = sc.linspace("two_theta", 76.0, 104.0, 201, unit="deg").to(
    unit="rad"
)
wavelength_bins = sc.linspace("wavelength", 0.8, 3.2, 241, unit="angstrom")

unfiltered_map = unfiltered_detector.hist(
    two_theta=two_theta_bins,
    wavelength=wavelength_bins,
    dim=unfiltered_detector.dims,
)
filtered_map = filtered_detector.hist(
    two_theta=two_theta_bins,
    wavelength=wavelength_bins,
    dim=filtered_detector.dims,
)

[11]:

filtered_map.plot(title="Kept by threshold", norm="log")

[11]:

[12]:

(unfiltered_map - filtered_map).plot(title="Removed by threshold", norm="log")

[12]:

[13]:

fig = sc.plot(
    {
        "baseline": histogram,
        "detector relative error < 0.02": thresholded_histogram,
    },
    xlabel=r"$d$-spacing [$\AA$]",
    ylabel=f"Intensity [{histogram.unit}]",
    xmin=sc.scalar(0.6, unit="angstrom"),
    xmax=sc.scalar(2.0, unit="angstrom"),
)

ymin, ymax = fig.ax.get_ylim()
for peak in silicon_peaks.values:
    fig.ax.axvline(peak, color="black", alpha=0.25, linewidth=1)

fig.ax.set_ylim(ymin, ymax)
fig

[13]:

BEER McStas powder data reduction#

Configure the workflow#

Compute the empty-can-subtracted spectrum#

Reduce frame-overlap artifacts#