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.

[1]:

import scipp as sc
from ess.reduce import time_of_flight
from ess import odin
import ess.odin.data  # noqa: F401
from ess.imaging.types import *

Create and configure the workflow#

[2]:

wf = odin.OdinBraggEdgeWorkflow()

wf[Filename[SampleRun]] = odin.data.iron_simulation_sample_small()
wf[NeXusDetectorName] = "event_mode_detectors/timepix3"
wf[time_of_flight.TimeOfFlightLookupTableFilename] = odin.data.odin_tof_lookup_table()

Downloading file 'iron_simulation_sample_small.nxs' from 'https://public.esss.dk/groups/scipp/ess/odin/1/iron_simulation_sample_small.nxs' to '/home/runner/.cache/ess/odin/1'.
Downloading file 'ODIN-tof-lookup-table.h5' from 'https://public.esss.dk/groups/scipp/ess/odin/1/ODIN-tof-lookup-table.h5' to '/home/runner/.cache/ess/odin/1'.

First look at the data#

We load the raw detector data and perform a quick visualization of the event_time_offset spectrum.

[3]:

tmpx3 = wf.compute(DetectorData[SampleRun])
tmpx3

[3]:
Show/Hide data repr Show/Hide attributes
scipp.DataArray (74.89 MB out of 74.91 MB)
    • dim_0: 1024
    • dim_1: 1024
    • detector_number
      (dim_0, dim_1)
      int32
      1, 2, ..., 1048575, 1048576
      Values:
      array([[      1,       2,       3, ...,    1022,    1023,    1024],
       [   1025,    1026,    1027, ...,    2046,    2047,    2048],
       [   2049,    2050,    2051, ...,    3070,    3071,    3072],
       ...,
       [1045505, 1045506, 1045507, ..., 1046526, 1046527, 1046528],
       [1046529, 1046530, 1046531, ..., 1047550, 1047551, 1047552],
       [1047553, 1047554, 1047555, ..., 1048574, 1048575, 1048576]],
      shape=(1024, 1024), dtype=int32)
    • gravity
      ()
      vector3
      m/s^2
      [ 0. -9.80665 0. ]
      Values:
      array([ 0.     , -9.80665,  0.     ])
    • position
      (dim_0, dim_1)
      vector3
      m
      [-0.007033 0.007033 0. ], [-0.007019 0.007033 0. ], ..., [ 0.007019 -0.007033 0. ], [ 0.007033 -0.007033 0. ]
      Values:
      array([[[-0.007033,  0.007033,  0.      ],
        [-0.007019,  0.007033,  0.      ],
        [-0.007006,  0.007033,  0.      ],
        ...,
        [ 0.007006,  0.007033,  0.      ],
        [ 0.007019,  0.007033,  0.      ],
        [ 0.007033,  0.007033,  0.      ]],

       [[-0.007033,  0.007019,  0.      ],
        [-0.007019,  0.007019,  0.      ],
        [-0.007006,  0.007019,  0.      ],
        ...,
        [ 0.007006,  0.007019,  0.      ],
        [ 0.007019,  0.007019,  0.      ],
        [ 0.007033,  0.007019,  0.      ]],

       [[-0.007033,  0.007006,  0.      ],
        [-0.007019,  0.007006,  0.      ],
        [-0.007006,  0.007006,  0.      ],
        ...,
        [ 0.007006,  0.007006,  0.      ],
        [ 0.007019,  0.007006,  0.      ],
        [ 0.007033,  0.007006,  0.      ]],

       ...,

       [[-0.007033, -0.007006,  0.      ],
        [-0.007019, -0.007006,  0.      ],
        [-0.007006, -0.007006,  0.      ],
        ...,
        [ 0.007006, -0.007006,  0.      ],
        [ 0.007019, -0.007006,  0.      ],
        [ 0.007033, -0.007006,  0.      ]],

       [[-0.007033, -0.007019,  0.      ],
        [-0.007019, -0.007019,  0.      ],
        [-0.007006, -0.007019,  0.      ],
        ...,
        [ 0.007006, -0.007019,  0.      ],
        [ 0.007019, -0.007019,  0.      ],
        [ 0.007033, -0.007019,  0.      ]],

       [[-0.007033, -0.007033,  0.      ],
        [-0.007019, -0.007033,  0.      ],
        [-0.007006, -0.007033,  0.      ],
        ...,
        [ 0.007006, -0.007033,  0.      ],
        [ 0.007019, -0.007033,  0.      ],
        [ 0.007033, -0.007033,  0.      ]]], shape=(1024, 1024, 3))
    • sample_position
      ()
      vector3
      m
      [0. 0. 0.]
      Values:
      array([0., 0., 0.])
    • source_position
      ()
      vector3
      m
      [ 0. 0. -60.5]
      Values:
      array([  0. ,   0. , -60.5])
    • x_pixel_offset
      (dim_0, dim_1)
      float32
      m
      -0.007033, -0.007019, ..., 0.007019, 0.007033
      Values:
      array([[-0.007033, -0.007019, -0.007006, ...,  0.007006,  0.007019,
         0.007033],
       [-0.007033, -0.007019, -0.007006, ...,  0.007006,  0.007019,
         0.007033],
       [-0.007033, -0.007019, -0.007006, ...,  0.007006,  0.007019,
         0.007033],
       ...,
       [-0.007033, -0.007019, -0.007006, ...,  0.007006,  0.007019,
         0.007033],
       [-0.007033, -0.007019, -0.007006, ...,  0.007006,  0.007019,
         0.007033],
       [-0.007033, -0.007019, -0.007006, ...,  0.007006,  0.007019,
         0.007033]], shape=(1024, 1024), dtype=float32)
    • y_pixel_offset
      (dim_0, dim_1)
      float32
      m
      0.007033, 0.007033, ..., -0.007033, -0.007033
      Values:
      array([[ 0.007033,  0.007033,  0.007033, ...,  0.007033,  0.007033,
         0.007033],
       [ 0.007019,  0.007019,  0.007019, ...,  0.007019,  0.007019,
         0.007019],
       [ 0.007006,  0.007006,  0.007006, ...,  0.007006,  0.007006,
         0.007006],
       ...,
       [-0.007006, -0.007006, -0.007006, ..., -0.007006, -0.007006,
        -0.007006],
       [-0.007019, -0.007019, -0.007019, ..., -0.007019, -0.007019,
        -0.007019],
       [-0.007033, -0.007033, -0.007033, ..., -0.007033, -0.007033,
        -0.007033]], shape=(1024, 1024), dtype=float32)
    • (dim_0, dim_1)
      float32
      counts
      binned data [len=2, len=0, ..., len=1, len=3]
      dim='event',
