Bragg-edge imaging with ODIN#

This notebook illustrates how to convert recorded events on the ODIN detector to a single wavelength spectrum, revealing a Bragg edge in the data. WFM mode was used in the chopper cascade.

Loading dataset#

Loader is not part of essimaging since McStas dataset format is not stabilized yet.

[1]:

import scipp as sc
import scippnexus as snx
import scipp.constants as scc
from typing import cast, NewType
from ess.reduce.nexus.types import FilePath


_DataPath = NewType('_DataPath', str)
_DefaultDataPath = _DataPath(
    "entry1/data/transmission_event_signal_dat_list_p_t_x_y_z_vx_vy_vz/events"
)
_FileLock = NewType('_FileLock', bool)
"""Lock the file to prevent concurrent access."""
_DefaultFileLock = _FileLock(True)
OdinSimulationRawData = NewType('OdinSimulationRawData', sc.DataArray)
ProbabilityToCountsScaleFactor = NewType('ProbabilityToCountsScaleFactor', sc.Variable)
"""Translate the probability to counts."""
DefaultProbabilityToCountsScaleFactor = ProbabilityToCountsScaleFactor(
    sc.scalar(1_000, unit='dimensionless')
)
DetectorStartX = NewType('DetectorStartX', sc.Variable)
"""Start of the detector in x direction."""
DefaultDetectorStartX = DetectorStartX(sc.scalar(-0.03, unit='m'))
DetectorStartY = NewType('DetectorStartY', sc.Variable)
"""Start of the detector in y direction."""
DefaultDetectorStartY = DetectorStartY(sc.scalar(-0.03, unit='m'))

DetectorEndX = NewType('DetectorEndX', sc.Variable)
"""End of the detector in x direction."""
DefaultDetectorEndX = DetectorEndX(sc.scalar(0.03, unit='m'))
DetectorEndY = NewType('DetectorEndY', sc.Variable)
"""End of the detector in y direction."""
DefaultDetectorEndY = DetectorEndY(sc.scalar(0.03, unit='m'))

McStasManualResolution = NewType('McStasManualResolution', tuple)
"""Manual resolution for McStas data (how many pixels per axis x, y)"""
DefaultMcStasManualResolution = McStasManualResolution((1024, 1024))

example_resolution = McStasManualResolution((128, 128))
# Small resolution for faster testing and documentation build.


def _nth_col_or_row_lookup(
    start: sc.Variable, stop: sc.Variable, resolution: int, dim: str
) -> sc.Lookup:
    """Lookup the nth column or row."""
    position = sc.linspace(
        dim, start=start, stop=stop, num=resolution + 1, unit=start.unit
    )
    nth_col_or_row = sc.arange(dim=dim, start=0, stop=resolution, unit='dimensionless')
    hist = sc.DataArray(data=nth_col_or_row, coords={dim: position})
    return sc.lookup(hist, dim)


def _position_to_pixel_id(
    *,
    x_pos: sc.Variable,
    y_pos: sc.Variable,
    detector_start_x: DetectorStartX = DefaultDetectorStartX,
    detector_start_y: DetectorStartY = DefaultDetectorStartY,
    detector_end_x: DetectorEndX = DefaultDetectorEndX,
    detector_end_y: DetectorEndY = DefaultDetectorEndY,
    resolution: McStasManualResolution = DefaultMcStasManualResolution,
) -> sc.Variable:
    """Hardcode pixel ids from positions."""
    x_position_lookup = _nth_col_or_row_lookup(
        detector_start_x, detector_end_x, resolution[0], 'x'
    )
    y_position_lookup = _nth_col_or_row_lookup(
        detector_start_y, detector_end_y, resolution[1], 'y'
    )
    n_cols = x_position_lookup[x_pos]
    n_rows = y_position_lookup[y_pos]
    return n_rows * resolution[0] + n_cols


McStasVelocities = NewType('McStasVelocities', sc.DataGroup)


def load_velocities(
    file_path: FilePath,
    _data_path: _DataPath = _DefaultDataPath,
    _file_lock: _FileLock = _DefaultFileLock,
) -> McStasVelocities:
    with snx.File(file_path, "r", locking=_file_lock) as f:
        data = f[_data_path][()].rename_dims({'dim_0': 'event'})
        velocities = data['dim_1', 5:8]
        vx = cast(sc.Variable, velocities['dim_1', 0].copy())
        vy = cast(sc.Variable, velocities['dim_1', 1].copy())
        vz = cast(sc.Variable, velocities['dim_1', 2].copy())
        for v_component in (vx, vy, vz):
            v_component.unit = 'm/s'
        # Add special tags if you want to use them as coordinates
        # for example, da.coords['vx_MC'] = vx
        # to distinguish them from the measurement
        return McStasVelocities(sc.DataGroup(vx=vx, vy=vy, vz=vz))


LoadTrueVelocities = NewType('LoadTrueVelocities', bool)
DefaultLoadTrueVelocities = LoadTrueVelocities(True)


def load_odin_simulation_data(
    file_path: FilePath,
    _data_path: _DataPath = _DefaultDataPath,
    _file_lock: _FileLock = _DefaultFileLock,
    detector_start_x: DetectorStartX = DefaultDetectorStartX,
    detector_start_y: DetectorStartY = DefaultDetectorStartY,
    detector_end_x: DetectorEndX = DefaultDetectorEndX,
    detector_end_y: DetectorEndY = DefaultDetectorEndY,
    resolution: McStasManualResolution = DefaultMcStasManualResolution,
    probability_scale_factor: ProbabilityToCountsScaleFactor = DefaultProbabilityToCountsScaleFactor,
    load_true_velocities: LoadTrueVelocities = DefaultLoadTrueVelocities,
) -> OdinSimulationRawData:
    with snx.File(file_path, "r", locking=_file_lock) as f:
        # The name p_t_x_y_z_vx_vy_vz represents
        # probability, time of arrival, position(x, y, z) and velocity(vx, vy, vz).
        # The name also represents the order of each field in the table.
        # For example, probability is the first field, so data['dim_1', 0] is the probability.
        data = f[_data_path][()].rename_dims({'dim_0': 'event'})
        probabilities = cast(sc.Variable, data['dim_1', 0].copy())
        probabilities.unit = 'dimensionless'
        time_of_arrival = cast(sc.Variable, data['dim_1', 1].copy())
        time_of_arrival.unit = 's'  # Hardcoded unit from the data.
        positions = data['dim_1', 2:5]
        counts = (probabilities / probabilities.max()) * probability_scale_factor
        counts.unit = 'counts'
        # Units are hardcoded from the data.
        x_pos = cast(sc.Variable, positions['dim_1', 0].copy())
        x_pos.unit = 'm'
        y_pos = cast(sc.Variable, positions['dim_1', 1].copy())
        y_pos.unit = 'm'
        pixel_id = _position_to_pixel_id(
            x_pos=x_pos,
            y_pos=y_pos,
            detector_start_x=detector_start_x,
            detector_start_y=detector_start_y,
            detector_end_x=detector_end_x,
            detector_end_y=detector_end_y,
            resolution=resolution,
        )
        da = sc.DataArray(
            data=counts.copy().astype(sc.DType.int32),
            coords={
                'time_of_arrival': time_of_arrival.to(unit='us'),
                'sample_position': sc.vector([0.0, 0.0, 60.5], unit='m'),
                # Hardcoded from the data.
                'source_position': sc.vector([0.0, 0.0, 0.0], unit="m"),
                # Hardcoded from the data.
                'pixel_id': pixel_id,
            },
        )
        if load_true_velocities:
            velocities = load_velocities(file_path, _data_path, _file_lock)
            speeds = sc.norm(
                sc.vectors(
                    dims=['event'],
                    values=sc.transpose(
                        sc.concat(list(velocities.values()), 'speed')
                    ).values,
                    unit='m/s',
                )
            )
            da.coords['sim_wavelength'] = (scc.h / scc.neutron_mass / speeds).to(
                unit='angstrom'
            )

        return OdinSimulationRawData(da.to(dtype=float))

