FREIA detector data#

Visualize detector images, arrival-time distributions, and wavelength spectra. Wavelengths are reconstructed using the WFM chopper settings.

[1]:

%matplotlib widget
import plopp as pp
import scipp as sc

from ess import freia
from ess.freia import data
from ess.reduce.nexus.types import DiskChoppers, Filename, RawDetector
from ess.reduce.unwrap import (
    ChopperFrameSequence,
    DistanceResolution,
    PulsePeriod,
    TimeResolution,
    WavelengthDetector,
)
from ess.reduce.unwrap.types import KeepEventTimeOffset
from ess.reflectometry.types import SampleRun

Select a run#

Download the example sample run and set the resolution used to reconstruct wavelengths. The file is cached locally by pooch (pip install pooch). To use your own data, replace the download function with a file path.

[2]:

freia_mcstas = freia.FreiaMcStasWorkflow(wavelength_from='analytical')
freia_mcstas[Filename[SampleRun]] = data.freia_mcstas_reference_run()
freia_mcstas[KeepEventTimeOffset] = True
freia_mcstas[TimeResolution] = sc.scalar(20.0, unit='us')
freia_mcstas[DistanceResolution] = sc.scalar(0.1, unit='m')

Inspect the WFM chopper cascade#

Inspect the chopper timing used to reconstruct wavelengths. For a non-WFM run, call freia_mcstas.insert(freia.mcstas.non_wfm_choppers) before computing wavelengths. Other configurations can be assigned to freia_mcstas[DiskChoppers[SampleRun]].

[3]:

choppers = freia_mcstas.compute(DiskChoppers[SampleRun])
sc.DataGroup(choppers)