content=DataArray(
          dims=(event: 1000000),
          data=float32[counts],
          coords={'event_time_offset':float64[ns], 'event_time_zero':datetime64[ns]})
[4]:

tmpx3.bins.concat().hist(event_time_offset=300).plot()

[4]:
../../_images/user-guide_odin_odin-data-reduction_6_0.svg

Compute neutron time-of-flight/wavelength#

We will now use the workflow to compute the neutron time-of-flight (equivalent to wavelength) using a lookup table built from the beamline chopper information.

[5]:

wf.visualize(DetectorTofData[SampleRun], graph_attr={"rankdir": "LR"})

[5]:
../../_images/user-guide_odin_odin-data-reduction_8_0.svg

Inspect the lookup table#

It is always a good idea to quickly plot the TOF lookup table, as a sanity check.

[6]:

table = wf.compute(time_of_flight.TimeOfFlightLookupTable)
table.plot(figsize=(9, 4))

[6]:
../../_images/user-guide_odin_odin-data-reduction_10_0.svg

Compute neutron wavelengths#

[7]:

sample_wavs = wf.compute(CountsWavelength[SampleRun])

sample_wavs.bins.concat().hist(wavelength=300).plot()

[7]:
../../_images/user-guide_odin_odin-data-reduction_12_0.svg

Process the open-beam run#

We now reuse the same workflow to process the open-beam run that will be used later for normalization.

[8]:

wf[Filename[OpenBeamRun]] = odin.data.iron_simulation_ob_small()
openbeam_wavs = wf.compute(CountsWavelength[OpenBeamRun])

openbeam_wavs.bins.concat().hist(wavelength=300).plot()

Downloading file 'iron_simulation_ob_small.nxs' from 'https://public.esss.dk/groups/scipp/ess/odin/1/iron_simulation_ob_small.nxs' to '/home/runner/.cache/ess/odin/1'.

[8]:
../../_images/user-guide_odin_odin-data-reduction_14_1.svg

Select region of interest by masking outer regions#

Making a 2D histogram of the data shows a dark square region in the centre of the detector panel; this is the region of interest, where the square sample has absorbed neutrons.

[9]:

sample_wavs.hist().plot(aspect='equal')

[9]:
../../_images/user-guide_odin_odin-data-reduction_16_0.svg

The brighter areas around the edges are regions where neutrons did not travel through the sample. We thus want to mask those out using masking rules based on the spatial coordinates of the data:

[10]:

wf[MaskingRules] = {
    'x_pixel_offset': lambda x: (x < sc.scalar(-5.8e-3, unit='m').to(unit=x.unit)) | (x > sc.scalar(5.8e-3, unit='m').to(unit=x.unit)),
    'y_pixel_offset': lambda y: (y < sc.scalar(-5.8e-3, unit='m').to(unit=y.unit)) | (y > sc.scalar(5.8e-3, unit='m').to(unit=y.unit))
}

[11]:

masked = wf.compute(CountsMasked[SampleRun])

masked.hist().plot(aspect='equal')

[11]:
../../_images/user-guide_odin_odin-data-reduction_19_0.svg

Normalize to open beam#

Finally, we use the masked sample and open-beam data to obtain a normalized signal, which reveals the Fe Bragg edges:

[12]:

wbins = sc.linspace('wavelength', 1.1, 9.4, 301, unit='angstrom')

normalized = (
    wf.compute(CountsMasked[SampleRun]).bins.concat().hist(wavelength=wbins) /
    wf.compute(CountsMasked[OpenBeamRun]).bins.concat().hist(wavelength=wbins)
)

normalized.plot()

[12]:
../../_images/user-guide_odin_odin-data-reduction_21_0.svg

Save the final result#

[13]:

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)