[1]:

%matplotlib ipympl

BEER (modulation) McStas data reduction#

[2]:

import scipp as sc
import scippneutron as scn

from ess.beer import BeerModMcStasWorkflow, BeerModMcStasWorkflowKnownPeaks
from ess.beer.data import mcstas_silicon_medium_resolution, mcstas_duplex, duplex_peaks_array, silicon_peaks_array
from ess.reduce.nexus.types import Filename, SampleRun
from ess.reduce.time_of_flight.types import DetectorTofData
from ess.beer.types import *

# Default bin edges for our d_hkl histograms
dspacing = sc.linspace('dspacing', 0.8, 2.2, 2000, unit='angstrom')

def ground_truth_peak_positions(p, arr):
    'Helper to display the true peak positions for comparison'
    p.ax.vlines(arr.values, 0, 20000, linestyle='--', color='black', lw=0.5, label='true $d_{hkl}$')
    return p

Beam modulation mode#

In the “modulation mode” the BEER instrument maximizes signal intensity while retaining wavelength resolution by utilizing choppers to create multiple pulses, each with well defined wavelenth, that reach the sample simultaneously. The pulses can be distinguished in the detector data under the assumption that the peaks in \(d_{hkl}\) are sufficiently well separated from each other.

Data reduction in modulation mode#

To read more about how this works, read the Beer modulation mode paper [{cite}][ROUIJAA20187].

Associating each event to the \(d_{hkl}\) of the scattering event can be done using two different methods. Which method can be used depends on if we have good prior estimates of the peak positions or not.

1. With prior peak position estimates#

Given a list of peak positions, and the geometry of the instrument/chopper system, we can figure out what regions of the plane (\(2\theta\), \(t\)) are going to be covered by events that were scattered from each particular peak (here \(2\theta\) is the scattering angle and \(t\) is the “time-of-arrival” in the detector).

From this we can associate each event with the \(d_{hkl}\) peak it scattered from, and using the relationship

\[\frac{\tau}{L} = \frac{m_n}{h} d_{hkl} \sin(\theta)\]

\(\tau\) (“time-of-flight”) can be estimated. From the time of flight a new \(d_{hkl}\) estimate can be computed for each event.

2. Without prior peak position estimates#

Without prior peak position estimates we first need to cluster the events so that each cluster contains events that were scattered from the same \(d_{hkl}\) peak and were produced in the same modulation pulse.

Within such a cluster

\[\frac{\tau}{L} = \frac{m_n}{h} d_{hkl} \sin(\theta)\]
\[\iff \frac{t}{L} = \frac{t_0}{L} + \frac{m_n}{h} d_{hkl} \sin(\theta)\]

the time of arrival \(t\) is an affine function of \(\sin(\theta)\), and \(t_0 := t - \tau\) (the origin time of the modulation pulse) can be found by linear regression within the cluster. Once the origin time of each modulation pulse is known time-of-flight can be estimated and \(d_{hkl}\) computed for each event as usual.

Difference between modulation modes#

There are different modulation modes depending on the desired trade off between resolution and intensity. Below is a table describing the available modes.

Diffraction, modulation

Description

Application

M0

Maximum intensity, modulation

diffraction, max. slit

M1

High flux, modulation

strain scanning

M2

Medium resolution, modulation

strain scanning

M3

High resolution, modulation

strain scanning

In the rest of the notebook you will find examples of simulated data for a few different samples measured in the different modulation modes.

The examples begin by displaying what the intensity distribution in the detector looked like in the experiment. Then they show the resulting \(d_{hkl}\) histograms, using the two different methods described above. Lastly is a diagram that shows what fraction of the data was excluded because it could not be associated with a unique modulation pulse time. The excluded data figure is only shown for one of the reduction methods, to keep the notebook from getting too long, but both methods might exclude data in situations such as when nearby peaks are overlapping.

Silicon sample (M2)#

Intensity distribution in detector#

[3]:

wf = BeerModMcStasWorkflowKnownPeaks()
wf[Filename[SampleRun]] = mcstas_silicon_medium_resolution()
wf.compute(DetectorData[SampleRun])['bank1'].hist(two_theta=1000, event_time_offset=1000).plot(norm='log')

Downloading file 'silicon-mode09.h5' from 'https://public.esss.dk/groups/scipp/ess/beer/1/silicon-mode09.h5' to '/home/runner/.cache/ess/beer/1'.

