Neutron Powder Diffraction#

In this tutorial demonstrates how neutron-scattering data can be loaded, visualized, and manipulated with generic functionality from scipp as well as neutron-specific functionality from scippneutron. It focuses on reducing data from the ORNL POWGEN neutron diffractometer.

[1]:

import numpy as np
import scipp as sc
import scippneutron as scn
import scippneutron.data

Loading Nexus files#

Loading Nexus files requires Mantid. See, e.g., Installation on how to install scipp and Mantid with conda.

We start by loading two files: the sample and the vanadium runs.

[2]:

raw_sample = scn.load_with_mantid(
    scn.data.get_path('PG3_4844_event.nxs'),
    load_pulse_times=False,
    mantid_args={'LoadMonitors': True},
)
raw_vanadium = scn.load_with_mantid(
    scn.data.get_path('PG3_4866_event.nxs'), load_pulse_times=False
)

FrameworkManager-[Notice] Welcome to Mantid 6.10.0
FrameworkManager-[Notice] Please cite: http://dx.doi.org/10.1016/j.nima.2014.07.029 and this release: http://dx.doi.org/10.5286/Software/Mantid6.10
CheckMantidVersion-[Notice] A new version of Mantid(6.11.0) is available for download from https://download.mantidproject.org
DownloadInstrument-[Notice] All instrument definitions up to date
Load-[Notice] Load started
Load-[Notice] Load successful, Duration 1.60 seconds
ExtractSpectra-[Notice] ExtractSpectra started
ExtractSpectra-[Notice] ExtractSpectra successful, Duration 0.00 seconds
DeleteWorkspace-[Notice] DeleteWorkspace started
DeleteWorkspace-[Notice] DeleteWorkspace successful, Duration 0.00 seconds
DeleteWorkspace-[Notice] DeleteWorkspace started
DeleteWorkspace-[Notice] DeleteWorkspace successful, Duration 0.00 seconds
DeleteWorkspace-[Notice] DeleteWorkspace started
DeleteWorkspace-[Notice] DeleteWorkspace successful, Duration 0.00 seconds
Load-[Notice] Load started
Load-[Notice] Load successful, Duration 4.20 seconds
DeleteWorkspace-[Notice] DeleteWorkspace started
DeleteWorkspace-[Notice] DeleteWorkspace successful, Duration 0.00 seconds

The optional mantid_args dict is forwarded to the Mantid algorithm used for loading the files – in this case LoadEventNexus – and can be used to control, e.g., which part of a file to load. Here we request loading monitors, which Mantid does not load by default. The resulting dataset looks as follows:

[3]:

raw_sample

[3]:

scipp.DataGroup

(spectrum: 24794, tof: None, time: None)

