Collimated data reduction for OFFSPEC#

This notebook implements a reduction workflow for reflectometry data collected from the ISIS instrument OFFSPEC using a collimated beam. This workflow implements the same procedure as the corresponding workflow in Mantid, see https://docs.mantidproject.org/nightly/techniques/ISIS_Reflectometry.html.

[1]:

from datetime import datetime
import platform

import scipp as sc
import scippneutron as scn
from orsopy import fileio

import ess
from ess import reflectometry
from ess.external import offspec

[2]:

logger = ess.logging.configure_workflow('offspec_reduction',
                                        filename='offspec.log')

Populating the ORSO header#

We will write the reduced data file following the ORSO `.ort`` standard <https://www.reflectometry.org/file_format/specification>`__, to enable a metadata rich header. We will create an empty header and then populate this.

The data source information#

[6]:

header = fileio.orso.Orso.empty()

header.data_source.owner = fileio.base.Person(
    name="Joshanial F. K. Cooper",
    affiliation="ISIS Neutron and Muon Source",
)
header.data_source.experiment = fileio.data_source.Experiment(
    title="OFFSPEC Sample Data",
    instrument="OFFSPEC",
    start_date="2020-12-14T10:34:02",
    probe="neutron",
    facility="RAL/ISIS/OFFSPEC",
)
header.data_source.sample = fileio.data_source.Sample(
    name="QCS sample",
    category="gas/solid",
    composition="Air | Si(790 A) | Cu(300 A) | SiO2",
)
header.data_source.measurement = fileio.data_source.Measurement(
    instrument_settings=fileio.data_source.InstrumentSettings(
        incident_angle=fileio.base.Value(
            sample.coords["theta"].value, sample.coords["theta"].unit
        ),
        wavelength=None,
        polarization="unpolarized",
    ),
    data_files=[
        offspec.data.sample().rsplit("/", 1)[-1],
        offspec.data.direct_beam().rsplit("/", 1)[-1],
    ],
    scheme="energy-dispersive",
)

The reduction details#

The reduction section can start to be populated also. Entries such as corrections will be filled up through the reduction process.

[7]:

header.reduction.software = fileio.reduction.Software(
    name="ess", version=ess.__version__, platform=platform.platform()
)
header.reduction.timestamp = datetime.now()
header.reduction.creator = fileio.base.Person(
    name="I. D. Scientist",
    affiliation="European Spallation Source",
    contact="i.d.scientist@ess.eu",
)
header.reduction.corrections = []
header.reduction.computer = platform.node()
header.reduction.script = "offspec_mantid.ipynb"

To ensure that the header object is carried through the process, we assign it to the sample scipp.DataArray. The direct beam header object will be overwritten at the normalisation step so we will keep this empty.

[8]:

sample.attrs['orso'] = sc.scalar(header)
direct_beam.attrs['orso'] = sc.scalar(fileio.orso.Orso.empty())

Correcting the position of detector pixels#

The pixel positions in the sample data must be modified to account for the transformation on the detector by rotating it around the sample. We can achieve this by understanding that the sample has been rotated by some amount and that sample measurement has the specular peak at the same pixel as the direct beam measurement has the direct beam. Therefore, we move the sample detector along the arc of the sample rotation by \(2\omega\) (in the OFFSPEC files, \(\omega\) is called 'Theta', which we stored as 'theta' earlier).

[9]:

from scipp.spatial import rotations_from_rotvecs

def pixel_position_correction(data: sc.DataArray):
    rotation = rotations_from_rotvecs(rotation_vectors=sc.vector(value=[-2.0 * data.coords['theta'].value, 0, 0], unit=sc.units.deg))
    return rotation * (data.coords['position'] - data.coords['sample_position'])

[10]:

logger.info("Correcting pixel positions in 'sample.nxs'")
sample.coords['position'] = pixel_position_correction(sample)
sample.attrs['orso'].value.data_source.measurement.comment = 'Pixel positions corrected'

We can visualize the data with a plot. In this plot of sample, we can see the specular intensity at around spectrum numbers 400-410. There is a more intense region, closer to spectrum number 300, which comes from the direct beam of neutrons traveling straight through our sample.

[11]:

