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

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

Process the run with a sample

Load the NeXus file data

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

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

scipp.DataArray (436.93 MB)

• Dimensions:
  • event: 7158661
• Coordinates: (6)
  • position
    (event)
    vector3
    m
    [0. 0.02443003 3. ], [0. 0.01506518 3. ], ..., [0. 0.33851329 3. ], [0. 0.33860203 3. ]

    Values:
    array([[0.        , 0.02443003, 3.        ],
           [0.        , 0.01506518, 3.        ],
           [0.        , 0.01905563, 3.        ],
           ...,
           [0.        , 0.33891775, 3.        ],
           [0.        , 0.33851329, 3.        ],
           [0.        , 0.33860203, 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
    250.244, 242.829, ..., 223.609, 223.613

    Values:
    array([250.24446464, 242.82938911, 242.84190568, ..., 223.60865317,
           223.60938236, 223.61310145])

  • 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.024, 0.015, ..., 0.339, 0.339

    Values:
    array([0.02443003, 0.01506518, 0.01905563, ..., 0.33891775, 0.33851329,
           0.33860203])

• Data: (1)
  • (event)
    float64
    counts
    0.005, 5.749e-06, ..., 3.606e-31, 3.407e-35

    σ = 0.005, 5.749e-06, ..., 3.606e-31, 3.407e-35

    Values:
    array([4.79528262e-03, 5.74904359e-06, 4.95983440e-10, ...,
           3.81613715e-27, 3.60572915e-31, 3.40692228e-35])
    
    Variances (σ²):
    array([2.29947354e-05, 3.30515022e-11, 2.45999573e-19, ...,
           1.45629028e-53, 1.30012827e-61, 1.16071194e-69])

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)

Histogramming along y only gives a 1D plot:

sample.hist(y=200).plot(norm="log")

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)

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

scipp.DataArray (764.63 MB)

• Dimensions:
  • event: 7158661
• Coordinates: (12)
  • L1
    ()
    float64
    m
    150.0

    Values:
    array(150.)

  • L2
    (event)
    float64
    m
    3.000, 3.000, ..., 3.019, 3.019

    Values:
    array([3.00009947, 3.00003783, 3.00006052, ..., 3.01908351, 3.01903813,
           3.01904808])

  • Ltotal
    (event)
    float64
    m
    153.000, 153.000, ..., 153.019, 153.019

    Values:
    array([153.00009947, 153.00003783, 153.00006052, ..., 153.01908351,
           153.01903813, 153.01904808])

  • incident_beam
    ()
    vector3
    m
    [ 0. 0. 150.]

    Values:
    array([  0.,   0., 150.])

  • position
    (event)
    vector3
    m
    [0. 0.02443003 3. ], [0. 0.01506518 3. ], ..., [0. 0.33851329 3. ], [0. 0.33860203 3. ]

    Values:
    array([[0.        , 0.02443003, 3.        ],
           [0.        , 0.01506518, 3.        ],
           [0.        , 0.01905563, 3.        ],
           ...,
           [0.        , 0.33891775, 3.        ],
           [0.        , 0.33851329, 3.        ],
           [0.        , 0.33860203, 3.        ]])

  • sample_position
    ()
    vector3
    m
    [0. 0. 0.]

    Values:
    array([0., 0., 0.])

  • scattered_beam
    (event)
    vector3
    m
    [0. 0.02443003 3. ], [0. 0.01506518 3. ], ..., [0. 0.33851329 3. ], [0. 0.33860203 3. ]

    Values:
    array([[0.        , 0.02443003, 3.        ],
           [0.        , 0.01506518, 3.        ],
           [0.        , 0.01905563, 3.        ],
           ...,
           [0.        , 0.33891775, 3.        ],
           [0.        , 0.33851329, 3.        ],
           [0.        , 0.33860203, 3.        ]])

  • source_position
    ()
    vector3
    m
    [ 0. 0. -150.]

    Values:
    array([   0.,    0., -150.])

  • tof
    (event)
    float64
    ms
    250.244, 242.829, ..., 223.609, 223.613

    Values:
    array([250.24446464, 242.82938911, 242.84190568, ..., 223.60865317,
           223.60938236, 223.61310145])

  • wavelength
    (event)
    float64
    Å
    6.470, 6.279, ..., 5.781, 5.781

    Values:
    array([6.47042463, 6.27869989, 6.27902259, ..., 5.78100075, 5.78102131,
           5.78111709])

  • 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.024, 0.015, ..., 0.339, 0.339

    Values:
    array([0.02443003, 0.01506518, 0.01905563, ..., 0.33891775, 0.33851329,
           0.33860203])

• Data: (1)
  • (event)
    float64
    counts
    0.005, 5.749e-06, ..., 3.606e-31, 3.407e-35

    σ = 0.005, 5.749e-06, ..., 3.606e-31, 3.407e-35

    Values:
    array([4.79528262e-03, 5.74904359e-06, 4.95983440e-10, ...,
           3.81613715e-27, 3.60572915e-31, 3.40692228e-35])
    
    Variances (σ²):
    array([2.29947354e-05, 3.30515022e-11, 2.45999573e-19, ...,
           1.45629028e-53, 1.30012827e-61, 1.16071194e-69])

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

sample_wav.hist(wavelength=200).plot()

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

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

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

scipp.DataArray (873.86 MB)

• Dimensions:
  • event: 7158661
• Coordinates: (14)
  • L1
    ()
    float64
    m
    150.0

    Values:
    array(150.)

  • L2
    (event)
    float64
    m
    3.000, 3.000, ..., 3.019, 3.019

    Values:
    array([3.00009947, 3.00003783, 3.00006052, ..., 3.01908351, 3.01903813,
           3.01904808])

  • Ltotal
    (event)
    float64
    m
    153.000, 153.000, ..., 153.019, 153.019

    Values:
    array([153.00009947, 153.00003783, 153.00006052, ..., 153.01908351,
           153.01903813, 153.01904808])

  • Q
    (event)
    float64
    1/Å
    0.008, 0.005, ..., 0.122, 0.122

    Values:
    array([0.0079075 , 0.00502527, 0.00635599, ..., 0.12220343, 0.12205854,
           0.12208821])

  • incident_beam
    ()
    vector3
    m
    [ 0. 0. 150.]

    Values:
    array([  0.,   0., 150.])

  • position
    (event)
    vector3
    m
    [0. 0.02443003 3. ], [0. 0.01506518 3. ], ..., [0. 0.33851329 3. ], [0. 0.33860203 3. ]

    Values:
    array([[0.        , 0.02443003, 3.        ],
           [0.        , 0.01506518, 3.        ],
           [0.        , 0.01905563, 3.        ],
           ...,
           [0.        , 0.33891775, 3.        ],
           [0.        , 0.33851329, 3.        ],
           [0.        , 0.33860203, 3.        ]])

  • sample_position
    ()
    vector3
    m
    [0. 0. 0.]

    Values:
    array([0., 0., 0.])

  • scattered_beam
    (event)
    vector3
    m
    [0. 0.02443003 3. ], [0. 0.01506518 3. ], ..., [0. 0.33851329 3. ], [0. 0.33860203 3. ]

    Values:
    array([[0.        , 0.02443003, 3.        ],
           [0.        , 0.01506518, 3.        ],
           [0.        , 0.01905563, 3.        ],
           ...,
           [0.        , 0.33891775, 3.        ],
           [0.        , 0.33851329, 3.        ],
           [0.        , 0.33860203, 3.        ]])

  • source_position
    ()
    vector3
    m
    [ 0. 0. -150.]

    Values:
    array([   0.,    0., -150.])

  • tof
    (event)
    float64
    ms
    250.244, 242.829, ..., 223.609, 223.613

    Values:
    array([250.24446464, 242.82938911, 242.84190568, ..., 223.60865317,
           223.60938236, 223.61310145])

  • two_theta
    (event)
    float64
    rad
    0.008, 0.005, ..., 0.112, 0.112

    Values:
    array([0.00814316, 0.00502169, 0.00635179, ..., 0.11249561, 0.11236249,
           0.1123917 ])

  • wavelength
    (event)
    float64
    Å
    6.470, 6.279, ..., 5.781, 5.781

    Values:
    array([6.47042463, 6.27869989, 6.27902259, ..., 5.78100075, 5.78102131,
           5.78111709])

  • 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.024, 0.015, ..., 0.339, 0.339

    Values:
    array([0.02443003, 0.01506518, 0.01905563, ..., 0.33891775, 0.33851329,
           0.33860203])

• Data: (1)
  • (event)
    float64
    counts
    0.005, 5.749e-06, ..., 3.606e-31, 3.407e-35

    σ = 0.005, 5.749e-06, ..., 3.606e-31, 3.407e-35

    Values:
    array([4.79528262e-03, 5.74904359e-06, 4.95983440e-10, ...,
           3.81613715e-27, 3.60572915e-31, 3.40692228e-35])
    
    Variances (σ²):
    array([2.29947354e-05, 3.30515022e-11, 2.45999573e-19, ...,
           1.45629028e-53, 1.30012827e-61, 1.16071194e-69])

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)

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).

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