[2]:

from ess.imaging.data import get_mcstas_ob_images_path, get_mcstas_sample_images_path

ob_file_path = FilePath(get_mcstas_ob_images_path())
sample_file_path = FilePath(get_mcstas_sample_images_path())
ob_da = load_odin_simulation_data(ob_file_path, resolution=example_resolution)
sample_da = load_odin_simulation_data(sample_file_path, resolution=example_resolution)
sample_da

Downloading file 'small_mcstas_ob_images.h5' from 'https://public.esss.dk/groups/scipp/ess/imaging/1/small_mcstas_ob_images.h5' to '/home/runner/.cache/essimaging/1'.
Downloading file 'small_mcstas_sample_images.h5' from 'https://public.esss.dk/groups/scipp/ess/imaging/1/small_mcstas_sample_images.h5' to '/home/runner/.cache/essimaging/1'.

[2]:
Show/Hide data repr Show/Hide attributes
scipp.DataArray (11.86 MB)
    • event: 388566
    • pixel_id
      (event)
      int64
      𝟙
      8122, 5408, ..., 15774, 969
      Values:
      array([ 8122,  5408,  4970, ...,  8714, 15774,   969], shape=(388566,))
    • sample_position
      ()
      vector3
      m
      [ 0. 0. 60.5]
      Values:
      array([ 0. ,  0. , 60.5])
    • sim_wavelength
      (event)
      float64
      Å
      5.727, 5.673, ..., 8.380, 6.448
      Values:
      array([5.72678489, 5.67272213, 3.1759751 , ..., 7.56736601, 8.38023051,
       6.44810453], shape=(388566,))
    • source_position
      ()
      vector3
      m
      [0. 0. 0.]
      Values:
      array([0., 0., 0.])
    • time_of_arrival
      (event)
      float64
      µs
      8.994e+04, 8.894e+04, ..., 1.304e+05, 9.970e+04
      Values:
      array([ 89942.25084908,  88942.61380286,  50284.40890975, ...,
       117180.5024966 , 130417.13999026,  99700.22439706], shape=(388566,))
    • (event)
      float64
      counts
      22.0, 0.0, ..., 3.0, 0.0
      Values:
      array([ 22.,   0., 198., ...,   0.,   3.,   0.], shape=(388566,))
[3]:

def _pixel_ids_to_x(
    *,
    pixel_id: sc.Variable,
    resolution: McStasManualResolution = DefaultMcStasManualResolution,
    detector_start_x: DetectorStartX = DefaultDetectorStartX,
    detector_end_x: DetectorEndX = DefaultDetectorEndX,
) -> sc.Variable:
    n_col = pixel_id % resolution[0]
    x_interval = (detector_end_x - detector_start_x) / resolution[0]
    return (
        detector_start_x + n_col * x_interval
    ) + x_interval / 2  # Center of the pixel|


def _pixel_ids_to_y(
    *,
    pixel_id: sc.Variable,
    resolution: McStasManualResolution = DefaultMcStasManualResolution,
    detector_start_y: DetectorStartY = DefaultDetectorStartY,
    detector_end_y: DetectorEndY = DefaultDetectorEndY,
) -> sc.Variable:
    n_row = pixel_id // resolution[0]
    y_interval = (detector_end_y - detector_start_y) / resolution[1]
    return (
        detector_start_y + n_row * y_interval
    ) + y_interval / 2  # Center of the pixel


def _pixel_ids_to_position(
    *, x: sc.Variable, y: sc.Variable, z_pos: sc.Variable
) -> sc.Variable:
    z = sc.zeros_like(x) + z_pos
    var = (
        sc.concat([x, y, z], 'event')
        .fold('event', dims=['pos', 'event'], shape=[3, len(x)])
        .transpose(dims=['event', 'pos'])
        .values
    )
    return sc.vectors(dims=['event'], values=var, unit='m')

[4]:

import scipp as sc
from scippneutron.conversion import graph

plane_graph = {**graph.beamline.beamline(False), **graph.tof.kinematic("tof")}

# TODO: Replace this with actual WFM stitching method
plane_graph['tof'] = lambda time_of_arrival: time_of_arrival
plane_graph['x'] = lambda pixel_id: _pixel_ids_to_x(
    pixel_id=pixel_id, resolution=example_resolution
)
plane_graph['y'] = lambda pixel_id: _pixel_ids_to_y(
    pixel_id=pixel_id, resolution=example_resolution
)
plane_graph['position'] = lambda x, y: _pixel_ids_to_position(
    x=x,
    y=y,
    z_pos=sc.scalar(60.5, unit='m'),  # Hardcoded from the data.
)

sc.show_graph(plane_graph, simplified=True)

[4]:
../_images/user-guide_odin_simulation_5_0.svg 
[5]:

coords = ["tof", "position", "x", "y", "sim_wavelength", "Ltotal"]

sample_da = sample_da.transform_coords(
    coords,
    graph=plane_graph,
    keep_intermediate=False,
)
ob_da = ob_da.transform_coords(
    coords,
    graph=plane_graph,
    keep_intermediate=False,
)

sample_da