sample.hist(tof=50).plot(norm='log')

[11]:

../../../_images/instruments_external_offspec_offspec_reduction_17_0.svg

A region of interest is then defined for the detector. This is defined as twenty-five pixels around the specular peak or the direct beam. The scipp.DataArray is concatenated along the 'spectrum' coordinate at this stage, essentially collapsing all of the events onto a single pixel.

[12]:

sample_roi = sample['spectrum', 389:415].bins.concat('spectrum')
direct_beam_roi = direct_beam['spectrum', 389:415].bins.concat('spectrum')

sample_roi.attrs['orso'].value.reduction.corrections += ['region of interest defined as spectrum 389:415']

The position of these events is then defined as the centre of the region of interest.

[13]:

sample_roi.coords['position'] = sample.coords['position'][401]
direct_beam_roi.coords['position'] = direct_beam.coords['position'][401]

Coordinate transform graph#

To compute the wavelength \(\lambda\), the scattering angle \(\theta\), and the \(Q\) vector for our data we can use a coordinate transform graph. The reflectometry graph is discussed in detail in the Amor reduction notebook and the one used here is nearly identical. The only difference is the Amor instrument uses choppers to define the pulse of neutrons, which is not the case here. The OFFSPEC graph is the standard reflectometry graph, shown below.

[14]:

graph = {**reflectometry.conversions.specular_reflection()}
sc.show_graph(graph, simplified=True)

[14]:

../../../_images/instruments_external_offspec_offspec_reduction_23_0.svg

Computing the wavelength#

The neutron wavelengths can be computed with transform_coords and the graph shown above. We will only use neutrons in the wavelength range of 2 Å to 15.0 Å.

[15]:

wavelength_edges = sc.linspace(dim='wavelength', start=2, stop=15, num=2, unit='angstrom')
sample_wav = reflectometry.conversions.tof_to_wavelength(sample_roi, wavelength_edges,graph=graph)

Since the direct beam measurement is not a reflectometry measurement, we use the no_scatter_graph to convert this to wavelength.

[16]:

no_scatter_graph = {**scn.conversion.graph.beamline.beamline(scatter=False),
                    **scn.conversion.graph.tof.elastic_wavelength(start='tof')}
sc.show_graph(no_scatter_graph, simplified=True)

direct_beam_wav = direct_beam_roi.transform_coords('wavelength', graph=no_scatter_graph)
direct_beam_wav = direct_beam_wav.bin(wavelength=wavelength_edges)

Normalization by monitor#

It is necessary to normalize the sample and direct beam measurements by the summed monitor counts, which accounts for different lengths of measurement and long-timescale natural variation in the pulse. This will ensure that the final data has the correct scaling when the reflectivity data is normalized. First, we convert the data to wavelength, using the no_scatter_graph used previously for the direct beam.

The most reliable monitor for the OFFSPEC instrument is 'monitor2' in the file, therefore this is used.

[17]:

sample_mon_wav = sample_full["monitors"]["monitor2"]["data"].transform_coords(
    "wavelength", graph=no_scatter_graph
)
direct_beam_mon_wav = direct_beam_full["monitors"]["monitor2"]["data"].transform_coords(
    "wavelength", graph=no_scatter_graph
)

A background subtraction is then performed on the monitor data, where the background is taken as any counts at wavelengths greater than 15 Å.

Note

We need to drop the variances of the monitor (using sc.values) because the monitor gets broadcast across all detector pixels. This would introduce correlations in the results and is thus not allowed by Scipp. See Heybrock et al. (2023).

[18]:

wav_min = 2 * sc.Unit('angstrom')
wav_max = 15 * sc.Unit('angstrom')

[19]:

sample_mon_wav -= sc.values(sample_mon_wav['wavelength', wav_max:].mean())
direct_beam_mon_wav -= sc.values(direct_beam_mon_wav['wavelength', wav_max:].mean())
sample_wav.attrs['orso'].value.reduction.corrections += ['monitor background subtraction, background above 15 Å']

The monitors are then summed along the 'wavelength' and this value is used to normalise the data.

[20]:

sample_mon_wav_sum = sample_mon_wav['wavelength', wav_min:wav_max].sum()
direct_beam_mon_wav_sum = direct_beam_mon_wav['wavelength', wav_min:wav_max].sum()
sample_norm = sample_wav / sc.values(sample_mon_wav_sum)
direct_beam_norm = direct_beam_wav / sc.values(direct_beam_mon_wav_sum)
sample_wav.attrs['orso'].value.reduction.corrections += ['normalisation by summed monitor']

Normalisation of sample by direct beam#

The sample and direct beam measurements (which have been normalised by monitor counts) are then histogrammed in wavelength to 100 geometrically spaced points. The histogrammed direct beam is then used to normalised the sample.

Importantly, some relevant metadata (including the ORSO header object) is carried over.

[21]:

edges = sc.geomspace(
    dim="wavelength", start=2, stop=14, num=100, unit=sc.units.angstrom
)
sample_norm_hist = sample_norm.hist(wavelength=edges)
sample_norm_hist.coords.set_aligned('theta', False)
direct_beam_norm_hist = direct_beam_norm.hist(wavelength=edges)
direct_beam_norm_hist.coords.set_aligned('theta', False)

norm_wav = sample_norm_hist / direct_beam_norm_hist
norm_wav.attrs["orso"] = sample_wav.attrs["orso"]
norm_wav.coords["theta"] = sample_wav.coords["theta"]

norm_wav.attrs["orso"].value.reduction.corrections += ["normalisation by direct beam"]

Conversion to \(Q\)#

This normalised data can then be used to compute the reflectivity as a function of the scattering vector \(Q\).

[22]:

norm_q = reflectometry.conversions.theta_to_q(norm_wav, graph=graph)

Which we can visualise.

[23]:

norm_q.plot(norm='log')

[23]:

../../../_images/instruments_external_offspec_offspec_reduction_40_0.svg

Saving the scipp-reduced data as .ort#

We constructed the ORSO header through the reduction process. We can now make use of this when we save our .ort file.

First, we will assume a 3 % of \(Q\) resolution function to be included in our file.

[24]:

norm_q.coords['sigma_Q'] = sc.midpoints(norm_q.coords['Q']) * 0.03

Then, due a bug in orsopy, we need to overwrite the incident angle and wavelength that have been out-populated by the reduction.

[25]:

incident_angle = norm_q.attrs['orso'].value.data_source.measurement.instrument_settings.incident_angle
wavelength = norm_q.attrs['orso'].value.data_source.measurement.instrument_settings.wavelength

norm_q.attrs['orso'].value.data_source.measurement.instrument_settings.wavelength = fileio.base.ValueRange(min=float(wavelength.min), max=float(wavelength.max), unit=wavelength.unit)
norm_q.attrs['orso'].value.data_source.measurement.instrument_settings.incident_angle = fileio.base.Value(magnitude=float(incident_angle.magnitude), unit=incident_angle.unit)

And it is necessary to add the column for our uncertainties, which details the meaning of the uncertainty values we have given.

[26]:

norm_q.attrs['orso'].value.columns.append(fileio.base.ErrorColumn(error_of='R', error_type='uncertainty', value_is='sigma'))
norm_q.attrs['orso'].value.columns.append(fileio.base.ErrorColumn(error_of='Q', error_type='resolution', value_is='sigma'))

Finally, we can save the file.

[27]:

reflectometry.io.save_ort(norm_q, 'offspec.ort')

[28]:

!head offspec.ort

# # ORSO reflectivity data file | 1.0 standard | YAML encoding | https://www.reflectometry.org/
# data_source:
#   owner:
#     name: Joshanial F. K. Cooper
#     affiliation: ISIS Neutron and Muon Source
#   experiment:
#     title: OFFSPEC Sample Data
#     instrument: OFFSPEC
#     start_date: 2020-12-14T10:34:02
#     probe: neutron

Collimated data reduction for OFFSPEC

Contents

Collimated data reduction for OFFSPEC#

Loading the data#

Populating the ORSO header#

The data source information#

The reduction details#

Correcting the position of detector pixels#

Coordinate transform graph#

Computing the wavelength#

Normalization by monitor#

Normalisation of sample by direct beam#

Conversion to \(Q\)#

Saving the scipp-reduced data as .ort#