Powder diffraction data reduction#
This notebook will guide you through the data reduction for the powder diffraction experiment that you simulated with McStas yesterday.
The following is a basic outline of what this notebook will cover:
Loading the NeXus files that contain the data
Inspect/visualize the data contents
How to convert the raw event
time-of-arrival
coordinate to something more useful (\(\lambda\), \(d\), …)Improve the estimate of neutron
time-of-flight
Save the reduced data to disk, in a state that is ready for analysis/fitting.
import numpy as np
import scipp as sc
import plopp as pp
import utils
Process the run with a Si sample#
Load the NeXus file data#
folder = "../3-mcstas/powder_Si_high_resolution"
sample_si = utils.load_powder(folder)
The first way to inspect the data is to view the HTML representation of the loaded object.
Try to explore what is inside the data, and familiarize yourself with the different sections (Dimensions
, Coordinates
, Data
).
sample_si
- events: 29174
- position(events)vector3m[2.16682882e+00 1.12500000e-01 1.63387609e+02], [ 3.19116715 -0.2625 159.20148402], ..., [-2.23803651e+00 -1.12500000e-01 1.57948054e+02], [ -1.8157064 0.5625 157.64680845]
Values:
array([[ 2.16682882e+00, 1.12500000e-01, 1.63387609e+02], [ 3.19116715e+00, -2.62500000e-01, 1.59201484e+02], [ 2.96004738e+00, -3.75000000e-02, 1.62506651e+02], ..., [-1.03789551e+00, 1.87500000e-01, 1.57296430e+02], [-2.23803651e+00, -1.12500000e-01, 1.57948054e+02], [-1.81570640e+00, 5.62500000e-01, 1.57646808e+02]], shape=(29174, 3)) - sample_position()vector3m[ 0. 0. 160.639]
Values:
array([ 0. , 0. , 160.639]) - sim_source_time(events)float64s0.002, 0.002, ..., 0.001, 0.001
Values:
array([0.00225332, 0.00150513, 0.00030688, ..., 0.00134982, 0.0006419 , 0.00097281], shape=(29174,)) - sim_speed(events)float64m/s1943.875, 1685.779, ..., 1343.720, 1521.287
Values:
array([1943.87495082, 1685.77870347, 1288.80471417, ..., 1609.50209365, 1343.72030251, 1521.28722408], shape=(29174,)) - sim_wavelength(events)float64Å2.035, 2.347, ..., 2.944, 2.600
Values:
array([2.03512783, 2.34671016, 3.0695372 , ..., 2.45792412, 2.94409038, 2.60045174], shape=(29174,)) - source_position()vector3m[0. 0. 0.]
Values:
array([0., 0., 0.]) - time_origin()float64s0.0
Values:
array(0.) - toa(events)float64s0.087, 0.099, ..., 0.123, 0.109
Values:
array([0.08667913, 0.09890457, 0.1277125 , ..., 0.10331841, 0.12277178, 0.10892322], shape=(29174,)) - x(events)float64m2.167, 3.191, ..., -2.238, -1.816
Values:
array([ 2.16682882, 3.19116715, 2.96004738, ..., -1.03789551, -2.23803651, -1.8157064 ], shape=(29174,)) - y(events)float64m0.113, -0.263, ..., -0.113, 0.562
Values:
array([ 0.1125, -0.2625, -0.0375, ..., 0.1875, -0.1125, 0.5625], shape=(29174,)) - z(events)float64m163.388, 159.201, ..., 157.948, 157.647
Values:
array([163.38760926, 159.20148402, 162.50665081, ..., 157.29643019, 157.94805359, 157.64680845], shape=(29174,))
- (events)float64counts593.365, 2.753, ..., 1.474e+04, 1072.312σ = 593.365, 2.753, ..., 1.474e+04, 1072.312
Values:
array([5.93364766e+02, 2.75336983e+00, 3.71286602e+02, ..., 1.55654788e+02, 1.47423092e+04, 1.07231238e+03], shape=(29174,))
Variances (σ²):
array([3.52081746e+05, 7.58104545e+00, 1.37853741e+05, ..., 2.42284130e+04, 2.17335681e+08, 1.14985384e+06], shape=(29174,))
Visualize the data#
Here is a 2D visualization of the neutron counts, histogrammed along the x
and toa
dimensions:
sample_si.hist(x=200, toa=200).plot(norm="log", vmin=1.0e-2)
We can also visualize the events in 3D, which show the shape of the detector panels:
pp.scatter3d(sample_si[::10], pos='position', cbar=True, vmin=0.01, pixel_size=0.07)
Coordinate transformations#
The first step of this data reduction workflow is to convert the raw event coordinates (position, time-of-flight) to something physically meaningful such as wavelength (\(\lambda\)) or d-spacing (\(d\)).
Scipp has a dedicated method for this called transform_coords
(see docs here).
We begin with a standard graph which describes how to compute e.g. the wavelength from the other coordinates in the raw data.
from scippneutron.conversion.graph.beamline import beamline
from scippneutron.conversion.graph.tof import kinematic
graph = {**beamline(scatter=True), **kinematic("tof")}
def compute_tof(toa, time_origin):
return toa - time_origin
graph['tof'] = compute_tof
sc.show_graph(graph, simplified=True)
To compute the wavelength of all the events, we simply call transform_coords
on our loaded data,
requesting the name of the coordinate we want in the output ("wavelength"
),
as well as providing it the graph to be used to compute it (i.e. the one we defined just above).
This yields
si_wav = sample_si.transform_coords(["wavelength", "two_theta"], graph=graph)
si_wav
- events: 29174
- L1()float64m160.63899999999998
Values:
array(160.639) - L2(events)float64m3.502, 3.510, ..., 3.502, 3.545
Values:
array([3.50180757, 3.50982995, 3.50020089, ..., 3.50501872, 3.50180757, 3.54491273], shape=(29174,)) - Ltotal(events)float64m164.141, 164.149, ..., 164.141, 164.184
Values:
array([164.14080757, 164.14882995, 164.13920089, ..., 164.14401872, 164.14080757, 164.18391273], shape=(29174,)) - incident_beam()vector3m[ 0. 0. 160.639]
Values:
array([ 0. , 0. , 160.639]) - position(events)vector3m[2.16682882e+00 1.12500000e-01 1.63387609e+02], [ 3.19116715 -0.2625 159.20148402], ..., [-2.23803651e+00 -1.12500000e-01 1.57948054e+02], [ -1.8157064 0.5625 157.64680845]
Values:
array([[ 2.16682882e+00, 1.12500000e-01, 1.63387609e+02], [ 3.19116715e+00, -2.62500000e-01, 1.59201484e+02], [ 2.96004738e+00, -3.75000000e-02, 1.62506651e+02], ..., [-1.03789551e+00, 1.87500000e-01, 1.57296430e+02], [-2.23803651e+00, -1.12500000e-01, 1.57948054e+02], [-1.81570640e+00, 5.62500000e-01, 1.57646808e+02]], shape=(29174, 3)) - sample_position()vector3m[ 0. 0. 160.639]
Values:
array([ 0. , 0. , 160.639]) - scattered_beam(events)vector3m[2.16682882 0.1125 2.74860926], [ 3.19116715 -0.2625 -1.43751598], ..., [-2.23803651 -0.1125 -2.69094641], [-1.8157064 0.5625 -2.99219155]
Values:
array([[ 2.16682882, 0.1125 , 2.74860926], [ 3.19116715, -0.2625 , -1.43751598], [ 2.96004738, -0.0375 , 1.86765081], ..., [-1.03789551, 0.1875 , -3.34256981], [-2.23803651, -0.1125 , -2.69094641], [-1.8157064 , 0.5625 , -2.99219155]], shape=(29174, 3)) - sim_source_time(events)float64s0.002, 0.002, ..., 0.001, 0.001
Values:
array([0.00225332, 0.00150513, 0.00030688, ..., 0.00134982, 0.0006419 , 0.00097281], shape=(29174,)) - sim_speed(events)float64m/s1943.875, 1685.779, ..., 1343.720, 1521.287
Values:
array([1943.87495082, 1685.77870347, 1288.80471417, ..., 1609.50209365, 1343.72030251, 1521.28722408], shape=(29174,)) - sim_wavelength(events)float64Å2.035, 2.347, ..., 2.944, 2.600
Values:
array([2.03512783, 2.34671016, 3.0695372 , ..., 2.45792412, 2.94409038, 2.60045174], shape=(29174,)) - source_position()vector3m[0. 0. 0.]
Values:
array([0., 0., 0.]) - time_origin()float64s0.0
Values:
array(0.) - toa(events)float64s0.087, 0.099, ..., 0.123, 0.109
Values:
array([0.08667913, 0.09890457, 0.1277125 , ..., 0.10331841, 0.12277178, 0.10892322], shape=(29174,)) - tof(events)float64s0.087, 0.099, ..., 0.123, 0.109
Values:
array([0.08667913, 0.09890457, 0.1277125 , ..., 0.10331841, 0.12277178, 0.10892322], shape=(29174,)) - two_theta(events)float64rad0.668, 1.993, ..., 2.447, 2.576
Values:
array([0.66824294, 1.99277741, 1.00796385, ..., 2.83594572, 2.44720387, 2.5756443 ], shape=(29174,)) - wavelength(events)float64Å2.089, 2.384, ..., 2.959, 2.625
Values:
array([2.08909412, 2.3836285 , 3.07808865, ..., 2.49007626, 2.95897989, 2.62451988], shape=(29174,)) - x(events)float64m2.167, 3.191, ..., -2.238, -1.816
Values:
array([ 2.16682882, 3.19116715, 2.96004738, ..., -1.03789551, -2.23803651, -1.8157064 ], shape=(29174,)) - y(events)float64m0.113, -0.263, ..., -0.113, 0.562
Values:
array([ 0.1125, -0.2625, -0.0375, ..., 0.1875, -0.1125, 0.5625], shape=(29174,)) - z(events)float64m163.388, 159.201, ..., 157.948, 157.647
Values:
array([163.38760926, 159.20148402, 162.50665081, ..., 157.29643019, 157.94805359, 157.64680845], shape=(29174,))
- (events)float64counts593.365, 2.753, ..., 1.474e+04, 1072.312σ = 593.365, 2.753, ..., 1.474e+04, 1072.312
Values:
array([5.93364766e+02, 2.75336983e+00, 3.71286602e+02, ..., 1.55654788e+02, 1.47423092e+04, 1.07231238e+03], shape=(29174,))
Variances (σ²):
array([3.52081746e+05, 7.58104545e+00, 1.37853741e+05, ..., 2.42284130e+04, 2.17335681e+08, 1.14985384e+06], shape=(29174,))
The result has a wavelength
coordinate. We can also plot the result:
si_wav.hist(two_theta=120, wavelength=200).plot(norm='log', vmin=1.0e-3)
Exercise 1: convert raw data to d-spacing#
Instead of wavelength as in the example above, the task is now to convert the raw data to interplanar lattice spacing \(d\).
The transformation graph is missing the computation for \(d\) so you will have to add it in yourself. As a reminder, \(d\) is computed as follows
You have to:
create a function that computes \(d\)
add it to the graph with name “dspacing”
call
transform_coords
using the new graph
Note that the graph already contains the necessary components to compute the scattering angle \(2 \theta\) (called two_theta
in code).
Solution:
Show code cell content
def compute_d(two_theta, wavelength):
return wavelength / (2 * sc.sin(two_theta / 2))
graph["dspacing"] = compute_d
# Show the updated graph
display(sc.show_graph(graph, simplified=True))
# Run the coordinate transformation
si_dspacing = sample_si.transform_coords("dspacing", graph=graph)
si_dspacing
- events: 29174
- L1()float64m160.63899999999998
Values:
array(160.639) - L2(events)float64m3.502, 3.510, ..., 3.502, 3.545
Values:
array([3.50180757, 3.50982995, 3.50020089, ..., 3.50501872, 3.50180757, 3.54491273], shape=(29174,)) - Ltotal(events)float64m164.141, 164.149, ..., 164.141, 164.184
Values:
array([164.14080757, 164.14882995, 164.13920089, ..., 164.14401872, 164.14080757, 164.18391273], shape=(29174,)) - dspacing(events)float64Å3.185, 1.420, ..., 1.573, 1.367
Values:
array([3.18518365, 1.41964733, 3.18698012, ..., 1.25971991, 1.57337136, 1.36661108], shape=(29174,)) - incident_beam()vector3m[ 0. 0. 160.639]
Values:
array([ 0. , 0. , 160.639]) - position(events)vector3m[2.16682882e+00 1.12500000e-01 1.63387609e+02], [ 3.19116715 -0.2625 159.20148402], ..., [-2.23803651e+00 -1.12500000e-01 1.57948054e+02], [ -1.8157064 0.5625 157.64680845]
Values:
array([[ 2.16682882e+00, 1.12500000e-01, 1.63387609e+02], [ 3.19116715e+00, -2.62500000e-01, 1.59201484e+02], [ 2.96004738e+00, -3.75000000e-02, 1.62506651e+02], ..., [-1.03789551e+00, 1.87500000e-01, 1.57296430e+02], [-2.23803651e+00, -1.12500000e-01, 1.57948054e+02], [-1.81570640e+00, 5.62500000e-01, 1.57646808e+02]], shape=(29174, 3)) - sample_position()vector3m[ 0. 0. 160.639]
Values:
array([ 0. , 0. , 160.639]) - scattered_beam(events)vector3m[2.16682882 0.1125 2.74860926], [ 3.19116715 -0.2625 -1.43751598], ..., [-2.23803651 -0.1125 -2.69094641], [-1.8157064 0.5625 -2.99219155]
Values:
array([[ 2.16682882, 0.1125 , 2.74860926], [ 3.19116715, -0.2625 , -1.43751598], [ 2.96004738, -0.0375 , 1.86765081], ..., [-1.03789551, 0.1875 , -3.34256981], [-2.23803651, -0.1125 , -2.69094641], [-1.8157064 , 0.5625 , -2.99219155]], shape=(29174, 3)) - sim_source_time(events)float64s0.002, 0.002, ..., 0.001, 0.001
Values:
array([0.00225332, 0.00150513, 0.00030688, ..., 0.00134982, 0.0006419 , 0.00097281], shape=(29174,)) - sim_speed(events)float64m/s1943.875, 1685.779, ..., 1343.720, 1521.287
Values:
array([1943.87495082, 1685.77870347, 1288.80471417, ..., 1609.50209365, 1343.72030251, 1521.28722408], shape=(29174,)) - sim_wavelength(events)float64Å2.035, 2.347, ..., 2.944, 2.600
Values:
array([2.03512783, 2.34671016, 3.0695372 , ..., 2.45792412, 2.94409038, 2.60045174], shape=(29174,)) - source_position()vector3m[0. 0. 0.]
Values:
array([0., 0., 0.]) - time_origin()float64s0.0
Values:
array(0.) - toa(events)float64s0.087, 0.099, ..., 0.123, 0.109
Values:
array([0.08667913, 0.09890457, 0.1277125 , ..., 0.10331841, 0.12277178, 0.10892322], shape=(29174,)) - tof(events)float64s0.087, 0.099, ..., 0.123, 0.109
Values:
array([0.08667913, 0.09890457, 0.1277125 , ..., 0.10331841, 0.12277178, 0.10892322], shape=(29174,)) - two_theta(events)float64rad0.668, 1.993, ..., 2.447, 2.576
Values:
array([0.66824294, 1.99277741, 1.00796385, ..., 2.83594572, 2.44720387, 2.5756443 ], shape=(29174,)) - wavelength(events)float64Å2.089, 2.384, ..., 2.959, 2.625
Values:
array([2.08909412, 2.3836285 , 3.07808865, ..., 2.49007626, 2.95897989, 2.62451988], shape=(29174,)) - x(events)float64m2.167, 3.191, ..., -2.238, -1.816
Values:
array([ 2.16682882, 3.19116715, 2.96004738, ..., -1.03789551, -2.23803651, -1.8157064 ], shape=(29174,)) - y(events)float64m0.113, -0.263, ..., -0.113, 0.562
Values:
array([ 0.1125, -0.2625, -0.0375, ..., 0.1875, -0.1125, 0.5625], shape=(29174,)) - z(events)float64m163.388, 159.201, ..., 157.948, 157.647
Values:
array([163.38760926, 159.20148402, 162.50665081, ..., 157.29643019, 157.94805359, 157.64680845], shape=(29174,))
- (events)float64counts593.365, 2.753, ..., 1.474e+04, 1072.312σ = 593.365, 2.753, ..., 1.474e+04, 1072.312
Values:
array([5.93364766e+02, 2.75336983e+00, 3.71286602e+02, ..., 1.55654788e+02, 1.47423092e+04, 1.07231238e+03], shape=(29174,))
Variances (σ²):
array([3.52081746e+05, 7.58104545e+00, 1.37853741e+05, ..., 2.42284130e+04, 2.17335681e+08, 1.14985384e+06], shape=(29174,))
Histogram the data in d-spacing#
The final step in processing the sample run is to histogram the data into \(d\) bins.
dbins = sc.linspace('dspacing', 0.5, 2.0, 201, unit='angstrom')
si_dspacing_hist = si_dspacing.hist(two_theta=120, dspacing=dbins)
si_dspacing_hist.plot(norm='log', vmin=1.0e-3)
si_dspacing_hist.sum('two_theta').plot()
Exercise 2: Correct the time-of-flight#
Time-of-flight at a long-pulsed source#
The time-of-flight of a neutron is typically not known from experimental measurements.
The only two parameters that we are able to record are usually the time at which the pulse was emitted (named time_origin
in our dataset),
and the timestamp when a neutron was detected by a detector pixel (named toa
, standing for time-of-arrival).
The difference between these two times is commonly used to compute an approximate time-of-flight. This approximation performs well at neutron sources with a short pulse, where it is satisfactory to assume that all neutrons left the source at the same time.
At a long-pulsed source (lasting a few milliseconds) such as ESS, with high-precision instruments, this approximation is no longer valid.
This is in fact the reason why in the \(2\theta\) vs d-spacing plot, the bright lines appear curved while they should in reality be vertical.
In the following exercise, we will add a correction to our time-of-flight computation to improve our results.
Finding a new time-distance origin#
One way to look at the problem is to visualize the paths of each neutron that made it to the detector, on a time-distance diagram.
The neutrons begin at the source (bottom left), travel along the guide, go through the chopper, and finally arrive at the detectors (top).
from powder_utils import time_distance_diagram
time_distance_diagram(si_wav)