[5]:
Show/Hide data repr Show/Hide attributes
scipp.DataArray (32.61 MB)
    • event: 388566
    • Ltotal
      (event)
      float64
      m
      60.500, 60.500, ..., 60.500, 60.500
      Values:
      array([60.50000006, 60.50000264, 60.50000446, ..., 60.50000523,
       60.50000847, 60.50000596], shape=(388566,))
    • pixel_id
      (event)
      int64
      𝟙
      8122, 5408, ..., 15774, 969
      Values:
      array([ 8122,  5408,  4970, ...,  8714, 15774,   969], shape=(388566,))
    • position
      (event)
      vector3
      m
      [-2.578125e-03 -2.343750e-04 6.050000e+01], [-1.4765625e-02 -1.0078125e-02 6.0500000e+01], ..., [-1.5703125e-02 2.7890625e-02 6.0500000e+01], [ 4.4531250e-03 -2.6484375e-02 6.0500000e+01]
      Values:
      array([[-2.5781250e-03, -2.3437500e-04,  6.0500000e+01],
       [-1.4765625e-02, -1.0078125e-02,  6.0500000e+01],
       [ 1.9921875e-02, -1.1953125e-02,  6.0500000e+01],
       ...,
       [-2.5078125e-02,  2.1093750e-03,  6.0500000e+01],
       [-1.5703125e-02,  2.7890625e-02,  6.0500000e+01],
       [ 4.4531250e-03, -2.6484375e-02,  6.0500000e+01]],
      shape=(388566, 3))
    • sample_position
      ()
      vector3
      m
      [ 0. 0. 60.5]
      Values:
      array([ 0. ,  0. , 60.5])
    • sim_wavelength
      (event)
      float64
      Å
      5.727, 5.673, ..., 8.380, 6.448
      Values:
      array([5.72678489, 5.67272213, 3.1759751 , ..., 7.56736601, 8.38023051,
       6.44810453], shape=(388566,))
    • source_position
      ()
      vector3
      m
      [0. 0. 0.]
      Values:
      array([0., 0., 0.])
    • time_of_arrival
      (event)
      float64
      µs
      8.994e+04, 8.894e+04, ..., 1.304e+05, 9.970e+04
      Values:
      array([ 89942.25084908,  88942.61380286,  50284.40890975, ...,
       117180.5024966 , 130417.13999026,  99700.22439706], shape=(388566,))
    • tof
      (event)
      float64
      µs
      8.994e+04, 8.894e+04, ..., 1.304e+05, 9.970e+04
      Values:
      array([ 89942.25084908,  88942.61380286,  50284.40890975, ...,
       117180.5024966 , 130417.13999026,  99700.22439706], shape=(388566,))
    • x
      (event)
      float64
      m
      -0.003, -0.015, ..., -0.016, 0.004
      Values:
      array([-0.00257812, -0.01476562,  0.01992187, ..., -0.02507812,
       -0.01570313,  0.00445313], shape=(388566,))
    • y
      (event)
      float64
      m
      -0.000, -0.010, ..., 0.028, -0.026
      Values:
      array([-0.00023438, -0.01007812, -0.01195313, ...,  0.00210938,
        0.02789062, -0.02648438], shape=(388566,))
    • (event)
      float64
      counts
      22.0, 0.0, ..., 3.0, 0.0
      Values:
      array([ 22.,   0., 198., ...,   0.,   3.,   0.], shape=(388566,))
[6]:

sample_da.hist(tof=300).plot()

[6]:
../_images/user-guide_odin_simulation_7_0.svg

Convert McStas raw data to NeXus#