Solution:

Show code cell content

Hide code cell content

folder = "../3-mcstas/SANS_without_sample"
# Load file
flat = utils.load_sans(folder)

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

# Histogram
flat_h = flat_q.hist(Q=200)
flat_h.plot()

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:

Show code cell content

Hide code cell content

normed = sample_h / flat_h

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

DatasetError: Mismatch in coordinate 'Q' in operation 'divide':
(Q: 201)    float64           [1/Å]  [0.00452346, 0.00545678, ..., 0.190253, 0.191186]
vs
(Q: 201)    float64           [1/Å]  [0.00457172, 0.00550184, ..., 0.189665, 0.190595]

Show code cell content

Hide code cell content

# Make common bins, restricting the range to avoid NaNs
qbins = sc.linspace("Q", 5.0e-3, 0.19, 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)

# Normalize
normed = sample_h / flat_h

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

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)

Sciline workflow

Pause here if you have not yet been introduced to Sciline!

Below we build a Sciline workflow for reducing the SANS data. We wrap the operations that were carried out throughout this notebook into function that have the correct type-annotations, building a pipeline of connected steps.

The custom types have been created for you in the sans_utils module, and are imported from there.

Take some time to read through the code and understand the different parts. Ask questions if anything is unclear.

import sciline as sl
from sans_utils import *