[3]:
  • WBC1
    scippneutron
    DiskChopper
    ()
    DiskChopper(axle_position=<scipp.Variable> () vector3 [m] (-0.0...
  • PSC1
    scippneutron
    DiskChopper
    ()
    DiskChopper(axle_position=<scipp.Variable> () vector3 [m] (-0.0...
  • PSC2
    scippneutron
    DiskChopper
    ()
    DiskChopper(axle_position=<scipp.Variable> () vector3 [m] (-0.0...
  • PSC3
    scippneutron
    DiskChopper
    ()
    DiskChopper(axle_position=<scipp.Variable> () vector3 [m] (-0.0...
  • WBC2
    scippneutron
    DiskChopper
    ()
    DiskChopper(axle_position=<scipp.Variable> () vector3 [m] (-0.0...
  • PSC4
    scippneutron
    DiskChopper
    ()
    DiskChopper(axle_position=<scipp.Variable> () vector3 [m] (-0.1...
  • PSC5
    scippneutron
    DiskChopper
    ()
    DiskChopper(axle_position=<scipp.Variable> () vector3 [m] (-0.1...
  • WBC3
    scippneutron
    DiskChopper
    ()
    DiskChopper(axle_position=<scipp.Variable> () vector3 [m] (-0.1...
[4]:

frames = freia_mcstas.compute(ChopperFrameSequence[SampleRun])
frames.draw()

[4]:

(<Figure size 640x480 with 1 Axes>,
 <Axes: title={'center': 'Frame propagation through chopper cascade'}, xlabel='ms', ylabel='Å'>)

Load and unwrap the detector events#

Compute wavelengths from the event arrival times and chopper transmission bands.

[5]:

results = freia_mcstas.compute((RawDetector[SampleRun], WavelengthDetector[SampleRun]))
raw = results[RawDetector[SampleRun]]
unwrapped = results[WavelengthDetector[SampleRun]]
raw

[5]:
Show/Hide data repr Show/Hide attributes
scipp.DataArray (29.13 MB out of 29.14 MB)
    • pixel_id: 131072
    • height
      (pixel_id)
      float64
      m
      -0.148, -0.148, ..., 0.148, 0.148
      Values:
      array([-0.14765625, -0.14765625, -0.14765625, ...,  0.14765625,
        0.14765625,  0.14765625], shape=(131072,))
    • longitude
      (pixel_id)
      float64
      deg
      -3.998, -3.994, ..., 3.994, 3.998
      Values:
      array([-3.99804688, -3.99414062, -3.99023438, ...,  3.99023438,
        3.99414062,  3.99804688], shape=(131072,))
    • pixel_id
      (pixel_id)
      int64
      0, 1, ..., 131070, 131071
      Values:
      array([     0,      1,      2, ..., 131069, 131070, 131071],
      shape=(131072,))
    • position
      (pixel_id)
      vector3
      m
      [ 0.10790306 -0.98356552 25.75584296], [ 0.10790311 -0.98336149 25.75585721], ..., [-0.18750685 -0.56550938 25.7580494 ], [-0.18750701 -0.56530535 25.75803515]
      Values:
      array([[ 0.10790306, -0.98356552, 25.75584296],
       [ 0.10790311, -0.98336149, 25.75585721],
       [ 0.10790317, -0.98315745, 25.75587145],
       ...,
       [-0.18750669, -0.56571342, 25.75806364],
       [-0.18750685, -0.56550938, 25.7580494 ],
       [-0.18750701, -0.56530535, 25.75803515]], shape=(131072, 3))
    • (pixel_id)
      float64
      counts
      binned data [len=0, len=0, ..., len=0, len=0]
      dim='event',
content=DataArray(
          dims=(event: 553856),
          data=float64[counts],
          coords={'event_time_offset':float64[s], 'event_time_zero':datetime64[ns],
                  'wavelength_from_mcstas':float64[Å]})

Detector image and arrival times#

The image uses longitude and height in the detector’s local frame. Intensities are sums of event weights.

[6]:

detector_image = raw.hist(longitude=120, height=64, dim=raw.dims)
pp.plot(detector_image, norm='log', title='FREIA detector', vmin=1e-1)

[6]:

[7]:

pulse_period = freia_mcstas.compute(PulsePeriod).to(unit='s').value
arrival_times = raw.hist(
    event_time_offset=sc.linspace(
        'event_time_offset', 0.0, pulse_period, 501, unit='s'
    ),
    dim=raw.dims,
)
arrival_times.coords['event_time_offset'] = arrival_times.coords[
    'event_time_offset'
].to(unit='ms')
pp.plot(arrival_times, title='Arrival time within the source period')

[7]:

Wavelengths from analytical frame unwrapping#

The chopper cascade and arrival times determine the wavelengths. Events outside the modeled transmission bands or above the workflow’s relative wavelength uncertainty threshold have NaN wavelengths and do not contribute to the wavelength histograms. Inspect the assigned-event count when assessing the result.

[8]:

valid = sc.isfinite(unwrapped.bins.coords['wavelength'])
print(f'{sc.sum(valid).value:,} / {valid.bins.size().sum().value:,} events have an assigned wavelength')

wavelength_bins = sc.linspace('wavelength', 1.0, 12.0, 441, unit='angstrom')
spectrum = unwrapped.hist(wavelength=wavelength_bins, dim=unwrapped.dims)
pp.plot(spectrum, title='FREIA wavelength spectrum')

548,132 / 553,856 events have an assigned wavelength

[8]:

[9]:

wavelength_image = unwrapped.hist(
    wavelength=wavelength_bins, longitude=120, dim=unwrapped.dims
)
pp.plot(wavelength_image, norm='log', title='Wavelength across the detector', vmin=1e0)

[9]:

Error of the wavelength estimate#

McStas provides the true wavelengths of the neutron events in addition to the time of arrival. We can use that to compare our wavelength estimates to the true wavelength values to figure out the accuracy of the estimates.

It is expected that the wavelength resolution for FREIA in WFM mode should be <2%.

[10]:

unwrapped.bins.coords['wavelength_error'] = (
    (unwrapped.bins.coords['wavelength'] - unwrapped.bins.coords['wavelength_from_mcstas']) /
    unwrapped.bins.coords['wavelength_from_mcstas']
)
unwrapped.hist(wavelength_error=sc.linspace('wavelength_error', -0.02, 0.02, 100), dim=unwrapped.dims).plot()

[10]: