Data reduction for POWGEN#

This notebook shows a basic reduction workflow for powder diffraction for the SNS POWGEN instrument. It serves mainly to develop and present routines for powder diffraction and will eventually be removed in favor of a workflow for DREAM at ESS.

Note that we load functions from external modules. These modules will be removed when their ESS counterparts exist.

[1]:

import scipp as sc
import scippneutron as scn
import plopp as pp

import ess
from ess.diffraction import powder
from ess import diffraction
from ess.external import powgen

Initialize the logging system to write logs to a file powgen.log. This also displays a widget which shows log messages emitted by the functions that we call.

[2]:

ess.logging.configure_workflow('powgen_reduction', filename='powgen.log')

[2]:

<Logger scipp.ess.powgen_reduction (INFO)>

Inspect the raw data#

We can plot the data array to get an idea of its contents.

[6]:

sample.hist(spectrum=500, tof=400).plot()

[6]:

../../../_images/instruments_external_powgen_powgen_reduction_13_0.svg

We can see how that data maps onto the detector by using POWGEN’s instrument view.

[7]:

scn.instrument_view(sample.hist())

[7]:

Filter out invalid events#

The file contains events that cannot be associated with a specific pulse. We can get a range of valid time-of-flight values from the instrument characterization file associated with the current run. There is currently no mechanism in scippneutron or ess to load such a file as it is not clear if ESS will use this approach. The values used below are taken from PG3_characterization_2011_08_31-HR.txt which is part of the sample files of Mantid. See, e.g., PowderDiffractionReduction.

We remove all events that have a time-of-flight value outside the valid range:

[8]:

sample = sample.bin(tof=sc.array(dims=['tof'], values=[0.0, 16666.67], unit='us'))

Normalize by proton charge#

Next, we normalize the data by the proton charge.

[9]:

sample /= sample.coords['gd_prtn_chrg']

We can check the unit of the event weight to see that the data was indeed divided by a charge.

[10]:

sample.data.values[0].unit

[10]:

counts/uAh

Inspect d-spacing#

We need to histogram the events in order to normalize our sample data by vanadium. For consistency, we use these bin edges for both vanadium and the sample data.

[22]:

d = vana_dspacing.coords['dspacing']
dspacing_edges = sc.linspace('dspacing', d.min().value, d.max().value, 200, unit=d.unit)

All spectra combined#

We start simple by combining all spectra using data.bins.concat('spectrum'). Then, we can normalize the same data by vanadium to get a d-spacing distribution.

Note that because we removed the variances on the Vanadium data, after the following cell, the standard deviations on the result are underestimated.

[23]:

all_spectra = diffraction.normalize_by_vanadium(
    sample_dspacing.bins.concat('spectrum'),
    vanadium=vana_dspacing.bins.concat('spectrum'),
    edges=dspacing_edges,
)

[24]:

all_spectra.hist(dspacing=dspacing_edges).plot()

[24]:

../../../_images/instruments_external_powgen_powgen_reduction_51_0.svg

Group into \(2\theta\) bins#

For a better resolution, we now group the sample and vanadium data into a number of bins in the scattering angle \(2\theta\) (see here) and normalize each individually.

[25]:

two_theta = sc.linspace(dim='two_theta', unit='deg', start=25.0, stop=90.0, num=16)
sample_by_two_theta = diffraction.group_by_two_theta(sample_dspacing, edges=two_theta)
vana_by_two_theta = diffraction.group_by_two_theta(vana_dspacing, edges=two_theta)

[26]:

normalized = diffraction.normalize_by_vanadium(
    sample_by_two_theta, vanadium=vana_by_two_theta, edges=dspacing_edges
)

Histogram the results in order to get a useful binning in the following plots.

[27]:

normalized = normalized.hist(dspacing=dspacing_edges)

Now we can inspect the d-spacing distribution as a function of \(2\theta\).

[28]:

normalized.plot()

[28]:

../../../_images/instruments_external_powgen_powgen_reduction_59_0.svg

In order to get 1-dimensional plots, we can select some ranges of scattering angles.

[29]:

angle = sc.midpoints(normalized.coords['two_theta']).values
results = {
    f'{round(angle[group], 3)} rad': normalized['two_theta', group]
    for group in range(2, 6)
}
sc.plot(results)

[29]:

../../../_images/instruments_external_powgen_powgen_reduction_61_0.svg

Or interactively by plotting with a 1d projection.

[30]:

%matplotlib widget
pp.superplot(normalized)

[30]:

Data reduction for POWGEN

Contents

Data reduction for POWGEN#

Load data#

Inspect the raw data#

Filter out invalid events#

Normalize by proton charge#

Compute d-spacing#

Vanadium correction#

Removing the variances of the Vanadium data#

Conversion to d-spacing#

Inspect d-spacing#

All spectra combined#

Group into \(2\theta\) bins#