def load(folder: Foldername[RunType]) -> RawData[RunType]:
    """Load raw data from file"""
    return RawData[RunType](utils.load_sans(folder))


def to_wavelength(
    data: RawData[RunType], graph: CoordTransformGraph
) -> WavelengthData[RunType]:
    """Compute wavelength for events"""
    return WavelengthData[RunType](data.transform_coords("wavelength", graph=graph))


def to_Q(data: WavelengthData[RunType], graph: CoordTransformGraph) -> QData[RunType]:
    """Compute Q for events"""
    return QData[RunType](data.transform_coords("Q", graph=graph))


def to_histogram(events: QData[RunType], qbins: QBins) -> QHistogram[RunType]:
    """Histogram data in Q bins"""
    return QHistogram[RunType](events.hist(Q=qbins))


def normalize(
    sample: QHistogram[SampleRun], flat: QHistogram[FlatFieldRun]
) -> NormalizedQ:
    """Normalized sample by flat-field"""
    return NormalizedQ(sample / flat)


wf = sl.Pipeline(
    (load, to_wavelength, to_Q, to_histogram, normalize),
    constraints={RunType: [SampleRun, FlatFieldRun]},
)

wf.visualize()

Exercise 4: Set the workflow parameters

In the visualization above, the red boxes indicate missing values; parameters required by the computation that have yet to be set on the workflow.

Using the wf[TYPE] = VALUE syntax, set the missing parameters on the workflow.

Once you have done so, use wf.compute(NormalizedQ) to compute the final normalized result.

Solution:

Show code cell content

Hide code cell content

wf[Foldername[SampleRun]] = "../3-mcstas/SANS_with_sample"
wf[Foldername[FlatFieldRun]] = "../3-mcstas/SANS_without_sample"

wf[CoordTransformGraph] = graph
wf[QBins] = qbins

iofq = wf.compute(NormalizedQ)
iofq.plot(norm="log")

Exercise 5: Make a \((2\theta, Q)\) map

We wish to make a 2d plot of counts as a function of \(2\theta\) and \(Q\) for the SampleRun. You have to:

• identify one of the intermediate results that would contain the necessary information

• compute that intermediate result

• histogram the result (200 bins in two_theta and 200 bins in Q will do nicely)

• plot it on a figure

Hint: Using norm='log', vmin=1.0e-4 when plotting makes for a better-looking figure.

Solution:

Show code cell content

Hide code cell content

da = wf.compute(QData[SampleRun])
da.hist(two_theta=200, Q=200).plot(norm="log", vmin=1.0e-4)