As the figure above illustrates,
the current time-of-flight computation (tof = toa - time_origin
) assumes that all neutrons were born at \(t = 0\) (black dot).
This is clearly not the case from looking at the neutron tracks.
Exercise: find a new origin#
Inspect the time-distance diagram and try to find a point which would be a better representation of a common origin for all tracks.
Set the origin as coordinates on the original data (hint: you will need to update the
source_position
andtime_origin
coordinates).Compute wavelength/d-spacing again.
Solution:
Show code cell content
better = sample_si.copy(deep=False)
# Read off the plot that the lines cross as they go through the chopper:
# time=5.53, distance=6.5 (approximately)
better.coords['source_position'] = sc.vector([0, 0, 6.5], unit='m')
better.coords['time_origin'] = sc.scalar(5.53, unit='ms').to(unit='s')
better_dspacing = better.transform_coords("dspacing", graph=graph)
better_dspacing
- events: 29174
- L1()float64m154.13899999999998
Values:
array(154.139) - L2(events)float64m3.502, 3.510, ..., 3.502, 3.545
Values:
array([3.50180757, 3.50982995, 3.50020089, ..., 3.50501872, 3.50180757, 3.54491273], shape=(29174,)) - Ltotal(events)float64m157.641, 157.649, ..., 157.641, 157.684
Values:
array([157.64080757, 157.64882995, 157.63920089, ..., 157.64401872, 157.64080757, 157.68391273], shape=(29174,)) - dspacing(events)float64Å3.105, 1.396, ..., 1.564, 1.351
Values:
array([3.10492932, 1.39553189, 3.17470259, ..., 1.2414557 , 1.56445475, 1.35070261], shape=(29174,)) - incident_beam()vector3m[ 0. 0. 154.139]
Values:
array([ 0. , 0. , 154.139]) - position(events)vector3m[2.16682882e+00 1.12500000e-01 1.63387609e+02], [ 3.19116715 -0.2625 159.20148402], ..., [-2.23803651e+00 -1.12500000e-01 1.57948054e+02], [ -1.8157064 0.5625 157.64680845]
Values:
array([[ 2.16682882e+00, 1.12500000e-01, 1.63387609e+02], [ 3.19116715e+00, -2.62500000e-01, 1.59201484e+02], [ 2.96004738e+00, -3.75000000e-02, 1.62506651e+02], ..., [-1.03789551e+00, 1.87500000e-01, 1.57296430e+02], [-2.23803651e+00, -1.12500000e-01, 1.57948054e+02], [-1.81570640e+00, 5.62500000e-01, 1.57646808e+02]], shape=(29174, 3)) - sample_position()vector3m[ 0. 0. 160.639]
Values:
array([ 0. , 0. , 160.639]) - scattered_beam(events)vector3m[2.16682882 0.1125 2.74860926], [ 3.19116715 -0.2625 -1.43751598], ..., [-2.23803651 -0.1125 -2.69094641], [-1.8157064 0.5625 -2.99219155]
Values:
array([[ 2.16682882, 0.1125 , 2.74860926], [ 3.19116715, -0.2625 , -1.43751598], [ 2.96004738, -0.0375 , 1.86765081], ..., [-1.03789551, 0.1875 , -3.34256981], [-2.23803651, -0.1125 , -2.69094641], [-1.8157064 , 0.5625 , -2.99219155]], shape=(29174, 3)) - sim_source_time(events)float64s0.002, 0.002, ..., 0.001, 0.001
Values:
array([0.00225332, 0.00150513, 0.00030688, ..., 0.00134982, 0.0006419 , 0.00097281], shape=(29174,)) - sim_speed(events)float64m/s1943.875, 1685.779, ..., 1343.720, 1521.287
Values:
array([1943.87495082, 1685.77870347, 1288.80471417, ..., 1609.50209365, 1343.72030251, 1521.28722408], shape=(29174,)) - sim_wavelength(events)float64Å2.035, 2.347, ..., 2.944, 2.600
Values:
array([2.03512783, 2.34671016, 3.0695372 , ..., 2.45792412, 2.94409038, 2.60045174], shape=(29174,)) - source_position()vector3m[0. 0. 6.5]
Values:
array([0. , 0. , 6.5]) - time_origin()float64s0.00553
Values:
array(0.00553) - toa(events)float64s0.087, 0.099, ..., 0.123, 0.109
Values:
array([0.08667913, 0.09890457, 0.1277125 , ..., 0.10331841, 0.12277178, 0.10892322], shape=(29174,)) - tof(events)float64s0.081, 0.093, ..., 0.117, 0.103
Values:
array([0.08114913, 0.09337457, 0.1221825 , ..., 0.09778841, 0.11724178, 0.10339322], shape=(29174,)) - two_theta(events)float64rad0.668, 1.993, ..., 2.447, 2.576
Values:
array([0.66824294, 1.99277741, 1.00796385, ..., 2.83594572, 2.44720387, 2.5756443 ], shape=(29174,)) - wavelength(events)float64Å2.036, 2.343, ..., 2.942, 2.594
Values:
array([2.03645701, 2.34313799, 3.06623061, ..., 2.45397358, 2.94221076, 2.59396833], shape=(29174,)) - x(events)float64m2.167, 3.191, ..., -2.238, -1.816
Values:
array([ 2.16682882, 3.19116715, 2.96004738, ..., -1.03789551, -2.23803651, -1.8157064 ], shape=(29174,)) - y(events)float64m0.113, -0.263, ..., -0.113, 0.562
Values:
array([ 0.1125, -0.2625, -0.0375, ..., 0.1875, -0.1125, 0.5625], shape=(29174,)) - z(events)float64m163.388, 159.201, ..., 157.948, 157.647
Values:
array([163.38760926, 159.20148402, 162.50665081, ..., 157.29643019, 157.94805359, 157.64680845], shape=(29174,))
- (events)float64counts593.365, 2.753, ..., 1.474e+04, 1072.312σ = 593.365, 2.753, ..., 1.474e+04, 1072.312
Values:
array([5.93364766e+02, 2.75336983e+00, 3.71286602e+02, ..., 1.55654788e+02, 1.47423092e+04, 1.07231238e+03], shape=(29174,))
Variances (σ²):
array([3.52081746e+05, 7.58104545e+00, 1.37853741e+05, ..., 2.42284130e+04, 2.17335681e+08, 1.14985384e+06], shape=(29174,))
better_dspacing.hist(two_theta=120, dspacing=dbins).plot(norm='log', vmin=1.0e-3)
pp.plot({'original': si_dspacing_hist.sum('two_theta'),
'improved': better_dspacing.hist(two_theta=120, dspacing=dbins).sum('two_theta')})
Process the Vanadium run#
folder = "../3-mcstas/powder_vanadium"
sample_van = utils.load_powder(folder)
# Read off the plot that the lines cross as they go through the chopper:
# time=5.53, distance=6.5 (approximately)
sample_van.coords['source_position'] = sc.vector([0, 0, 6.5], unit='m')
sample_van.coords['time_origin'] = sc.scalar(5.53, unit='ms').to(unit='s')
van_dspacing = sample_van.transform_coords("dspacing", graph=graph)
van_dspacing.hist(two_theta=120, dspacing=dbins).plot(norm='log', vmin=1.0e-3)
Normalize by Vanadium#
num = better_dspacing.hist(two_theta=120, dspacing=dbins)
den = van_dspacing.hist(two_theta=num.coords['two_theta'], dspacing=dbins)
normed = (num.sum('two_theta') / den.sum('two_theta'))
normed.plot()
Save to disk#
average_l = better_dspacing.coords["Ltotal"].mean()
average_two_theta = better_dspacing.coords["two_theta"].mean()
difc = sc.to_unit(
2
* sc.constants.m_n
/ sc.constants.h
* average_l
* sc.sin(0.5 * average_two_theta),
unit='us / angstrom',
)
difc
- ()float6410N/mW67376.59706577176
Values:
array(67376.59706577)
from scippneutron.io import save_xye
result_si = normed.copy(deep=False)
result_si.coords["tof"] = (sc.midpoints(result_si.coords["dspacing"]) * difc).to(unit='us')
save_xye("powder_reduced_Si.xye", result_si.rename_dims(dspacing='tof'),
header=f"DIFC = {difc.to(unit='us/angstrom').value} [µ/Å]\ntof [µs] Y [counts] E [counts]")
Process and normalize unknown sample#
folder = "../3-mcstas/powder_sample_2"
unknown = utils.load_powder(folder)
# Read off the plot that the lines cross as they go through the chopper:
# time=5.53, distance=6.5 (approximately)
unknown.coords['source_position'] = sc.vector([0, 0, 6.5], unit='m')
unknown.coords['time_origin'] = sc.scalar(5.53, unit='ms').to(unit='s')
unknown_dspacing = unknown.transform_coords("dspacing", graph=graph)
unknown_dspacing.hist(two_theta=120, dspacing=dbins).plot(norm='log', vmin=1.0e-3)
# Normalize
unknown_num = unknown_dspacing.hist(two_theta=den.coords['two_theta'], dspacing=dbins)
unknown_normed = (unknown_num.sum('two_theta') / den.sum('two_theta'))
unknown_normed.plot()
Save to disk#
result_unknown = unknown_normed.copy(deep=False)
result_unknown.coords["tof"] = (sc.midpoints(result_unknown.coords["dspacing"]) * difc).to(unit='us')
save_xye("powder_reduced_unknown.xye", result_unknown.rename_dims(dspacing='tof'),
header=f"DIFC = {difc.to(unit='us/angstrom').value} [µ/Å]\ntof [µs] Y [counts] E [counts]")