SANS data reduction

SANS data reduction#

This notebook will guide you through the data reduction for the SANS 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 time-of-flight coordinate to something more useful (\(\lambda\), \(Q\), …)

  • Normalize to a flat field run

  • Write the results to file

  • Process more than one pulse

import numpy as np
import scipp as sc
import plopp as pp
import utils

Process the run with a sample#

Load the NeXus file data#

folder = "../3-mcstas/SANS_with_sample_1_pulse"

sample = utils.load_sans(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
Show/Hide data repr Show/Hide attributes
scipp.DataArray (423.36 MB)
    • event: 7927147
    • position
      (event)
      vector3
      m
      [0. 0.11814271 3. ], [0. 0.11937109 3. ], ..., [0. 0.19366142 3. ], [0. 0.19253724 3. ]
      Values:
      array([[0.        , 0.11814271, 3.        ],
       [0.        , 0.11937109, 3.        ],
       [0.        , 0.11951049, 3.        ],
       ...,
       [0.        , 0.20169624, 3.        ],
       [0.        , 0.19366142, 3.        ],
       [0.        , 0.19253724, 3.        ]])
    • sample_position
      ()
      vector3
      m
      [0. 0. 0.]
      Values:
      array([0., 0., 0.])
    • source_position
      ()
      vector3
      m
      [ 0. 0. -150.]
      Values:
      array([   0.,    0., -150.])
    • tof
      (event)
      float64
      ms
      253.831, 253.838, ..., 235.612, 235.615
      Values:
      array([253.83145495, 253.83804838, 253.83940871, ..., 257.95434162,
       235.61182539, 235.61501943])
    • x
      (event)
      float64
      m
      0.0, 0.0, ..., 0.0, 0.0
      Values:
      array([0., 0., 0., ..., 0., 0., 0.])
    • y
      (event)
      float64
      m
      0.118, 0.119, ..., 0.194, 0.193
      Values:
      array([0.11814271, 0.11937109, 0.11951049, ..., 0.20169624, 0.19366142,
       0.19253724])
    • (event)
      float64
      counts
      7.985e-09, 6.581e-13, ..., 4.828e-14, 4.341e-18
      Values:
      array([7.98463800e-09, 6.58130913e-13, 5.42462035e-17, ...,
       3.67754741e-28, 4.82836503e-14, 4.34070359e-18])

Visualize the data#

Here is a 2D visualization of the neutron counts, histogrammed along the tof and y dimensions:

sample.hist(tof=200, y=200).plot(norm="log", vmin=1.0e-2)
../_images/8b0ab67091f2c10f975d169ec2e8288d9a78f1ad17e417589a0fdfe0f3997eb8.svg

Histogramming along y only gives a 1D plot:

sample.hist(y=200).plot(norm="log")
../_images/3a2385bca27acb119d6075aeacdbbd4c9d4410a7e50b1ac68fd0fa5b33bd927f.svg

Coordinate transformations#

The first step in the data reduction is to convert the raw event coordinates (position, time-of-flight) to something physically meaningful such as wavelength (\(\lambda\)) or momentum transfer (\(Q\)).

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")}
sc.show_graph(graph, simplified=True)
../_images/b2ba5ca46c8f506c34751deff73cd9cba2c2f875b2e18fb014d5ba4d773aef4c.svg

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

sample_wav = sample.transform_coords("wavelength", graph=graph)
sample_wav
Show/Hide data repr Show/Hide attributes
scipp.DataArray (786.24 MB)
    • event: 7927147
    • L1
      ()
      float64
      m
      150.0
      Values:
      array(150.)
    • L2
      (event)
      float64
      m
      3.002, 3.002, ..., 3.006, 3.006
      Values:
      array([3.00232538, 3.00237397, 3.00237952, ..., 3.00677258, 3.00624429,
       3.00617208])
    • Ltotal
      (event)
      float64
      m
      153.002, 153.002, ..., 153.006, 153.006
      Values:
      array([153.00232538, 153.00237397, 153.00237952, ..., 153.00677258,
       153.00624429, 153.00617208])
    • incident_beam
      ()
      vector3
      m
      [ 0. 0. 150.]
      Values:
      array([  0.,   0., 150.])
    • position
      (event)
      vector3
      m
      [0. 0.11814271 3. ], [0. 0.11937109 3. ], ..., [0. 0.19366142 3. ], [0. 0.19253724 3. ]
      Values:
      array([[0.        , 0.11814271, 3.        ],
       [0.        , 0.11937109, 3.        ],
       [0.        , 0.11951049, 3.        ],
       ...,
       [0.        , 0.20169624, 3.        ],
       [0.        , 0.19366142, 3.        ],
       [0.        , 0.19253724, 3.        ]])
    • sample_position
      ()
      vector3
      m
      [0. 0. 0.]
      Values:
      array([0., 0., 0.])
    • scattered_beam
      (event)
      vector3
      m
      [0. 0.11814271 3. ], [0. 0.11937109 3. ], ..., [0. 0.19366142 3. ], [0. 0.19253724 3. ]
      Values:
      array([[0.        , 0.11814271, 3.        ],
       [0.        , 0.11937109, 3.        ],
       [0.        , 0.11951049, 3.        ],
       ...,
       [0.        , 0.20169624, 3.        ],
       [0.        , 0.19366142, 3.        ],
       [0.        , 0.19253724, 3.        ]])
    • source_position
      ()
      vector3
      m
      [ 0. 0. -150.]
      Values:
      array([   0.,    0., -150.])
    • tof
      (event)
      float64
      ms
      253.831, 253.838, ..., 235.612, 235.615
      Values:
      array([253.83145495, 253.83804838, 253.83940871, ..., 257.95434162,
       235.61182539, 235.61501943])
    • wavelength
      (event)
      float64
      Å
      6.563, 6.563, ..., 6.092, 6.092
      Values:
      array([6.56307587, 6.56324426, 6.5632792 , ..., 6.66948353, 6.09183239,
       6.09191785])
    • x
      (event)
      float64
      m
      0.0, 0.0, ..., 0.0, 0.0
      Values:
      array([0., 0., 0., ..., 0., 0., 0.])
    • y
      (event)
      float64
      m
      0.118, 0.119, ..., 0.194, 0.193
      Values:
      array([0.11814271, 0.11937109, 0.11951049, ..., 0.20169624, 0.19366142,
       0.19253724])
    • (event)
      float64
      counts
      7.985e-09, 6.581e-13, ..., 4.828e-14, 4.341e-18
      Values:
      array([7.98463800e-09, 6.58130913e-13, 5.42462035e-17, ...,
       3.67754741e-28, 4.82836503e-14, 4.34070359e-18])

The result has a wavelength coordinate. We can also plot the result:

sample_wav.hist(wavelength=200).plot()
../_images/741bba0b5b4599d661b0d97cdac71a7185db2b697926bdf47f305a35ac912f2b.svg

We can see that the range of observed wavelengths agrees with the range set in the McStas model (5.25 - 6.75 Å)

Exercise 1: convert raw data to \(Q\)#

Instead of wavelength as in the example above, the task is now to convert the raw data to momentum-transfer \(Q\).

The transformation graph is missing the computation for \(Q\) so you will have to add it in yourself. As a reminder, \(Q\) is computed as follows

\[Q = \frac{4\pi \sin(\theta)}{\lambda}\]

You have to:

  • create a function that computes \(Q\)

  • add it to the graph

  • call transform_coords using the new graph

Solution:

Hide code cell content 
def compute_q(two_theta, wavelength):
    return (4.0 * np.pi) * sc.sin(two_theta / 2) / wavelength


graph["Q"] = compute_q

# Show the updated graph
display(sc.show_graph(graph, simplified=True))

# Run the coordinate transformation
sample_q = sample.transform_coords("Q", graph=graph)
sample_q
../_images/7c13132c96f1db9aeb205740074adcef49924b68b608f3fb36bbc2aa7850be8d.svg
Show/Hide data repr Show/Hide attributes
scipp.DataArray (907.19 MB)
    • event: 7927147
    • L1
      ()
      float64
      m
      150.0
      Values:
      array(150.)
    • L2
      (event)
      float64
      m
      3.002, 3.002, ..., 3.006, 3.006
      Values:
      array([3.00232538, 3.00237397, 3.00237952, ..., 3.00677258, 3.00624429,
       3.00617208])
    • Ltotal
      (event)
      float64
      m
      153.002, 153.002, ..., 153.006, 153.006
      Values:
      array([153.00232538, 153.00237397, 153.00237952, ..., 153.00677258,
       153.00624429, 153.00617208])
    • Q
      (event)
      float64
      1/Å
      0.038, 0.038, ..., 0.066, 0.066
      Values:
      array([0.03767955, 0.03806988, 0.03811408, ..., 0.06323093, 0.06647776,
       0.06609213])
    • incident_beam
      ()
      vector3
      m
      [ 0. 0. 150.]
      Values:
      array([  0.,   0., 150.])
    • position
      (event)
      vector3
      m
      [0. 0.11814271 3. ], [0. 0.11937109 3. ], ..., [0. 0.19366142 3. ], [0. 0.19253724 3. ]
      Values:
      array([[0.        , 0.11814271, 3.        ],
       [0.        , 0.11937109, 3.        ],
       [0.        , 0.11951049, 3.        ],
       ...,
       [0.        , 0.20169624, 3.        ],
       [0.        , 0.19366142, 3.        ],
       [0.        , 0.19253724, 3.        ]])
    • sample_position
      ()
      vector3
      m
      [0. 0. 0.]
      Values:
      array([0., 0., 0.])
    • scattered_beam
      (event)
      vector3
      m
      [0. 0.11814271 3. ], [0. 0.11937109 3. ], ..., [0. 0.19366142 3. ], [0. 0.19253724 3. ]
      Values:
      array([[0.        , 0.11814271, 3.        ],
       [0.        , 0.11937109, 3.        ],
       [0.        , 0.11951049, 3.        ],
       ...,
       [0.        , 0.20169624, 3.        ],
       [0.        , 0.19366142, 3.        ],
       [0.        , 0.19253724, 3.        ]])
    • source_position
      ()
      vector3
      m
      [ 0. 0. -150.]
      Values:
      array([   0.,    0., -150.])
    • tof
      (event)
      float64
      ms
      253.831, 253.838, ..., 235.612, 235.615
      Values:
      array([253.83145495, 253.83804838, 253.83940871, ..., 257.95434162,
       235.61182539, 235.61501943])
    • two_theta
      (event)
      float64
      rad
      0.039, 0.040, ..., 0.064, 0.064
      Values:
      array([0.03936057, 0.03976938, 0.03981578, ..., 0.06713105, 0.06446436,
       0.06409118])
    • wavelength
      (event)
      float64
      Å
      6.563, 6.563, ..., 6.092, 6.092
      Values:
      array([6.56307587, 6.56324426, 6.5632792 , ..., 6.66948353, 6.09183239,
       6.09191785])
    • x
      (event)
      float64
      m
      0.0, 0.0, ..., 0.0, 0.0
      Values:
      array([0., 0., 0., ..., 0., 0., 0.])
    • y
      (event)
      float64
      m
      0.118, 0.119, ..., 0.194, 0.193
      Values:
      array([0.11814271, 0.11937109, 0.11951049, ..., 0.20169624, 0.19366142,
       0.19253724])
    • (event)
      float64
      counts
      7.985e-09, 6.581e-13, ..., 4.828e-14, 4.341e-18
      Values:
      array([7.98463800e-09, 6.58130913e-13, 5.42462035e-17, ...,
       3.67754741e-28, 4.82836503e-14, 4.34070359e-18])

Histogram the data in \(Q\)#

The final step in processing the sample run is to histogram the data into \(Q\) bins.

sample_h = sample_q.hist(Q=200)
sample_h.plot(norm="log", vmin=1)
../_images/de7962026f94a306130a1b5ab07400bf59a3b5d1cddafe28a91f47a8b25058a3.svg

The histogrammed data currently has no standard deviations on the counts. This needs to be added after we have performed the histogramming operation.

When dealing with neutron events, we assume the data has a Poisson distribution. This means that the variance in a bin is equal to the counts in that bin (the standard deviation is then \(\sigma = \sqrt{\mathrm{counts}}\)).

We provide a helper function that will add Poisson variances to any given input:

utils.add_variances(sample_h)
sample_h.data
Show/Hide data repr Show/Hide attributes
scipp.Variable (3.38 KB)
    • (Q: 200)
      float64
      counts
      1203.199, 3838.421, ..., 0.107, 0.023
      σ = 34.687, 61.955, ..., 0.327, 0.152
      Values:
      array([1.20319921e+03, 3.83842052e+03, 4.16243475e+03, 3.80845060e+03,
       3.81889154e+03, 3.75953094e+03, 3.60715341e+03, 3.42202908e+03,
       3.41696459e+03, 3.04297520e+03, 2.98142937e+03, 2.77637369e+03,
       2.74380377e+03, 2.73597054e+03, 2.26035746e+03, 2.32372813e+03,
       2.17605644e+03, 2.15799149e+03, 1.73330252e+03, 1.67568229e+03,
       1.48492875e+03, 1.21852036e+03, 1.03449461e+03, 7.78452867e+02,
       4.89107567e+02, 4.16770967e+02, 2.81629075e+02, 3.07739170e+02,
       2.08926256e+02, 2.39772921e+02, 2.18546584e+02, 2.69402391e+02,
       2.42882128e+02, 2.26130124e+02, 2.20815427e+02, 2.08557827e+02,
       1.77100469e+02, 1.57605198e+02, 1.27170703e+02, 1.18295339e+02,
       8.44413374e+01, 7.24746483e+01, 5.22761562e+01, 4.04924693e+01,
       3.18159047e+01, 2.50157750e+01, 1.82234131e+01, 1.84665801e+01,
       1.79167641e+01, 1.57478223e+01, 2.19445884e+01, 2.31367855e+01,
       2.53466474e+01, 2.71640795e+01, 3.32310554e+01, 3.55565496e+01,
       3.68226899e+01, 3.90213554e+01, 4.20055217e+01, 4.65200876e+01,
       5.46723424e+01, 4.03857914e+01, 5.43462339e+01, 4.64065524e+01,
       4.74948881e+01, 4.94723753e+01, 4.62599716e+01, 4.61894402e+01,
       3.82319237e+01, 3.81788928e+01, 4.41858574e+01, 4.55374986e+01,
       2.63960116e+01, 4.04167809e+01, 2.82209614e+01, 2.57355011e+01,
       2.75838775e+01, 2.23910343e+01, 2.82767942e+01, 2.37990215e+01,
       1.87417033e+01, 1.81573345e+01, 1.91120313e+01, 1.67962592e+01,
       1.34354821e+01, 1.68201779e+01, 1.64625731e+01, 1.60895670e+01,
       1.42252180e+01, 1.20387849e+01, 1.64472661e+01, 1.72378536e+01,
       1.31960515e+01, 1.89912527e+01, 1.62530594e+01, 1.69286625e+01,
       2.24863204e+01, 1.81086714e+01, 2.11759979e+01, 2.47523526e+01,
       2.06744436e+01, 1.65472230e+01, 2.63489840e+01, 1.60743098e+01,
       2.20246085e+01, 1.99097352e+01, 2.25028989e+01, 1.81039938e+01,
       2.20971456e+01, 2.24417055e+01, 1.62910917e+01, 1.61656492e+01,
       1.80399523e+01, 1.94254774e+01, 1.72362368e+01, 1.40404499e+01,
       1.63209550e+01, 1.87918057e+01, 2.13742064e+01, 1.38328883e+01,
       2.02892640e+01, 1.68071718e+01, 1.79179080e+01, 1.57740789e+01,
       1.66757702e+01, 1.40609284e+01, 1.32936411e+01, 1.53559844e+01,
       2.11765968e+01, 2.33383054e+01, 1.33368682e+01, 2.01282136e+01,
       1.62903603e+01, 1.61440343e+01, 1.54520338e+01, 1.51022983e+01,
       1.69442429e+01, 1.65260148e+01, 1.65947838e+01, 1.62734246e+01,
       1.67368227e+01, 1.61661796e+01, 1.74394232e+01, 1.81259503e+01,
       1.50705837e+01, 1.56838866e+01, 1.73243290e+01, 1.64787133e+01,
       1.73517386e+01, 2.22572410e+01, 1.28226622e+01, 1.59654627e+01,
       2.14451934e+01, 1.34985819e+01, 1.34816771e+01, 2.08227519e+01,
       1.30541216e+01, 1.25606898e+01, 1.15144506e+01, 1.68134008e+01,
       1.59001770e+01, 2.03105312e+01, 1.63786669e+01, 1.37252214e+01,
       9.58447034e+00, 1.02095384e+01, 1.95574749e+01, 1.18401535e+01,
       1.19573046e+01, 1.13625853e+01, 9.44852584e+00, 1.51063660e+01,
       7.68975156e+00, 1.46107078e+01, 1.11963946e+01, 6.61584061e+00,
       9.07996074e+00, 8.91277717e+00, 6.29280673e+00, 9.19561981e+00,
       8.97817944e+00, 3.23734259e+00, 6.18888271e+00, 6.62608654e+00,
       4.84884112e+00, 6.27434973e+00, 4.89488596e+00, 5.53116022e+00,
       1.95896835e+00, 2.18534305e+00, 3.59895518e+00, 2.78997525e+00,
       4.47214725e+00, 2.65703239e+00, 6.25249261e-01, 3.26856043e+00,
       2.55131145e+00, 1.69336707e-01, 1.06709480e-01, 2.31543986e-02])

      Variances (σ²):
      array([1.20319921e+03, 3.83842052e+03, 4.16243475e+03, 3.80845060e+03,
       3.81889154e+03, 3.75953094e+03, 3.60715341e+03, 3.42202908e+03,
       3.41696459e+03, 3.04297520e+03, 2.98142937e+03, 2.77637369e+03,
       2.74380377e+03, 2.73597054e+03, 2.26035746e+03, 2.32372813e+03,
       2.17605644e+03, 2.15799149e+03, 1.73330252e+03, 1.67568229e+03,
       1.48492875e+03, 1.21852036e+03, 1.03449461e+03, 7.78452867e+02,
       4.89107567e+02, 4.16770967e+02, 2.81629075e+02, 3.07739170e+02,
       2.08926256e+02, 2.39772921e+02, 2.18546584e+02, 2.69402391e+02,
       2.42882128e+02, 2.26130124e+02, 2.20815427e+02, 2.08557827e+02,
       1.77100469e+02, 1.57605198e+02, 1.27170703e+02, 1.18295339e+02,
       8.44413374e+01, 7.24746483e+01, 5.22761562e+01, 4.04924693e+01,
       3.18159047e+01, 2.50157750e+01, 1.82234131e+01, 1.84665801e+01,
       1.79167641e+01, 1.57478223e+01, 2.19445884e+01, 2.31367855e+01,
       2.53466474e+01, 2.71640795e+01, 3.32310554e+01, 3.55565496e+01,
       3.68226899e+01, 3.90213554e+01, 4.20055217e+01, 4.65200876e+01,
       5.46723424e+01, 4.03857914e+01, 5.43462339e+01, 4.64065524e+01,
       4.74948881e+01, 4.94723753e+01, 4.62599716e+01, 4.61894402e+01,
       3.82319237e+01, 3.81788928e+01, 4.41858574e+01, 4.55374986e+01,
       2.63960116e+01, 4.04167809e+01, 2.82209614e+01, 2.57355011e+01,
       2.75838775e+01, 2.23910343e+01, 2.82767942e+01, 2.37990215e+01,
       1.87417033e+01, 1.81573345e+01, 1.91120313e+01, 1.67962592e+01,
       1.34354821e+01, 1.68201779e+01, 1.64625731e+01, 1.60895670e+01,
       1.42252180e+01, 1.20387849e+01, 1.64472661e+01, 1.72378536e+01,
       1.31960515e+01, 1.89912527e+01, 1.62530594e+01, 1.69286625e+01,
       2.24863204e+01, 1.81086714e+01, 2.11759979e+01, 2.47523526e+01,
       2.06744436e+01, 1.65472230e+01, 2.63489840e+01, 1.60743098e+01,
       2.20246085e+01, 1.99097352e+01, 2.25028989e+01, 1.81039938e+01,
       2.20971456e+01, 2.24417055e+01, 1.62910917e+01, 1.61656492e+01,
       1.80399523e+01, 1.94254774e+01, 1.72362368e+01, 1.40404499e+01,
       1.63209550e+01, 1.87918057e+01, 2.13742064e+01, 1.38328883e+01,
       2.02892640e+01, 1.68071718e+01, 1.79179080e+01, 1.57740789e+01,
       1.66757702e+01, 1.40609284e+01, 1.32936411e+01, 1.53559844e+01,
       2.11765968e+01, 2.33383054e+01, 1.33368682e+01, 2.01282136e+01,
       1.62903603e+01, 1.61440343e+01, 1.54520338e+01, 1.51022983e+01,
       1.69442429e+01, 1.65260148e+01, 1.65947838e+01, 1.62734246e+01,
       1.67368227e+01, 1.61661796e+01, 1.74394232e+01, 1.81259503e+01,
       1.50705837e+01, 1.56838866e+01, 1.73243290e+01, 1.64787133e+01,
       1.73517386e+01, 2.22572410e+01, 1.28226622e+01, 1.59654627e+01,
       2.14451934e+01, 1.34985819e+01, 1.34816771e+01, 2.08227519e+01,
       1.30541216e+01, 1.25606898e+01, 1.15144506e+01, 1.68134008e+01,
       1.59001770e+01, 2.03105312e+01, 1.63786669e+01, 1.37252214e+01,
       9.58447034e+00, 1.02095384e+01, 1.95574749e+01, 1.18401535e+01,
       1.19573046e+01, 1.13625853e+01, 9.44852584e+00, 1.51063660e+01,
       7.68975156e+00, 1.46107078e+01, 1.11963946e+01, 6.61584061e+00,
       9.07996074e+00, 8.91277717e+00, 6.29280673e+00, 9.19561981e+00,
       8.97817944e+00, 3.23734259e+00, 6.18888271e+00, 6.62608654e+00,
       4.84884112e+00, 6.27434973e+00, 4.89488596e+00, 5.53116022e+00,
       1.95896835e+00, 2.18534305e+00, 3.59895518e+00, 2.78997525e+00,
       4.47214725e+00, 2.65703239e+00, 6.25249261e-01, 3.26856043e+00,
       2.55131145e+00, 1.69336707e-01, 1.06709480e-01, 2.31543986e-02])
sample_h.plot(norm="log", vmin=1)
../_images/0aa26aae4b7585fc1fea8ccc513c621489095689470556bdabcdc07d497c2127.svg

Exercise 2: process flat-field run#

Repeat the step carried out above for the run that contained no sample (also known as a “flat-field” run).

Solution:

Hide code cell content 
folder = "../3-mcstas/SANS_without_sample_1_pulse"
# Load file
flat = utils.load_sans(folder)

# Convert to Q
flat_q = flat.transform_coords("Q", graph=graph)

# Histogram and add variances
flat_h = flat_q.hist(Q=200)
utils.add_variances(flat_h)
flat_h.plot()
../_images/174f6b5ead8ef613ba97db9ff1ec4c50c7f2d09e0785a6294194a67c41ce9891.svg

Bonus question: can you explain why the counts in the flat-field data drop at high \(Q\)?

Exercise 3: Normalize the sample run#

The flat-field run is giving a measure of the efficiency of each detector pixel. This efficiency now needs to be used to correct the counts in the sample run to yield a realistic \(I(Q)\) profile.

In particular, this should remove any unwanted artifacts in the data, such as the drop in counts around 0.03 Å-1 due to the air bubble inside the detector tube.

Normalizing is essentially just dividing the sample run by the flat field run.

Hint: you may run into an error like "Mismatch in coordinate 'Q'". Why is this happening? How can you get around it?

Solution:

Hide code cell content 
normed = sample_h / flat_h

---------------------------------------------------------------------------
DatasetError                              Traceback (most recent call last)
Cell In[14], line 1
----> 1 normed = sample_h / flat_h

DatasetError: Mismatch in coordinate 'Q' in operation 'divide':
(Q: 201)    float64           [1/Å]  [0.00451612, 0.0054472, ..., 0.189801, 0.190732]
vs
(Q: 201)    float64           [1/Å]  [0.00454553, 0.00547781, ..., 0.190068, 0.191001]
Hide code cell content 
# Make common bins
qbins = sc.linspace("Q", 5.0e-3, 0.2, 201, unit="1/angstrom")

# Histogram both sample and flat-field with same bins
sample_h = sample_q.hist(Q=qbins)
flat_h = flat_q.hist(Q=qbins)

# Add variances to both
utils.add_variances(sample_h, flat_h)

# Normalize
normed = sample_h / flat_h

normed.plot(norm="log", vmin=1.0e-3, vmax=10.0)
../_images/a0aaccd5e8a819a76dbff4d48cddeee814b90b3a10360b21a5e83eb3b87a0af8.svg

Save result to disk#

Finally, we need to save our results to disk, so that the reduced data can be forwarded to the next step in the pipeline (data analysis).

We will use a simple text file for this:

from scippneutron.io import save_xye

# The simple file format does not support bin-edge coordinates.
# So we convert to bin-centers first.
data = normed.copy()
data.coords["Q"] = sc.midpoints(data.coords["Q"])

save_xye("sans_iofq.dat", data)

Process data from 3 pulses#

We now want to repeat the reduction, but using more than a single pulse to improve our statistics.

We begin by loading the run with 3 pulses.

folder = "../3-mcstas/SANS_with_sample_3_pulse"
sample = utils.load_sans(folder)
sample
Show/Hide data repr Show/Hide attributes
scipp.DataArray (1.24 GB)
    • event: 23790295
    • position
      (event)
      vector3
      m
      [0. 0.35775631 3. ], [0. 0.19948785 3. ], ..., [0. 0.4845898 3. ], [0. 0.48404781 3. ]
      Values:
      array([[0.        , 0.35775631, 3.        ],
       [0.        , 0.19948785, 3.        ],
       [0.        , 0.20249984, 3.        ],
       ...,
       [0.        , 0.4809241 , 3.        ],
       [0.        , 0.4845898 , 3.        ],
       [0.        , 0.48404781, 3.        ]])
    • sample_position
      ()
      vector3
      m
      [0. 0. 0.]
      Values:
      array([0., 0., 0.])
    • source_position
      ()
      vector3
      m
      [ 0. 0. -150.]
      Values:
      array([   0.,    0., -150.])
    • tof
      (event)
      float64
      ms
      318.646, 331.640, ..., 333.437, 333.438
      Values:
      array([318.64622973, 331.63962524, 331.64618968, ..., 333.42472814,
       333.43705074, 333.43799127])
    • x
      (event)
      float64
      m
      0.0, 0.0, ..., 0.0, 0.0
      Values:
      array([0., 0., 0., ..., 0., 0., 0.])
    • y
      (event)
      float64
      m
      0.358, 0.199, ..., 0.485, 0.484
      Values:
      array([0.35775631, 0.19948785, 0.20249984, ..., 0.4809241 , 0.4845898 ,
       0.48404781])
    • (event)
      float64
      counts
      2.393e-08, 6.859e-14, ..., 1.500e-43, 1.205e-47
      Values:
      array([2.39281130e-08, 6.85931760e-14, 5.67098902e-18, ...,
       1.86665026e-39, 1.49958695e-43, 1.20470400e-47])
sample.hist(tof=200, y=200).plot(norm="log", vmin=1.0e-2)
../_images/91f88eaf00c12672dda9633c6581673941e9c59924c120828b2891d6b273a2fc.svg

We can see that we have 3 distinct pulses (3 horizontal bands of events).

In order to combine the data from the 3 pulses (also known as ‘pulse folding’), we need to compute a neutron time-of-flight relative to the start of each pulses, instead of absolute time.

Exercise 4: fold the pulses#

You are provided with a function that will perform the folding correctly, assuming you supply it the correct inputs.

help(utils.fold_pulses)

Help on function fold_pulses in module utils:

fold_pulses(data, tof_edges, offsets)
    Fold the data into a single pulse.
    
    Parameters
    ----------
    data
        Data to fold.
    tof_edges
        Edges of the time-of-flight bins.
    offsets
        Time offset to apply to each of the time-of-flight bins.

Instead of reading the bin edges off the figure above, you should try to compute your bin edges in a way that is based on some known parameters of your simulation (e.g. wavelength range) as well as the physical layout of the instrument.

Hints:

  • You have 3 pulses, so you need to supply 4 time-of-flight edges

  • Assume the first pulse has a zero time offset

  • The frequency of the ESS source is 14 Hz

  • To convert neutron wavelength to speed, you can use: \(~~~~v = \displaystyle\frac{h}{m_{\mathrm{n}}\lambda} \)

Solution:

Hide code cell content 
AA = sc.Unit("angstrom")
wmin = 5.25 * AA
wmax = 6.75 * AA
h = sc.constants.h
m_n = sc.constants.m_n

# Make 3 pulse offsets, starting at 3, using 14 pulses per second
pulse_offsets = sc.arange("tof", 3.0) * sc.scalar(1.0 / 14.0, unit="s").to(unit="ms")

# Compute minimum and maximum neutron speeds
vmin = (h / (m_n * wmax)).to(unit="m/s")
vmax = (h / (m_n * wmin)).to(unit="m/s")

# Compute pixel mean distance
mean_distance = sc.norm(
    sample.coords["position"].mean() - sample.coords["source_position"]
)

# Compute min and max times for neutrons to arrive
tmin = (mean_distance / vmax).to(unit="ms")
tmax = (mean_distance / vmin).to(unit="ms") + pulse_offsets[-1]

# Make evenly-spaced edges
edges = sc.linspace("tof", tmin.value, tmax.value, 4, unit="ms")

# Fold the data
folded = utils.fold_pulses(sample, edges, pulse_offsets)

# Convert to Q and histogram
folded_q = folded.transform_coords("Q", graph=graph)
folded_h = folded_q.hist(Q=qbins)
utils.add_variances(folded_h)

folded.hist(tof=200, y=200).plot(norm="log", vmin=1.0e-2)
../_images/c0b6d578a493af3f68504b36548d555604ff457772be77c70a269124cf640b02.svg 
folded_h.plot(norm="log", vmin=1.0)
../_images/4a6a5799673707fbfa5eca9d1d62d8c34e6a2cf36ccaa71199c5d4c2291b89df.svg
Hide code cell content 
folder = "../3-mcstas/SANS_without_sample_3_pulse"
flat = utils.load_sans(folder)

flat_folded = utils.fold_pulses(flat, edges, pulse_offsets)
flat_folded_q = flat_folded.transform_coords("Q", graph=graph)
flat_folded_h = flat_folded_q.hist(Q=qbins)
utils.add_variances(flat_folded_h)
Hide code cell content 
folded_normed = folded_h / flat_folded_h
folded_normed.plot(norm="log", vmin=1.0e-3, vmax=10.0)
../_images/8bc3d7d3e78eeafab8d80c51f8817a7285aba0aa32a92263b75bc82224f89ef8.svg 
pp.plot({"1-pulse": normed, "3-pulses": folded_normed}, norm="log", vmin=1.0e-3, vmax=10.0)
../_images/c286dbc45c522ad7491693c2a3d00de62a15eb99f8af90bcbdfbae577906ac2a.svg

Save results to disk#

Once again, we need to save our results to disk:

data = folded_normed.copy()
data.coords["Q"] = sc.midpoints(data.coords["Q"])

save_xye("sans_iofq_3pulses.dat", data)
Contents