The raw McStas data looks different from what data in a NeXus file would look like. The time-of-flight recorded by the McStas monitor is a unwrapped time of arrival (see https://scipp.github.io/scippneutron/user-guide/chopper/frame-unwrapping.html); the tof coordinate has values beyond 71ms, as can be seen in the plot above.

The workflow that computes wavelengths from the WFM chopper cascade expects data in the NeXus format, so we transform the data here.

[7]:

def to_nexus(da):
    unit = da.coords['tof'].unit
    period = (1.0 / sc.scalar(14., unit='Hz')).to(unit=unit)
    # Bin the data into bins with a 71ms period
    n = int(sample_da.coords['tof'].max() / period)
    da = da.bin(tof=sc.arange('tof', n + 2) * period)
    # Add a event_time_zero coord for each bin, but not as bin edges, as all events in the same pulse have the same event_time_zero, hence the `[:2]`
    da.coords['event_time_zero'] = (sc.scalar(1730450434078980000, unit='ns').to(unit=unit) + da.coords['tof'])[:-1]
    # Remove the meaningless tof coord at the top level
    del da.coords['tof']
    # Remove the original (wrong) event_time_zero event coord inside the bins and rename the dim
    del da.bins.coords['time_of_arrival']
    del da.bins.coords['Ltotal']
    da = da.rename_dims(tof='event_time_zero')
    # Compute a proper event_time_offset as tof % period
    da.bins.coords['event_time_offset'] = (da.bins.coords.pop('tof') % period)#.to(unit=)
    return da

def add_positions(da):
    temp = da.bins.concat('event_time_zero').copy()
    out = da.copy()
    out.coords['position'] =  temp.bins.coords['position'].bins.mean()
    del out.bins.coords['position']
    return out.transform_coords(
        "Ltotal",
        graph=plane_graph,
        keep_intermediate=True,
    )

sample_nexus = add_positions(to_nexus(sample_da).group('pixel_id'))
ob_nexus = add_positions(to_nexus(ob_da).group('pixel_id'))

sample_nexus

[7]:
Show/Hide data repr Show/Hide attributes
scipp.DataArray (16.20 MB)
    • event_time_zero: 3
    • pixel_id: 16384
    • Ltotal
      (pixel_id)
      float64
      m
      60.500, 60.500, ..., 60.500, 60.500
      Values:
      array([60.50001464, 60.50001442, 60.50001419, ..., 60.50001419,
       60.50001442, 60.50001464], shape=(16384,))
    • event_time_zero
      (event_time_zero)
      float64
      µs
      1.730e+15, 1.730e+15, 1.730e+15
      Values:
      array([1.73045043e+15, 1.73045043e+15, 1.73045043e+15])
    • pixel_id
      (pixel_id)
      int64
      𝟙
      0, 1, ..., 16382, 16383
      Values:
      array([    0,     1,     2, ..., 16381, 16382, 16383], shape=(16384,))
    • position
      (pixel_id)
      vector3
      m
      [-2.9765625e-02 -2.9765625e-02 6.0500000e+01], [-2.9296875e-02 -2.9765625e-02 6.0500000e+01], ..., [2.9296875e-02 2.9765625e-02 6.0500000e+01], [2.9765625e-02 2.9765625e-02 6.0500000e+01]
      Values:
      array([[-2.9765625e-02, -2.9765625e-02,  6.0500000e+01],
       [-2.9296875e-02, -2.9765625e-02,  6.0500000e+01],
       [-2.8828125e-02, -2.9765625e-02,  6.0500000e+01],
       ...,
       [ 2.8828125e-02,  2.9765625e-02,  6.0500000e+01],
       [ 2.9296875e-02,  2.9765625e-02,  6.0500000e+01],
       [ 2.9765625e-02,  2.9765625e-02,  6.0500000e+01]], shape=(16384, 3))
    • sample_position
      ()
      vector3
      m
      [ 0. 0. 60.5]
      Values:
      array([ 0. ,  0. , 60.5])
    • source_position
      ()
      vector3
      m
      [0. 0. 0.]
      Values:
      array([0., 0., 0.])
    • (event_time_zero, pixel_id)
      DataArrayView
      binned data [len=5, len=2, ..., len=0, len=2]
      dim='event',
content=DataArray(
          dims=(event: 388566),
          data=float64[counts],
          coords={'sim_wavelength':float64[Å], 'y':float64[m], 'x':float64[m],
                  'event_time_offset':float64[µs]})
[8]:

# Visualize
fig_nexus = sample_nexus.bins.concat().hist(event_time_offset=300).plot(title='McStas simulation: sample')
fig_nexus + ob_nexus.bins.concat().hist(event_time_offset=300).plot(title='McStas simulation: open beam')

[8]:
../_images/user-guide_odin_simulation_10_0.svg 
[9]:

import plopp as pp

pp.scatter3d(sample_nexus.sum('event_time_zero'), pos='position', cbar=True, pixel_size=0.0005)

[9]:

Choppers#

To accurately compute the wavelengths of the neutrons from their time-of-arrival, we need the parameters of the choppers in the beamline.

[10]:

import sciline as sl
from scippneutron.chopper import DiskChopper
from scippneutron.tof import unwrap
from scippneutron.tof import chopper_cascade

Hz = sc.Unit("Hz")
deg = sc.Unit("deg")
meter = sc.Unit("m")

parameters = {
    "WFMC_1": {
        "frequency": 56.0,
        "phase": 93.244,
        "distance": 6.85,
        "open": [-1.9419, 49.5756, 98.9315, 146.2165, 191.5176, 234.9179],
        "close": [1.9419, 55.7157, 107.2332, 156.5891, 203.8741, 249.1752]
    },
    "WFMC_2": {
        "frequency": 56.0,
        "phase": 97.128,
        "distance": 7.15,
        "open": [-1.9419, 51.8318, 103.3493, 152.7052, 199.9903, 245.2914],
        "close": [1.9419, 57.9719, 111.6510, 163.0778, 212.3468, 259.5486]
    },
    "FOC_1": {
        "frequency": 42.0,
        "phase": 81.303297,
        "distance": 8.4,
        "open": [-5.1362, 42.5536, 88.2425, 132.0144, 173.9497, 216.7867],
        "close": [5.1362, 54.2095, 101.2237, 146.2653, 189.417, 230.7582]
    },
    "BP_1": {
        "frequency": 7.0,
        "phase": 31.080,
        "distance": 8.45,
        "open": [-23.6029],
        "close": [23.6029]
    },
    "FOC_2": {
        "frequency": 42.0,
        "phase": 107.013442,
        "distance": 12.2,
        "open": [-16.3227, 53.7401, 120.8633, 185.1701, 246.7787, 307.0165],
        "close": [16.3227, 86.8303, 154.3794, 218.7551, 280.7508, 340.3188]
    },
    "BP_2": {
        "frequency": 7.0,
        "phase": 44.224,
        "distance": 12.25,
        "open": [-34.4663],
        "close": [34.4663]
    },
    "T0_alpha": {
        "frequency": 14.0,
        "phase": 179.672,
        "distance": 13.5,
        "open": [-167.8986],
        "close": [167.8986]
    },
    "T0_beta": {
        "frequency": 14.0,
        "phase": 179.672,
        "distance": 13.7,
        "open": [-167.8986],
        "close": [167.8986]
    },
    "FOC_3": {
        "frequency": 28.0,
        "phase": 92.993,
        "distance": 17.0,
        "open": [-20.302, 45.247, 108.0457, 168.2095, 225.8489, 282.2199],
        "close": [20.302, 85.357, 147.6824, 207.3927, 264.5977, 319.4024]
    },
    "FOC_4": {
        "frequency": 14.0,
        "phase": 61.584,
        "distance": 23.69,
        "open": [-16.7157, 29.1882, 73.1661, 115.2988, 155.6636, 195.5254],
        "close": [16.7157, 61.8217, 105.0352, 146.4355, 186.0987, 224.0978]
    },
    "FOC_5": {
        "frequency": 14.0,
        "phase": 82.581,
        "distance": 33.0,
        "open": [-25.8514, 38.3239, 99.8064, 160.1254, 217.4321, 272.5426],
        "close": [25.8514, 88.4621, 147.4729, 204.0245, 257.7603, 313.7139]
    },

}

disk_choppers = {key: DiskChopper(
    frequency=-ch["frequency"] * Hz,
    beam_position=sc.scalar(0.0, unit="deg"),
    phase=-ch["phase"] * deg,
    axle_position=sc.vector(value=[0, 0, ch["distance"]], unit="m"),
    slit_begin=sc.array(dims=["cutout"], values=ch["open"], unit="deg"),
    slit_end=sc.array(dims=["cutout"], values=ch["close"], unit="deg")
) for key, ch in parameters.items() }

[11]:

disk_choppers["WFMC_1"]

[11]:
  • axle_position
    scipp
    Variable
    ()
    vector3
    m
    [0. 0. 6.85]
  • frequency
    scipp
    Variable
    ()
    float64
    Hz
    -56.0
  • beam_position
    scipp
    Variable
    ()
    float64
    deg
    0.0
  • phase
    scipp
    Variable
    ()
    float64
    deg
    -93.244
  • slit_begin
    scipp
    Variable
    (cutout: 6)
    float64
    deg
    -1.942, 49.576, ..., 191.518, 234.918
  • slit_end
    scipp
    Variable
    (cutout: 6)
    float64
    deg
    1.942, 55.716, ..., 203.874, 249.175
  • slit_height
    NoneType
    ()
    None
  • radius
    NoneType
    ()
    None
begin0 end0 begin1 end1 begin2 end2 begin3 end3 begin4 end4 begin5 end5 TDC beam position 
[12]:

choppers = {
    key: chopper_cascade.Chopper.from_disk_chopper(
        chop,
        pulse_frequency=sc.scalar(14.0, unit="Hz"),
        npulses=1,
    )
    for key, chop in disk_choppers.items()
}

Check that the chopper settings make sense with a quick tof run#

As useful sanity check is to run a basic simulation, propagating neutrons through the chopper cascade, using the Tof package.

[13]:

from scippneutron.tof.fakes import FakeBeamlineEss

Ltotal = sample_da.coords['Ltotal'].mean()
ess_beamline = FakeBeamlineEss(
    choppers=choppers,
    monitors={"detector": Ltotal},
    run_length=sc.scalar(1 / 14, unit="s") * 8,
    events_per_pulse=100_000,
)

ess_beamline.model_result.plot()

[13]:

Plot(ax=<Axes: xlabel='Time-of-flight (us)', ylabel='Distance (m)'>, fig=<Figure size 1200x480 with 2 Axes>)
../_images/user-guide_odin_simulation_17_1.png

We observe that the WFM choppers make 6 distinct frames at the detector, and that the other choppers skip every other pulse to maximize wavelength coverage.

We can now compare the counts on the detector to our raw data, to make sure they broadly resemble each other.

[14]:

raw_data = ess_beamline.get_monitor("detector")[0]

# Visualize
fig_nexus + raw_data.hist(event_time_offset=300).sum("pulse").plot(title='Tof simulation')

[14]:
../_images/user-guide_odin_simulation_19_0.svg

Use WFM workflow#

We now set up the workflow which will convert the raw neutron arrival times to a real time-of-flight, and thus a wavelength.

Chopper cascade#

[15]:

one_pulse = ess_beamline.source.data["pulse", 0]
pulse_tmin = one_pulse.coords["time"].min()
pulse_tmax = one_pulse.coords["time"].max()
pulse_wmin = one_pulse.coords["wavelength"].min()
pulse_wmax = one_pulse.coords["wavelength"].max()

frames = chopper_cascade.FrameSequence.from_source_pulse(
    time_min=pulse_tmin,
    time_max=pulse_tmax,
    wavelength_min=pulse_wmin,
    wavelength_max=pulse_wmax,
)

# Chop the frames
chopped = frames.chop(choppers.values())

# Propagate the neutrons to the detector
at_sample = chopped.propagate_to(Ltotal)

# Visualize the results
cascade_fig, cascade_ax = at_sample.draw()
../_images/user-guide_odin_simulation_21_0.png

Pipeline#

[16]:

workflow = sl.Pipeline(unwrap.providers(), params=unwrap.params())

workflow[unwrap.PulsePeriod] = sc.reciprocal(ess_beamline.source.frequency)
workflow[unwrap.PulseStride] = 2  # Need for pulse-skipping
workflow[unwrap.SourceTimeRange] = pulse_tmin, pulse_tmax
workflow[unwrap.SourceWavelengthRange] = pulse_wmin, pulse_wmax
workflow[unwrap.Choppers] = choppers

workflow[unwrap.Ltotal] = sample_nexus.coords['Ltotal']
workflow[unwrap.RawData] = sample_nexus

workflow.visualize(unwrap.TofData)

[16]:
../_images/user-guide_odin_simulation_23_0.svg 
[17]:

sample_tofs = workflow.compute(unwrap.TofData)
sample_tofs

[17]:
Show/Hide data repr Show/Hide attributes
scipp.DataArray (19.17 MB)
    • event_time_zero: 3
    • pixel_id: 16384
    • Ltotal
      (pixel_id)
      float64
      m
      60.500, 60.500, ..., 60.500, 60.500
      Values:
      array([60.50001464, 60.50001442, 60.50001419, ..., 60.50001419,
       60.50001442, 60.50001464], shape=(16384,))
    • event_time_zero
      (event_time_zero)
      float64
      µs
      1.730e+15, 1.730e+15, 1.730e+15
      Values:
      array([1.73045043e+15, 1.73045043e+15, 1.73045043e+15])
    • pixel_id
      (pixel_id)
      int64
      𝟙
      0, 1, ..., 16382, 16383
      Values:
      array([    0,     1,     2, ..., 16381, 16382, 16383], shape=(16384,))
    • position
      (pixel_id)
      vector3
      m
      [-2.9765625e-02 -2.9765625e-02 6.0500000e+01], [-2.9296875e-02 -2.9765625e-02 6.0500000e+01], ..., [2.9296875e-02 2.9765625e-02 6.0500000e+01], [2.9765625e-02 2.9765625e-02 6.0500000e+01]
      Values:
      array([[-2.9765625e-02, -2.9765625e-02,  6.0500000e+01],
       [-2.9296875e-02, -2.9765625e-02,  6.0500000e+01],
       [-2.8828125e-02, -2.9765625e-02,  6.0500000e+01],
       ...,
       [ 2.8828125e-02,  2.9765625e-02,  6.0500000e+01],
       [ 2.9296875e-02,  2.9765625e-02,  6.0500000e+01],
       [ 2.9765625e-02,  2.9765625e-02,  6.0500000e+01]], shape=(16384, 3))
    • sample_position
      ()
      vector3
      m
      [ 0. 0. 60.5]
      Values:
      array([ 0. ,  0. , 60.5])
    • source_position
      ()
      vector3
      m
      [0. 0. 0.]
      Values:
      array([0., 0., 0.])
    • (event_time_zero, pixel_id)
      DataArrayView
      binned data [len=5, len=2, ..., len=0, len=2]
      dim='event',
content=DataArray(
          dims=(event: 388566),
          data=float64[counts],
          coords={'sim_wavelength':float64[Å], 'y':float64[m], 'x':float64[m],
                  'event_time_offset':float64[µs], 'tof':float64[µs]})
[18]:

sample_wavs = sample_tofs.transform_coords('wavelength', graph=plane_graph)
sample_wavs

[18]:
Show/Hide data repr Show/Hide attributes
scipp.DataArray (22.13 MB)
    • event_time_zero: 3
    • pixel_id: 16384
    • Ltotal
      (pixel_id)
      float64
      m
      60.500, 60.500, ..., 60.500, 60.500
      Values:
      array([60.50001464, 60.50001442, 60.50001419, ..., 60.50001419,
       60.50001442, 60.50001464], shape=(16384,))
    • event_time_zero
      (event_time_zero)
      float64
      µs
      1.730e+15, 1.730e+15, 1.730e+15
      Values:
      array([1.73045043e+15, 1.73045043e+15, 1.73045043e+15])
    • pixel_id
      (pixel_id)
      int64
      𝟙
      0, 1, ..., 16382, 16383
      Values:
      array([    0,     1,     2, ..., 16381, 16382, 16383], shape=(16384,))
    • position
      (pixel_id)
      vector3
      m
      [-2.9765625e-02 -2.9765625e-02 6.0500000e+01], [-2.9296875e-02 -2.9765625e-02 6.0500000e+01], ..., [2.9296875e-02 2.9765625e-02 6.0500000e+01], [2.9765625e-02 2.9765625e-02 6.0500000e+01]
      Values:
      array([[-2.9765625e-02, -2.9765625e-02,  6.0500000e+01],
       [-2.9296875e-02, -2.9765625e-02,  6.0500000e+01],
       [-2.8828125e-02, -2.9765625e-02,  6.0500000e+01],
       ...,
       [ 2.8828125e-02,  2.9765625e-02,  6.0500000e+01],
       [ 2.9296875e-02,  2.9765625e-02,  6.0500000e+01],
       [ 2.9765625e-02,  2.9765625e-02,  6.0500000e+01]], shape=(16384, 3))
    • sample_position
      ()
      vector3
      m
      [ 0. 0. 60.5]
      Values:
      array([ 0. ,  0. , 60.5])
    • source_position
      ()
      vector3
      m
      [0. 0. 0.]
      Values:
      array([0., 0., 0.])
    • (event_time_zero, pixel_id)
      DataArrayView
      binned data [len=5, len=2, ..., len=0, len=2]
      dim='event',
content=DataArray(
          dims=(event: 388566),
          data=float64[counts],
          coords={'sim_wavelength':float64[Å], 'y':float64[m], 'x':float64[m],
                  'event_time_offset':float64[µs], 'tof':float64[µs], 'wavelength':float64[Å]})

We can now compare our computed wavelengths to the true wavelengths of the neutrons in the McStas simulation:

[19]:

true_wavs = sample_da.hist(sim_wavelength=300).rename(sim_wavelength='wavelength')

pp.plot({
    'true': true_wavs,
    'wfm': sample_wavs.bins.concat().hist(wavelength=true_wavs.coords['wavelength'])
}, title="ODIN McStas simulation")

[19]:
../_images/user-guide_odin_simulation_27_0.svg

Region of interest#

Looking at the counts on the 2d detector panel, we see that there is a central rectangular darker region, surrounded by brighter edges.

[20]:

sample_folded = sample_wavs.bins.concat('event_time_zero').fold(dim='pixel_id', sizes={'y': 128, 'x': 128})
sample_folded.hist().plot(aspect='equal')

[20]:
../_images/user-guide_odin_simulation_29_0.svg

The dark region is where the beam was absorbed by the sample, and this is the region of interest. The brighter edges need to be discarded.

We crop the data using simple array slicing:

[21]:

sel = slice(11, 116, 1)
sample_cropped = sample_folded['y', sel]['x', sel]
sample_cropped.hist().plot(aspect='equal')

[21]:
../_images/user-guide_odin_simulation_31_0.svg

Repeat for the open-beam#

We repeat the conversion to wavelength and crop the edges of the open-beam measurement.

[22]:

# Give the same pixel positions to both sample and open beam.
# Note: this is only because of the way we computed the positions.
# In practice, the geometry should come from the nexus file and this won't be needed.
ob_nexus.coords.update({key: sample_nexus.coords[key] for key in ('position', 'Ltotal')})

workflow[unwrap.Ltotal] = ob_nexus.coords['Ltotal']
workflow[unwrap.RawData] = ob_nexus

ob_tofs = workflow.compute(unwrap.TofData)
ob_wavs = ob_tofs.transform_coords('wavelength', graph=plane_graph)
ob_folded = ob_wavs.bins.concat('event_time_zero').fold(dim='pixel_id', sizes={'y': 128, 'x': 128})
ob_cropped = ob_folded['y', sel]['x', sel]

Normalize the signal#

Finally, we are able to normalize our sample measurement to the open-beam data.

Here, we sum over all pixels before normalizing. There is no spatial structure in the signal, and we are only interested in the wavelength spectrum (where the Bragg edge is). So this is effectively like degrading the detector resolution to a single pixel.

[23]:

# Common set of bins
bins = sc.linspace('wavelength', 1.1, 9.5, 301, unit='angstrom')

num = sample_cropped.bins.concat().hist(wavelength=bins)
den = ob_cropped.bins.concat().hist(wavelength=bins)

# Add variances
num.variances = num.values
den.variances = den.values

normalized = num / den
normalized

[23]:
Show/Hide data repr Show/Hide attributes
scipp.DataArray (8.70 KB)
    • wavelength: 300
    • sample_position
      ()
      vector3
      m
      [ 0. 0. 60.5]
      Values:
      array([ 0. ,  0. , 60.5])
    • source_position
      ()
      vector3
      m
      [0. 0. 0.]
      Values:
      array([0., 0., 0.])
    • wavelength
      (wavelength [bin-edge])
      float64
      Å
      1.1, 1.128, ..., 9.472, 9.5
      Values:
      array([1.1  , 1.128, 1.156, 1.184, 1.212, 1.24 , 1.268, 1.296, 1.324,
       1.352, 1.38 , 1.408, 1.436, 1.464, 1.492, 1.52 , 1.548, 1.576,
       1.604, 1.632, 1.66 , 1.688, 1.716, 1.744, 1.772, 1.8  , 1.828,
       1.856, 1.884, 1.912, 1.94 , 1.968, 1.996, 2.024, 2.052, 2.08 ,
       2.108, 2.136, 2.164, 2.192, 2.22 , 2.248, 2.276, 2.304, 2.332,
       2.36 , 2.388, 2.416, 2.444, 2.472, 2.5  , 2.528, 2.556, 2.584,
       2.612, 2.64 , 2.668, 2.696, 2.724, 2.752, 2.78 , 2.808, 2.836,
       2.864, 2.892, 2.92 , 2.948, 2.976, 3.004, 3.032, 3.06 , 3.088,
       3.116, 3.144, 3.172, 3.2  , 3.228, 3.256, 3.284, 3.312, 3.34 ,
       3.368, 3.396, 3.424, 3.452, 3.48 , 3.508, 3.536, 3.564, 3.592,
       3.62 , 3.648, 3.676, 3.704, 3.732, 3.76 , 3.788, 3.816, 3.844,
       3.872, 3.9  , 3.928, 3.956, 3.984, 4.012, 4.04 , 4.068, 4.096,
       4.124, 4.152, 4.18 , 4.208, 4.236, 4.264, 4.292, 4.32 , 4.348,
       4.376, 4.404, 4.432, 4.46 , 4.488, 4.516, 4.544, 4.572, 4.6  ,
       4.628, 4.656, 4.684, 4.712, 4.74 , 4.768, 4.796, 4.824, 4.852,
       4.88 , 4.908, 4.936, 4.964, 4.992, 5.02 , 5.048, 5.076, 5.104,
       5.132, 5.16 , 5.188, 5.216, 5.244, 5.272, 5.3  , 5.328, 5.356,
       5.384, 5.412, 5.44 , 5.468, 5.496, 5.524, 5.552, 5.58 , 5.608,
       5.636, 5.664, 5.692, 5.72 , 5.748, 5.776, 5.804, 5.832, 5.86 ,
       5.888, 5.916, 5.944, 5.972, 6.   , 6.028, 6.056, 6.084, 6.112,
       6.14 , 6.168, 6.196, 6.224, 6.252, 6.28 , 6.308, 6.336, 6.364,
       6.392, 6.42 , 6.448, 6.476, 6.504, 6.532, 6.56 , 6.588, 6.616,
       6.644, 6.672, 6.7  , 6.728, 6.756, 6.784, 6.812, 6.84 , 6.868,
       6.896, 6.924, 6.952, 6.98 , 7.008, 7.036, 7.064, 7.092, 7.12 ,
       7.148, 7.176, 7.204, 7.232, 7.26 , 7.288, 7.316, 7.344, 7.372,
       7.4  , 7.428, 7.456, 7.484, 7.512, 7.54 , 7.568, 7.596, 7.624,
       7.652, 7.68 , 7.708, 7.736, 7.764, 7.792, 7.82 , 7.848, 7.876,
       7.904, 7.932, 7.96 , 7.988, 8.016, 8.044, 8.072, 8.1  , 8.128,
       8.156, 8.184, 8.212, 8.24 , 8.268, 8.296, 8.324, 8.352, 8.38 ,
       8.408, 8.436, 8.464, 8.492, 8.52 , 8.548, 8.576, 8.604, 8.632,
       8.66 , 8.688, 8.716, 8.744, 8.772, 8.8  , 8.828, 8.856, 8.884,
       8.912, 8.94 , 8.968, 8.996, 9.024, 9.052, 9.08 , 9.108, 9.136,
       9.164, 9.192, 9.22 , 9.248, 9.276, 9.304, 9.332, 9.36 , 9.388,
       9.416, 9.444, 9.472, 9.5  ])
    • (wavelength)
      float64
      𝟙
      0.203, 0.299, ..., nan, nan
      σ = 0.002, 0.002, ..., nan, nan
      Values:
      array([0.20305885, 0.298517  , 0.26657272, 0.23468619, 0.27309398,
       0.25859862, 0.20911155, 0.19287048, 0.22262896, 0.21579718,
       0.22066977, 0.20901057, 0.25778419, 0.23492068, 0.26978935,
       0.30174112, 0.25789055, 0.19465909, 0.22742941, 0.19387524,
       0.23701684, 0.24139434, 0.17726019, 0.23851864, 0.21214185,
       0.20058296, 0.23701192, 0.21973267, 0.2246839 , 0.22673686,
       0.23648984, 0.20615584, 0.19810347, 0.19144866, 0.23542659,
       0.19635902, 0.19977864, 0.22267813, 0.23223037, 0.22675592,
       0.13598439, 0.15963286, 0.1612627 , 0.17513632, 0.18322103,
       0.2877577 , 0.27020968, 0.18514289, 0.17197507, 0.19590591,
       0.22271499, 0.23195811, 0.21803925, 0.1929823 , 0.20716739,
       0.24946071, 0.24546904, 0.19065116, 0.17914902, 0.1868659 ,
       0.17505631, 0.1830747 , 0.18057887, 0.19849887, 0.18383637,
       0.23887873, 0.23452593, 0.1918844 , 0.17905533, 0.1910841 ,
       0.19864238, 0.19010148, 0.17789515, 0.18603461, 0.18273807,
       0.20076333, 0.17672305, 0.19126784, 0.16163583, 0.16926268,
       0.21341869, 0.18230094, 0.16424782, 0.16681085, 0.15389352,
       0.15473417, 0.14923241, 0.15741474, 0.16312388, 0.16086082,
       0.12666257, 0.16915116, 0.15223396, 0.14300318, 0.12823844,
       0.13796595, 0.14674089, 0.11887511, 0.14403659, 0.10191998,
       0.12120476, 0.10131311, 0.13067276, 0.1333081 , 0.1186684 ,
       0.13361134, 0.24023971, 0.32965852, 0.30455242, 0.31352362,
       0.32912617, 0.25780771, 0.29260109, 0.28582712, 0.31070556,
       0.29087978, 0.28538297, 0.32503014, 0.27752369, 0.28088129,
       0.31081308, 0.33021583, 0.30645791, 0.29337918, 0.29597272,
       0.28752214, 0.2888452 , 0.2668881 , 0.31362102, 0.27338785,
       0.2952788 , 0.28913268, 0.30910181, 0.29839627, 0.29444584,
       0.26028107, 0.26012592, 0.27539727, 0.26093925, 0.28724735,
       0.30257028, 0.28892286, 0.25253944, 0.2550615 , 0.27857276,
       0.27005181, 0.27486623, 0.31311787, 0.25343973, 0.24858144,
       0.2529135 , 0.30244067, 0.27809181, 0.27362763, 0.2405646 ,
       0.24111017, 0.28706148, 0.26004099, 0.30621664, 0.25846482,
       0.2524315 , 0.25019702, 0.25244483, 0.26831422, 0.22910217,
       0.28519051, 0.25028803, 0.24192345, 0.23548351, 0.23333206,
       0.24378836, 0.2486244 , 0.2438818 , 0.25757255, 0.24609881,
       0.25059398, 0.25083546, 0.26783868, 0.25518748, 0.24165695,
       0.25689221, 0.22934907, 0.20707947, 0.23493095, 0.23799512,
       0.24201602, 0.24119263, 0.22098713, 0.22910896, 0.22992514,
       0.22194994, 0.23683219, 0.22335878, 0.24688016, 0.2326519 ,
       0.23146235, 0.22675936, 0.21173573, 0.23648312, 0.21927844,
       0.21474325, 0.2131792 , 0.17967469, 0.20694922, 0.20539804,
       0.21448856, 0.2127286 , 0.20335412, 0.21594698, 0.22229512,
       0.20864877, 0.21768901, 0.20946439, 0.19372763, 0.20676283,
       0.21544889, 0.21878148, 0.23102135, 0.21072899, 0.20565189,
       0.21556784, 0.2012182 , 0.21137754, 0.22243013, 0.20368567,
       0.20037298, 0.1982842 , 0.21000778, 0.19987097, 0.20541848,
       0.18450061, 0.18039971, 0.20891713, 0.19252837, 0.19521733,
       0.17340786, 0.19945363, 0.1976459 , 0.19508912, 0.18514961,
       0.17993773, 0.18287986, 0.18815831, 0.19075495, 0.19692821,
       0.18647962, 0.17933189, 0.15515995, 0.19903651,        nan,
       0.19188055, 0.19590059, 0.16248811, 0.18427167, 0.1843403 ,
       0.17652679, 0.18536852, 0.19037717, 0.17309643, 0.1877581 ,
       0.19000313, 0.17676971, 0.1732694 , 0.17835733, 0.18041845,
       0.15968751, 0.17132655, 0.16133463, 0.15584889, 0.16554913,
       0.17727815, 0.17068043, 0.15456427, 0.17129402, 0.17579703,
       0.1560759 , 0.16721622, 0.17051381, 0.15455044, 0.1604986 ,
       0.16573348, 0.16172193, 0.17159208, 0.17024748, 0.15159817,
       0.16531245, 0.15625617, 0.16697528, 0.15485388, 0.15718571,
       0.15059453, 0.16038043, 0.14699499, 0.15503587, 0.13781138,
       0.14903435, 0.13130395,        nan,        nan,        nan])

      Variances (σ²):
      array([3.57878940e-06, 6.03840533e-06, 5.00828802e-06, 4.36535900e-06,
       4.50951138e-06, 4.68190314e-06, 2.95884459e-06, 2.60362703e-06,
       3.25620410e-06, 3.07237660e-06, 3.15737245e-06, 3.09282281e-06,
       4.36853278e-06, 3.91029104e-06, 4.48221435e-06, 5.60566315e-06,
       4.09712523e-06, 2.77736139e-06, 3.41221789e-06, 2.91698606e-06,
       3.86259080e-06, 3.95828031e-06, 2.64013266e-06, 4.09132159e-06,
       3.57479906e-06, 3.32670000e-06, 4.17650643e-06, 3.56811141e-06,
       3.81497821e-06, 3.85324456e-06, 4.40301282e-06, 3.65842849e-06,
       3.26660777e-06, 2.94826350e-06, 4.19882015e-06, 2.95142709e-06,
       2.94767439e-06, 3.34476270e-06, 4.24987852e-06, 3.95347157e-06,
       1.76316459e-06, 2.08219555e-06, 2.35803872e-06, 2.67830572e-06,
       3.00572575e-06, 5.60896965e-06, 5.12954452e-06, 2.98312501e-06,
       3.45132524e-06, 4.61145640e-06, 4.32867519e-06, 3.01679273e-06,
       2.59021924e-06, 2.09443481e-06, 2.05463213e-06, 2.40132014e-06,
       2.07446373e-06, 1.49744394e-06, 1.42942569e-06, 1.42053169e-06,
       1.36279997e-06, 1.39694701e-06, 1.52960813e-06, 1.43191172e-06,
       1.27546259e-06, 1.94690851e-06, 1.85706993e-06, 1.46170375e-06,
       1.24584636e-06, 1.28676880e-06, 1.47738434e-06, 1.45683120e-06,
       1.27952762e-06, 1.37659557e-06, 1.32224349e-06, 1.60014105e-06,
       1.34012626e-06, 1.52555794e-06, 1.04471818e-06, 1.23654375e-06,
       1.77694224e-06, 1.39072510e-06, 1.18499599e-06, 1.27563714e-06,
       1.13894757e-06, 1.15298438e-06, 1.09384408e-06, 1.22266458e-06,
       1.40771536e-06, 1.37062745e-06, 9.96884254e-07, 1.52275533e-06,
       1.24874804e-06, 1.28536229e-06, 1.09282548e-06, 1.32465323e-06,
       1.43678840e-06, 1.13106455e-06, 1.52861460e-06, 1.03418813e-06,
       1.43564571e-06, 1.41285569e-06, 2.37854587e-06, 1.42959074e-06,
       1.21375294e-06, 1.71168211e-06, 2.91760749e-06, 3.79904021e-06,
       3.27425459e-06, 3.18300745e-06, 3.44004751e-06, 2.48705761e-06,
       2.90379577e-06, 2.96211375e-06, 3.35987321e-06, 2.98202659e-06,
       2.78932114e-06, 3.43843844e-06, 2.71611846e-06, 2.85851526e-06,
       3.25437010e-06, 3.76206166e-06, 3.55162209e-06, 3.03414777e-06,
       3.29512710e-06, 2.96691689e-06, 3.18212454e-06, 2.82891321e-06,
       3.77782317e-06, 3.20203788e-06, 3.50827266e-06, 3.32178089e-06,
       3.84969921e-06, 3.76799558e-06, 3.83155761e-06, 3.26710657e-06,
       3.16037958e-06, 3.44323466e-06, 3.19584804e-06, 3.78990594e-06,
       4.12625294e-06, 4.09603578e-06, 3.24223413e-06, 3.39147432e-06,
       4.41275521e-06, 3.94910530e-06, 4.59576209e-06, 5.82331053e-06,
       4.21080333e-06, 4.74701628e-06, 4.93703628e-06, 7.38407808e-06,
       8.40630239e-06, 1.05229699e-05, 6.97426874e-06, 6.21186717e-06,
       7.06732797e-06, 5.28734895e-06, 5.35090192e-06, 3.94179310e-06,
       3.90216193e-06, 3.73505093e-06, 3.83142746e-06, 4.35531308e-06,
       3.44583231e-06, 4.92249605e-06, 3.83542249e-06, 3.93131059e-06,
       3.61667260e-06, 3.67346922e-06, 3.88550760e-06, 3.99103268e-06,
       4.01084331e-06, 4.56985892e-06, 4.31896001e-06, 4.45861112e-06,
       4.40554203e-06, 4.94598122e-06, 4.51519782e-06, 4.21621016e-06,
       4.98688463e-06, 4.25623559e-06, 3.55482929e-06, 4.46680572e-06,
       4.52167396e-06, 4.82631589e-06, 4.89625004e-06, 4.21841449e-06,
       4.75579075e-06, 4.83353338e-06, 4.74727313e-06, 5.51807846e-06,
       4.99119444e-06, 6.12682318e-06, 5.81029665e-06, 6.07172582e-06,
       6.10055408e-06, 5.79709321e-06, 6.93713996e-06, 7.36330124e-06,
       9.33468983e-06, 1.21039254e-05, 1.33626081e-05, 1.58952015e-05,
       8.65595875e-06, 7.11304409e-06, 5.88274859e-06, 5.10433704e-06,
       5.27375124e-06, 5.57091513e-06, 4.90790896e-06, 5.38589339e-06,
       5.28595001e-06, 4.57293756e-06, 5.11959490e-06, 5.33714706e-06,
       5.71430919e-06, 6.05257218e-06, 5.67914749e-06, 5.23675394e-06,
       5.61276032e-06, 5.22067816e-06, 5.75980761e-06, 6.44020108e-06,
       5.96848736e-06, 5.53765047e-06, 5.74192405e-06, 6.18154748e-06,
       5.73017729e-06, 5.98972504e-06, 5.22563030e-06, 5.12746850e-06,
       6.44574179e-06, 5.80328961e-06, 5.84867725e-06, 5.26558881e-06,
       6.73769620e-06, 6.82888984e-06, 6.82340364e-06, 6.38694807e-06,
       6.54506353e-06, 6.77306458e-06, 7.13571229e-06, 7.64171695e-06,
       9.35464238e-06, 1.08150486e-05, 1.34828391e-05, 1.61968696e-05,
       6.05101531e-05,            nan, 2.82204713e-05, 1.67317258e-05,
       9.45255979e-06, 9.75101519e-06, 8.75282206e-06, 8.26851271e-06,
       8.56413462e-06, 9.24039272e-06, 7.62776759e-06, 8.60681553e-06,
       8.84325399e-06, 7.60992278e-06, 7.73648757e-06, 8.02507393e-06,
       7.84474963e-06, 6.85651497e-06, 7.62343609e-06, 6.77307230e-06,
       6.56861758e-06, 7.43567046e-06, 8.30504168e-06, 7.63808244e-06,
       6.97060232e-06, 8.21603862e-06, 8.10562811e-06, 7.02986637e-06,
       8.04457539e-06, 8.53162201e-06, 7.52863929e-06, 8.00594905e-06,
       8.32117624e-06, 8.40720967e-06, 8.98805914e-06, 9.03053623e-06,
       7.97169770e-06, 9.48641657e-06, 8.91415814e-06, 1.06119169e-05,
       1.00119586e-05, 1.02956394e-05, 1.03016197e-05, 1.21231395e-05,
       1.26190039e-05, 1.66842441e-05, 1.94328111e-05, 2.92677477e-05,
       6.74896323e-05,            nan,            nan,            nan])
[24]:

normalized.plot()

[24]:
../_images/user-guide_odin_simulation_36_0.svg

Save the final result#

[25]:

from scippneutron.io import save_xye

to_disk = normalized.copy(deep=False)
to_disk.coords['wavelength'] = sc.midpoints(to_disk.coords['wavelength'])

save_xye('fe_bragg_edge.xye', to_disk)