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 ess
from ess.diffraction.external import load_calibration
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.

[5]:

sample.plot()

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

[6]:

powgen.instrument_view(sample)

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:

[7]:

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

Normalize by proton charge¶

Next, we normalize the data by the proton charge.

[8]:

sample /= sample.attrs['gd_prtn_chrg']

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

[9]:

sample.data.values[0].unit

[9]:

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.

[18]:

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.

[19]:

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

[20]:

sc.histogram(all_spectra, bins=dspacing_edges).plot()

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.

[21]:

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)

[22]:

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.

[23]:

normalized = sc.histogram(normalized, bins=dspacing_edges)

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

[24]:

normalized.plot()

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

[25]:

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)

Or interactively by plotting with a 1d projection.

[26]:

normalized.plot(projection='1d')

Data reduction for POWGEN

Contents