data
scipp
DataArray
(spectrum: 24794, tof: 1)
DataArrayView
binned data [len=0, len=0, ..., len=0, len=0]
sample
scipp
Variable
()
PyObject
<mantid.api._api.Sample object at 0x7fa61258ec70>
instrument_name
str
()
POWGEN
start_time
scipp
Variable
()
string
𝟙
2011-08-12T15:50:17
end_time
scipp
Variable
()
string
𝟙
2011-08-12T17:22:05
ChopperStatus1
scipp
DataArray
(time: 2)
float64
𝟙
4.0, 4.0
ChopperStatus2
scipp
DataArray
(time: 2)
float64
𝟙
4.0, 4.0
ChopperStatus3
scipp
DataArray
(time: 2)
float64
𝟙
4.0, 4.0
CurrentSP
scipp
DataArray
(time: 2)
float64
𝟙
300.0, 300.0
EnergyRequest
scipp
DataArray
(time: 2)
float64
meV
287.955, 287.955
LKSRampRate
scipp
DataArray
(time: 2)
float64
𝟙
0.0, 0.0
LambdaRequest
scipp
DataArray
(time: 2)
float64
Å
0.533, 0.533
Phase1
scipp
DataArray
(time: 1794)
float64
µs
8166.720, 8165.158, ..., 8163.853, 8163.853
Phase2
scipp
DataArray
(time: 1793)
float64
µs
8335.626, 8334.088, ..., 8332.859, 8332.859
Phase3
scipp
DataArray
(time: 1777)
float64
µs
1.400e+04, 1.400e+04, ..., 1.400e+04, 1.400e+04
PhaseRequest1
scipp
DataArray
(time: 2)
float64
µs
8164.075, 8164.075
PhaseRequest2
scipp
DataArray
(time: 2)
float64
µs
8332.893, 8332.893
PhaseRequest3
scipp
DataArray
(time: 2)
float64
µs
1.400e+04, 1.400e+04
SampleTemp
scipp
DataArray
(time: 467)
float64
𝟙
299.352, 299.446, ..., 300.0, 300.0
Speed1
scipp
DataArray
(time: 2)
float64
Hz
60.0, 60.0
Speed2
scipp
DataArray
(time: 2)
float64
Hz
60.0, 60.0
Speed3
scipp
DataArray
(time: 2)
float64
Hz
60.0, 60.0
SpeedRequest1
scipp
DataArray
(time: 2)
float64
Hz
60.0, 60.0
SpeedRequest2
scipp
DataArray
(time: 2)
float64
Hz
60.0, 60.0
SpeedRequest3
scipp
DataArray
(time: 2)
float64
Hz
60.0, 60.0
TolRequest
scipp
DataArray
(time: 2)
float64
𝟙
20.0, 20.0
currentsample
scipp
DataArray
(time: 2)
float64
𝟙
4.0, 4.0
fernsstatus
scipp
DataArray
(time: 2)
float64
𝟙
3.0, 3.0
frequency
scipp
DataArray
(time: 330473)
float64
Hz
0.0, 60.024, ..., 60.002, 59.999
proton_charge
scipp
DataArray
(time: 330473)
float64
pC
1.484e+07, 1.484e+07, ..., 1.487e+07, 1.484e+07
samplerequest
scipp
DataArray
(time: 2)
float64
𝟙
4.0, 4.0
S1HCenter
scipp
DataArray
(time: 2)
float64
mm
0.0, 0.0
S1HCenterOffset
scipp
DataArray
(time: 2)
float64
𝟙
0.0, 0.0
S1HWidth
scipp
DataArray
(time: 2)
float64
mm
10.0, 10.0
S1VCenter
scipp
DataArray
(time: 2)
float64
mm
5.0, 5.0
S1VCenterOffset
scipp
DataArray
(time: 2)
float64
𝟙
0.0, 0.0
S1VHeight
scipp
DataArray
(time: 2)
float64
mm
30.0, 30.0
commErrs
scipp
DataArray
(time: 2)
float64
𝟙
0.0, 0.0
guide
scipp
DataArray
(time: 2)
float64
mm
-55.463, -55.463
s1b
scipp
DataArray
(time: 2)
float64
mm
20.0, 20.0
s1l
scipp
DataArray
(time: 2)
float64
mm
5.0, 5.0
s1r
scipp
DataArray
(time: 2)
float64
mm
-5.0, -5.0
s1t
scipp
DataArray
(time: 2)
float64
mm
-10.0, -10.0
vGuide
scipp
DataArray
(time: 2)
float64
𝟙
2.0, 2.0
veto_pulse_time
scipp
DataArray
(time: 1)
float64
𝟙
0.0
gd_prtn_chrg
scipp
Variable
()
float64
µAh
1171.953902925
run_start
scipp
Variable
()
string
𝟙
2011-08-12T15:50:17
run_title
scipp
Variable
()
string
𝟙
diamond cw0.533 4.22e12 60Hz [10x30]
file_notes
scipp
Variable
()
string
𝟙
NONE
run_number
scipp
Variable
()
string
𝟙
4844
experiment_identifier
scipp
Variable
()
string
𝟙
IPTS-2767
duration
scipp
Variable
()
float64
s
5508.0
running
scipp
DataArray
(time: 1)
bool
𝟙
True
Filename
scipp
Variable
()
string
𝟙
/home/runner/.cache/scippneutron/5/PG3_4844_event.nxs
monitors
scipp
DataGroup
(tof: 200001)

Extract the actual events.

[4]:

sample = raw_sample['data']
vanadium = raw_vanadium['data']

[5]:

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

[5]:

$../_images/tutorials_powder-diffraction_8_0.svg$

Instrument view#

Scipp provides a simple 3D instrument view inspired by Mantid’s own instrument view, which can be used to take a quick look at the neutron counts on the detector panels in 3D space or using various cylindrical and spherical projections

[6]:

scn.instrument_view(sample.hist())

[6]:

From events to histogram#

Note: We are continuing here with data that has not been normalized to the monitors.

We histogram the event data:

[22]:

dspacing = sc.arange(dim='dspacing', start=0.3, stop=2.0, step=0.001, unit='Angstrom')
hist = sc.Dataset(
    data={
        'sample': dspacing_sample.hist(dspacing=dspacing, dim='dspacing'),
        'vanadium': dspacing_vanadium.hist(dspacing=dspacing, dim='dspacing'),
    }
)
sc.show(hist['spectrum', 0:3]['dspacing', 0:7])

[23]:

sc.plot(hist['sample'].hist(spectrum=500).transpose())

[23]:

$../_images/tutorials_powder-diffraction_38_0.svg$

Summing (focussing) and normalizing#

After conversion to 'dspacing', generic sum and / operations can be used to “focus” and normalize the diffraction data to the vanadium run:

[24]:

summed = sc.sum(hist, 'spectrum')
sc.plot(summed)

[24]:

$../_images/tutorials_powder-diffraction_40_0.svg$

[25]:

normalized = summed['sample'] / summed['vanadium']
sc.plot(normalized)

[25]:

$../_images/tutorials_powder-diffraction_41_0.svg$

Focussing with $\theta$ dependence in event-mode#

The approach used above combines reflections from all crystallographic planes and is therefore of limited use. We can use bin to focus each of multiple groups of spectra into a distinct output spectrum. Here we define groups based on a range of scattering angles – a simple $\theta$-dependent binning. This also demonstrates how we can postpone histogramming until after the focussing step.

[26]:

two_theta = sc.linspace(dim='two_theta', unit='rad', start=0.0, stop=np.pi, num=16)

focussed_sample = dspacing_sample.bin(two_theta=two_theta)
focussed_vanadium = dspacing_vanadium.bin(two_theta=two_theta)
norm = focussed_vanadium.hist(dspacing=dspacing)

Similarly as when normalizing by monitor, we have to drop the variances of the normalization term. Otherwise, the normalization would be broadcast to focussed_sample which would introduce correlations. In practice, we have to make sure that the uncertainties of the vanadium measurement can be neglected!

[27]:

norm = sc.values(norm)

[28]:

focussed_sample.bins /= sc.lookup(func=norm, dim='dspacing')
normalized = focussed_sample.hist(dspacing=dspacing)

The normalized output looks as follows:

[29]:

normalized.plot(vmin=sc.scalar(0), vmax=sc.scalar(0.4))

[29]:

$../_images/tutorials_powder-diffraction_48_0.svg$

As a bonus, we can use slicing and a dict-comprehension to quickly create of plot comparing the spectra for different scattering angle bins:

[30]:

# compute centers of theta bins
angle = normalized.coords['two_theta'].values
angle = 0.5 * (angle[1:] + angle[:-1])
results = {
    f'{round(angle[group], 3)} rad': normalized['dspacing', 300:500]['two_theta', group]
    for group in range(2, 6)
}
sc.plot(results)

[30]:

$../_images/tutorials_powder-diffraction_50_0.svg$

Neutron Powder Diffraction#

Loading Nexus files#

Instrument view#

Plot against scattering angle \(\theta\) using `bin`#

Coordinate transformation#

Neutron monitors#

From events to histogram#

Summing (focussing) and normalizing#

Focussing with \(\theta\) dependence in event-mode#

Neutron Powder Diffraction#

Loading Nexus files#

Instrument view#

Plot against scattering angle \(\theta\) using bin#

Coordinate transformation#

Neutron monitors#

From events to histogram#

Summing (focussing) and normalizing#

Focussing with \(\theta\) dependence in event-mode#

Plot against scattering angle \(\theta\) using `bin`#