[3]:

Known peaks workflow#

[4]:

wf[DHKLList] = silicon_peaks_array()
da = wf.compute(DetectorTofData[SampleRun])
da = da.transform_coords(('dspacing',), graph=scn.conversion.graph.tof.elastic('tof'),)
ground_truth_peak_positions(da.hist(dspacing=dspacing, dim=da.dims).plot(), silicon_peaks_array())

Downloading file 'silicon-dhkl.tab' from 'https://public.esss.dk/groups/scipp/ess/beer/1/silicon-dhkl.tab' to '/home/runner/.cache/ess/beer/1'.

[4]:

Automatic peak finder workflow#

[5]:

wf = BeerModMcStasWorkflow()
wf[Filename[SampleRun]] = mcstas_silicon_medium_resolution()
da = wf.compute(DetectorTofData[SampleRun])
da = da.transform_coords(('dspacing',), graph=scn.conversion.graph.tof.elastic('tof'),)
ground_truth_peak_positions(da.hist(dspacing=dspacing, dim=da.dims).plot(), silicon_peaks_array())

[5]:

Data that was excluded because tof could not be computed#

The most common reason this happens is overlap between nearby peaks.

[6]:

da = da['bank2'].copy()
da.masks.clear()
da.bins.masks['too_far_from_center'] = ~da.bins.masks.pop('too_far_from_center')

da.bins.concat().hist(two_theta=1000, event_time_offset=1000).plot(norm='log')

[6]:

Duplex sample (M1)#

Intensity distribution in detector#

[7]:

wf = BeerModMcStasWorkflowKnownPeaks()
wf[Filename[SampleRun]] = mcstas_duplex(8)
wf.compute(DetectorData[SampleRun])['bank1'].hist(two_theta=1000, event_time_offset=1000).plot(norm='log')

Downloading file 'duplex-mode08.h5' from 'https://public.esss.dk/groups/scipp/ess/beer/1/duplex-mode08.h5' to '/home/runner/.cache/ess/beer/1'.

[7]:

Known peaks workflow#

[8]:

wf[DHKLList] = duplex_peaks_array()
da = wf.compute(DetectorTofData[SampleRun])
da = da.transform_coords(('dspacing',), graph=scn.conversion.graph.tof.elastic('tof'),)
ground_truth_peak_positions(da.hist(dspacing=dspacing, dim=da.dims).plot(), duplex_peaks_array())

Downloading file 'duplex-dhkl.tab' from 'https://public.esss.dk/groups/scipp/ess/beer/1/duplex-dhkl.tab' to '/home/runner/.cache/ess/beer/1'.

[8]:

Automatic peak finder workflow#

[9]:

wf = BeerModMcStasWorkflow()
wf[Filename[SampleRun]] = mcstas_duplex(8)
da = wf.compute(DetectorTofData[SampleRun])
da = da.transform_coords(('dspacing',), graph=scn.conversion.graph.tof.elastic('tof'),)
ground_truth_peak_positions(da.hist(dspacing=dspacing, dim=da.dims).plot(), duplex_peaks_array())

[9]:

Data that was excluded because tof could not be computed#

The most common reason this happens is overlap between nearby peaks.

[10]:

da = da['bank2'].copy()
da.masks.clear()
da.bins.masks['too_far_from_center'] = ~da.bins.masks.pop('too_far_from_center')

da.bins.concat().hist(two_theta=1000, event_time_offset=1000).plot(norm='log')

[10]:

Duplex sample (M2)#

Intensity distribution in detector#

[11]:

wf = BeerModMcStasWorkflowKnownPeaks()
wf[Filename[SampleRun]] = mcstas_duplex(9)
wf.compute(DetectorData[SampleRun])['bank1'].hist(two_theta=1000, event_time_offset=1000).plot(norm='log')

Downloading file 'duplex-mode09.h5' from 'https://public.esss.dk/groups/scipp/ess/beer/1/duplex-mode09.h5' to '/home/runner/.cache/ess/beer/1'.

[11]:

Known peaks workflow#

[12]:

wf[DHKLList] = duplex_peaks_array()
da = wf.compute(DetectorTofData[SampleRun])
da = da.transform_coords(('dspacing',), graph=scn.conversion.graph.tof.elastic('tof'),)
ground_truth_peak_positions(da.hist(dspacing=dspacing, dim=da.dims).plot(), duplex_peaks_array())

[12]:

Automatic peak finder workflow#

[13]:

wf = BeerModMcStasWorkflow()
wf[Filename[SampleRun]] = mcstas_duplex(9)
da = wf.compute(DetectorTofData[SampleRun])
da = da.transform_coords(('dspacing',), graph=scn.conversion.graph.tof.elastic('tof'),)
ground_truth_peak_positions(da.hist(dspacing=dspacing, dim=da.dims).plot(), duplex_peaks_array())

[13]:

Data that was excluded because tof could not be computed#

The most common reason this happens is overlap between nearby peaks.

[14]:

da = da['bank2'].copy()
da.masks.clear()
da.bins.masks['too_far_from_center'] = ~da.bins.masks.pop('too_far_from_center')

da.bins.concat().hist(two_theta=1000, event_time_offset=1000).plot(norm='log')

[14]:

Duplex sample (M3)#

Intensity distribution in detector#

[15]:

wf = BeerModMcStasWorkflowKnownPeaks()
wf[Filename[SampleRun]] = mcstas_duplex(10)
wf.compute(DetectorData[SampleRun])['bank1'].hist(two_theta=1000, event_time_offset=1000).plot(norm='log')

Downloading file 'duplex-mode10.h5' from 'https://public.esss.dk/groups/scipp/ess/beer/1/duplex-mode10.h5' to '/home/runner/.cache/ess/beer/1'.

[15]:

Known peaks workflow#

[16]:

wf[DHKLList] = duplex_peaks_array()
da = wf.compute(DetectorTofData[SampleRun])
da = da.transform_coords(('dspacing',), graph=scn.conversion.graph.tof.elastic('tof'),)
ground_truth_peak_positions(da.hist(dspacing=dspacing, dim=da.dims).plot(), duplex_peaks_array())

[16]:

Automatic peak finder workflow#

[17]:

wf = BeerModMcStasWorkflow()
wf[Filename[SampleRun]] = mcstas_duplex(10)
da = wf.compute(DetectorTofData[SampleRun])
da = da.transform_coords(('dspacing',), graph=scn.conversion.graph.tof.elastic('tof'),)
ground_truth_peak_positions(da.hist(dspacing=dspacing, dim=da.dims).plot(), duplex_peaks_array())

[17]:

Data that was excluded because tof could not be computed#

The most common reason this happens is overlap between nearby peaks.

[18]:

da = da['bank2'].copy()
da.masks.clear()
da.bins.masks['too_far_from_center'] = ~da.bins.masks.pop('too_far_from_center')

da.bins.concat().hist(two_theta=1000, event_time_offset=1000).plot(norm='log')

[18]:

Duplex sample (M4)#

Intensity distribution in detector#

[19]:

wf = BeerModMcStasWorkflowKnownPeaks()
wf[Filename[SampleRun]] = mcstas_duplex(16)
wf.compute(DetectorData[SampleRun])['bank1'].hist(two_theta=1000, event_time_offset=1000).plot(norm='log')

Downloading file 'duplex-mode16.h5' from 'https://public.esss.dk/groups/scipp/ess/beer/1/duplex-mode16.h5' to '/home/runner/.cache/ess/beer/1'.

[19]:

Known peaks workflow#

[20]:

wf[DHKLList] = duplex_peaks_array()
da = wf.compute(DetectorTofData[SampleRun])
da = da.transform_coords(('dspacing',), graph=scn.conversion.graph.tof.elastic('tof'),)
ground_truth_peak_positions(da.hist(dspacing=dspacing, dim=da.dims).plot(), duplex_peaks_array())

[20]:

Automatic peak finder workflow#

[21]:

wf = BeerModMcStasWorkflow()
wf[Filename[SampleRun]] = mcstas_duplex(16)
da = wf.compute(DetectorTofData[SampleRun])
da = da.transform_coords(('dspacing',), graph=scn.conversion.graph.tof.elastic('tof'),)
ground_truth_peak_positions(da.hist(dspacing=dspacing, dim=da.dims).plot(), duplex_peaks_array())

[21]:

Data that was excluded because tof could not be computed#

The most common reason this happens is overlap between nearby peaks.

[22]:

da = da['bank2'].copy()
da.masks.clear()
da.bins.masks['too_far_from_center'] = ~da.bins.masks.pop('too_far_from_center')

da.bins.concat().hist(two_theta=1000, event_time_offset=1000).plot(norm='log')

[22]: