From 813e50b927f0e43133f2e6aadef95b4de07ea6ce Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Tue, 10 Oct 2023 14:14:32 +0200 Subject: [PATCH 1/7] Add reflectometry related code from ess --- src/essreflectometry/amor/__init__.py | 7 + src/essreflectometry/amor/beamline.py | 133 +++++++ src/essreflectometry/amor/calibrations.py | 89 +++++ src/essreflectometry/amor/conversions.py | 31 ++ src/essreflectometry/amor/data.py | 33 ++ src/essreflectometry/amor/instrument_view.py | 27 ++ src/essreflectometry/amor/load.py | 160 ++++++++ src/essreflectometry/amor/normalize.py | 44 +++ src/essreflectometry/amor/orso.py | 60 +++ src/essreflectometry/amor/resolution.py | 135 +++++++ src/essreflectometry/amor/tools.py | 112 ++++++ src/essreflectometry/choppers/__init__.py | 16 + src/essreflectometry/choppers/make_chopper.py | 63 +++ src/essreflectometry/choppers/utils.py | 140 +++++++ src/essreflectometry/logging.py | 366 ++++++++++++++++++ .../reflectometry/__init__.py | 6 + .../reflectometry/conversions.py | 251 ++++++++++++ .../reflectometry/corrections.py | 98 +++++ src/essreflectometry/reflectometry/io.py | 41 ++ src/essreflectometry/reflectometry/orso.py | 14 + 20 files changed, 1826 insertions(+) create mode 100644 src/essreflectometry/amor/__init__.py create mode 100644 src/essreflectometry/amor/beamline.py create mode 100644 src/essreflectometry/amor/calibrations.py create mode 100644 src/essreflectometry/amor/conversions.py create mode 100644 src/essreflectometry/amor/data.py create mode 100644 src/essreflectometry/amor/instrument_view.py create mode 100644 src/essreflectometry/amor/load.py create mode 100644 src/essreflectometry/amor/normalize.py create mode 100644 src/essreflectometry/amor/orso.py create mode 100644 src/essreflectometry/amor/resolution.py create mode 100644 src/essreflectometry/amor/tools.py create mode 100644 src/essreflectometry/choppers/__init__.py create mode 100644 src/essreflectometry/choppers/make_chopper.py create mode 100644 src/essreflectometry/choppers/utils.py create mode 100644 src/essreflectometry/logging.py create mode 100644 src/essreflectometry/reflectometry/__init__.py create mode 100644 src/essreflectometry/reflectometry/conversions.py create mode 100644 src/essreflectometry/reflectometry/corrections.py create mode 100644 src/essreflectometry/reflectometry/io.py create mode 100644 src/essreflectometry/reflectometry/orso.py diff --git a/src/essreflectometry/amor/__init__.py b/src/essreflectometry/amor/__init__.py new file mode 100644 index 0000000..a74676e --- /dev/null +++ b/src/essreflectometry/amor/__init__.py @@ -0,0 +1,7 @@ +# SPDX-License-Identifier: BSD-3-Clause +# Copyright (c) 2023 Scipp contributors (https://github.com/scipp) +# flake8: noqa: F401 +from . import calibrations, conversions, data, normalize, resolution, tools +from .beamline import instrument_view_components, make_beamline +from .instrument_view import instrument_view +from .load import load diff --git a/src/essreflectometry/amor/beamline.py b/src/essreflectometry/amor/beamline.py new file mode 100644 index 0000000..2599c79 --- /dev/null +++ b/src/essreflectometry/amor/beamline.py @@ -0,0 +1,133 @@ +# SPDX-License-Identifier: BSD-3-Clause +# Copyright (c) 2023 Scipp contributors (https://github.com/scipp) +import scipp as sc +from scipp.constants import g + +from ..choppers import make_chopper +from ..logging import log_call + + +@log_call( + instrument='amor', message='Constructing AMOR beamline from default parameters' +) +def make_beamline( + sample_rotation: sc.Variable, + beam_size: sc.Variable = None, + sample_size: sc.Variable = None, + detector_spatial_resolution: sc.Variable = None, + gravity: sc.Variable = None, + chopper_frequency: sc.Variable = None, + chopper_phase: sc.Variable = None, + chopper_1_position: sc.Variable = None, + chopper_2_position: sc.Variable = None, +) -> dict: + """ + Amor beamline components. + + Parameters + ---------- + sample_rotation: + Sample rotation (omega) angle. + beam_size: + Size of the beam perpendicular to the scattering surface. Default is `0.001 m`. + sample_size: + Size of the sample in direction of the beam. Default :code:`0.01 m`. + detector_spatial_resolution: + Spatial resolution of the detector. Default is `2.5 mm`. + gravity: + Vector representing the direction and magnitude of the Earth's gravitational + field. Default is `[0, -g, 0]`. + chopper_frequency: + Rotational frequency of the chopper. Default is `6.6666... Hz`. + chopper_phase: + Phase offset between chopper pulse and ToF zero. Default is `-8. degrees of + arc`. + chopper_position: + Position of the chopper. Default is `-15 m`. + + Returns + ------- + : + A dict. + """ + if beam_size is None: + beam_size = 2.0 * sc.units.mm + if sample_size is None: + sample_size = 10.0 * sc.units.mm + if detector_spatial_resolution is None: + detector_spatial_resolution = 0.0025 * sc.units.m + if gravity is None: + gravity = sc.vector(value=[0, -1, 0]) * g + if chopper_frequency is None: + chopper_frequency = sc.scalar(20 / 3, unit='Hz') + if chopper_phase is None: + chopper_phase = sc.scalar(-8.0, unit='deg') + if chopper_1_position is None: + chopper_1_position = sc.vector(value=[0, 0, -15.5], unit='m') + if chopper_2_position is None: + chopper_2_position = sc.vector(value=[0, 0, -14.5], unit='m') + beamline = { + 'sample_rotation': sample_rotation, + 'beam_size': beam_size, + 'sample_size': sample_size, + 'detector_spatial_resolution': detector_spatial_resolution, + 'gravity': gravity, + } + # TODO: in scn.load_nexus, the chopper parameters are stored as coordinates + # of a DataArray, and the data value is a string containing the name of the + # chopper. This does not allow storing e.g. chopper cutout angles. + # We should change this to be a Dataset, which is what we do here. + beamline["source_chopper_2"] = sc.scalar( + make_chopper( + frequency=chopper_frequency, + phase=chopper_phase, + position=chopper_2_position, + ) + ) + beamline["source_chopper_1"] = sc.scalar( + make_chopper( + frequency=chopper_frequency, + phase=chopper_phase, + position=chopper_1_position, + ) + ) + return beamline + + +@log_call(instrument='amor', level='DEBUG') +def instrument_view_components(da: sc.DataArray) -> dict: + """ + Create a dict of instrument view components, containing: + - the sample + - the source chopper + + Parameters + ---------- + da: + The DataArray containing the sample and source chopper coordinates. + + Returns + ------- + : + Dict of instrument view definitions. + """ + return { + "sample": { + 'center': da.meta['sample_position'], + 'color': 'red', + 'size': sc.vector(value=[0.2, 0.01, 0.2], unit=sc.units.m), + 'type': 'box', + }, + "source_chopper_2": { + 'center': da.meta['source_chopper_2'].value["position"].data, + 'color': 'blue', + 'size': sc.vector(value=[0.5, 0, 0], unit=sc.units.m), + 'type': 'disk', + }, + "source_chopper_1": { + 'center': da.meta['source_chopper_1'].value["position"].data, + 'color': 'blue', + 'size': sc.vector(value=[0.5, 0, 0], unit=sc.units.m), + 'type': 'disk', + }, + } diff --git a/src/essreflectometry/amor/calibrations.py b/src/essreflectometry/amor/calibrations.py new file mode 100644 index 0000000..fd3d149 --- /dev/null +++ b/src/essreflectometry/amor/calibrations.py @@ -0,0 +1,89 @@ +# SPDX-License-Identifier: BSD-3-Clause +# Copyright (c) 2023 Scipp contributors (https://github.com/scipp) +import scipp as sc + +from ..reflectometry import orso + + +def supermirror_calibration( + data_array: sc.DataArray, + m_value: sc.Variable = None, + critical_edge: sc.Variable = None, + alpha: sc.Variable = None, +) -> sc.Variable: + """ + Calibrate supermirror measurements + + Parameters + ---------- + data_array: + Data array to get q-bins/values from. + m_value: + m-value for the supermirror. Defaults to 5. + critical_edge: + Supermirror critical edge. Defaults to 0.022 1/angstrom. + alpha: + Supermirror alpha value. Defaults to 0.25 / 0.088 angstrom. + + Returns + ------- + : + Calibrated supermirror measurement. + """ + if m_value is None: + m_value = sc.scalar(5, unit=sc.units.dimensionless) + if critical_edge is None: + critical_edge = 0.022 * sc.Unit('1/angstrom') + if alpha is None: + alpha = sc.scalar(0.25 / 0.088, unit=sc.units.angstrom) + calibration = calibration_factor(data_array, m_value, critical_edge, alpha) + data_array_cal = data_array * calibration + try: + data_array_cal.attrs['orso'].value.reduction.corrections += [ + 'supermirror calibration' + ] + except KeyError: + orso.not_found_warning() + return data_array_cal + + +def calibration_factor( + data_array: sc.DataArray, + m_value: sc.Variable = None, + critical_edge: sc.Variable = None, + alpha: sc.Variable = None, +) -> sc.Variable: + """ + Return the calibration factor for the supermirror. + + Parameters + ---------- + data_array: + Data array to get q-bins/values from. + m_value: + m-value for the supermirror. Defaults to 5. + critical_edge: + Supermirror critical edge. Defaults to 0.022 1/angstrom. + alpha: + Supermirror alpha value. Defaults to 0.25 / 0.088 angstrom. + + Returns + ------- + : + Calibration factor at the midpoint of each Q-bin. + """ + if m_value is None: + m_value = sc.scalar(5, unit=sc.units.dimensionless) + if critical_edge is None: + critical_edge = 0.022 * sc.Unit('1/angstrom') + if alpha is None: + alpha = sc.scalar(0.25 / 0.088, unit=sc.units.angstrom) + q = data_array.coords['Q'] + if data_array.coords.is_edges('Q'): + q = sc.midpoints(q) + max_q = m_value * critical_edge + lim = (q < critical_edge).astype(float) + lim.unit = 'one' + nq = 1.0 / (1.0 - alpha * (q - critical_edge)) + calibration_factor = sc.where(q < max_q, lim + (1 - lim) * nq, sc.scalar(1.0)) + return calibration_factor diff --git a/src/essreflectometry/amor/conversions.py b/src/essreflectometry/amor/conversions.py new file mode 100644 index 0000000..adf16a6 --- /dev/null +++ b/src/essreflectometry/amor/conversions.py @@ -0,0 +1,31 @@ +# SPDX-License-Identifier: BSD-3-Clause +# Copyright (c) 2023 Scipp contributors (https://github.com/scipp) +import scipp as sc + +from ..reflectometry.conversions import specular_reflection as spec_relf_graph + + +def incident_beam( + *, + source_chopper_1: sc.Variable, + source_chopper_2: sc.Variable, + sample_position: sc.Variable +) -> sc.Variable: + """ + Compute the incident beam vector from the source chopper position vector, + instead of the source_position vector. + """ + chopper_midpoint = ( + source_chopper_1.value['position'].data + + source_chopper_2.value['position'].data + ) * sc.scalar(0.5) + return sample_position - chopper_midpoint + + +def specular_reflection() -> dict: + """ + Generate a coordinate transformation graph for Amor reflectometry. + """ + graph = spec_relf_graph() + graph['incident_beam'] = incident_beam + return graph diff --git a/src/essreflectometry/amor/data.py b/src/essreflectometry/amor/data.py new file mode 100644 index 0000000..a6f608f --- /dev/null +++ b/src/essreflectometry/amor/data.py @@ -0,0 +1,33 @@ +# SPDX-License-Identifier: BSD-3-Clause +# Copyright (c) 2023 Scipp contributors (https://github.com/scipp) +_version = '1' + +__all__ = ['get_path'] + + +def _make_pooch(): + import pooch + + return pooch.create( + path=pooch.os_cache('ess/amor'), + env='ESS_AMOR_DATA_DIR', + base_url='https://public.esss.dk/groups/scipp/ess/amor/{version}/', + version=_version, + registry={ + "reference.nxs": "md5:56d493c8051e1c5c86fb7a95f8ec643b", + "sample.nxs": "md5:4e07ccc87b5c6549e190bc372c298e83", + }, + ) + + +_pooch = _make_pooch() + + +def get_path(name: str) -> str: + """ + Return the path to a data file bundled with scippneutron. + + This function only works with example data and cannot handle + paths to custom files. + """ + return _pooch.fetch(name) diff --git a/src/essreflectometry/amor/instrument_view.py b/src/essreflectometry/amor/instrument_view.py new file mode 100644 index 0000000..ec12935 --- /dev/null +++ b/src/essreflectometry/amor/instrument_view.py @@ -0,0 +1,27 @@ +# SPDX-License-Identifier: BSD-3-Clause +# Copyright (c) 2023 Scipp contributors (https://github.com/scipp) +import scipp as sc +import scippneutron as scn + +from .beamline import instrument_view_components + + +def instrument_view( + da: sc.DataArray, components: dict = None, pixel_size: float = 0.0035, **kwargs +): + """ + Instrument view for the Amor instrument, which automatically populates a list of + additional beamline components, and sets the pixel size. + + :param da: The input data for which to display the instrument view. + :param components: A dict of additional components to display. By default, a + set of components defined in `beamline.instrument_view_components()` are added. + :param pixel_size: The detector pixel size. Default is 0.0035. + """ + default_components = instrument_view_components(da) + if components is not None: + default_components = {**default_components, **components} + + return scn.instrument_view( + da, components=default_components, pixel_size=pixel_size, **kwargs + ) diff --git a/src/essreflectometry/amor/load.py b/src/essreflectometry/amor/load.py new file mode 100644 index 0000000..42f27a5 --- /dev/null +++ b/src/essreflectometry/amor/load.py @@ -0,0 +1,160 @@ +# SPDX-License-Identifier: BSD-3-Clause +# Copyright (c) 2023 Scipp contributors (https://github.com/scipp) +from datetime import datetime +from pathlib import Path +from typing import Any, Optional, Union + +import scipp as sc +import scippnexus as snx + +from ..logging import get_logger +from .beamline import make_beamline + + +def _tof_correction(data: sc.DataArray, dim: str = 'tof') -> sc.DataArray: + """ + A correction for the presence of the chopper with respect to the "true" ToF. + Also fold the two pulses. + TODO: generalize mechanism to fold any number of pulses. + + Parameters + ---------- + data: + Input data array to correct. + dim: + Name of the time of flight dimension. + + Returns + ------- + : + ToF corrected data array. + """ + if 'orso' in data.attrs: + data.attrs['orso'].value.reduction.corrections += ['chopper ToF correction'] + tof_unit = data.bins.coords[dim].bins.unit + tau = sc.to_unit( + 1 / (2 * data.coords['source_chopper_2'].value['frequency'].data), + tof_unit, + ) + chopper_phase = data.coords['source_chopper_2'].value['phase'].data + tof_offset = tau * chopper_phase / (180.0 * sc.units.deg) + # Make 2 bins, one for each pulse + edges = sc.concat([-tof_offset, tau - tof_offset, 2 * tau - tof_offset], dim) + data = data.bin({dim: sc.to_unit(edges, tof_unit)}) + # Make one offset for each bin + offset = sc.concat([tof_offset, tof_offset - tau], dim) + # Apply the offset on both bins + data.bins.coords[dim] += offset + # Rebin to exclude second (empty) pulse range + return data.bin({dim: sc.concat([0.0 * sc.units.us, tau], dim)}) + + +def _assemble_event_data(dg: sc.DataGroup) -> sc.DataArray: + """Extract the events as a data array with all required coords. + + Parameters + ---------- + dg: + A data group with the structure of an Amor NeXus file. + + Returns + ------- + : + A data array with the events extracted from ``dg``. + """ + events = dg['instrument']['multiblade_detector'] + events.bins.coords['tof'] = events.bins.coords.pop('event_time_offset') + events.coords['position'] = sc.spatial.as_vectors( + events.coords.pop('x_pixel_offset'), + events.coords.pop('y_pixel_offset'), + events.coords.pop('z_pixel_offset'), + ) + events.coords['sample_position'] = sc.vector([0, 0, 0], unit='m') + return events + + +def _load_nexus_entry(filename: Union[str, Path]) -> sc.DataGroup: + """Load the single entry of a nexus file.""" + with snx.File(filename, 'r') as f: + if len(f.keys()) != 1: + raise snx.NexusStructureError( + f"Expected a single entry in file {filename}, got {len(f.keys())}" + ) + return f['entry'][()] + + +def load( + filename: Union[str, Path], + orso: Optional[Any] = None, + beamline: Optional[dict] = None, +) -> sc.DataArray: + """Load a single Amor data file. + + Parameters + ---------- + filename: + Path of the file to load. + orso: + The orso object to be populated by additional information from the loaded file. + beamline: + A dict defining the beamline parameters. + + Returns + ------- + : + Data array object for Amor dataset. + """ + get_logger('amor').info( + "Loading '%s' as an Amor NeXus file", + filename.filename if hasattr(filename, 'filename') else filename, + ) + full_data = _load_nexus_entry(filename) + data = _assemble_event_data(full_data) + + # Recent versions of scippnexus no longer add variances for events by default, so + # we add them here if they are missing. + if data.bins.constituents['data'].data.variances is None: + data.bins.constituents['data'].data.variances = data.bins.constituents[ + 'data' + ].data.values + + # Convert tof nanoseconds to microseconds for convenience + data.bins.coords['tof'] = data.bins.coords['tof'].to( + unit='us', dtype='float64', copy=False + ) + + # Add beamline parameters + beamline = make_beamline() if beamline is None else beamline + for key, value in beamline.items(): + data.coords[key] = value + + if orso is not None: + populate_orso(orso=orso, data=full_data, filename=filename) + data.attrs['orso'] = sc.scalar(orso) + + # Perform tof correction and fold two pulses + return _tof_correction(data) + + +def populate_orso(orso: Any, data: sc.DataGroup, filename: str) -> Any: + """ + Populate the Orso object, by calling the :code:`base_orso` and adding data from the + file. + + Parameters + ---------- + orso: + The orso object to be populated by additional information from the loaded file. + data: + Data group to source information from. + Should mimic the structure of the NeXus file. + filename: + Path of the file to load. + """ + orso.data_source.experiment.title = data['title'] + orso.data_source.experiment.instrument = data['name'] + orso.data_source.experiment.start_date = datetime.strftime( + datetime.strptime(data['start_time'][:-3], '%Y-%m-%dT%H:%M:%S.%f'), + '%Y-%m-%d', + ) + orso.data_source.measurement.data_files = [filename] diff --git a/src/essreflectometry/amor/normalize.py b/src/essreflectometry/amor/normalize.py new file mode 100644 index 0000000..ceead34 --- /dev/null +++ b/src/essreflectometry/amor/normalize.py @@ -0,0 +1,44 @@ +# SPDX-License-Identifier: BSD-3-Clause +# Copyright (c) 2023 Scipp contributors (https://github.com/scipp) +import scipp as sc + +from ..reflectometry import orso + + +def normalize_by_supermirror( + sample: sc.DataArray, supermirror: sc.DataArray +) -> sc.DataArray: + """ + Normalize the sample measurement by the (ideally calibrated) supermirror. + + Parameters + ---------- + sample: + Sample measurement with coords of 'Q' and 'detector_id'. + supermirror: + Supermirror measurement with coords of 'Q' and 'detector_id', ideally + calibrated. + + Returns + ------- + : + normalized sample. + """ + normalized = sample / supermirror + normalized.masks['no_reference_neutrons'] = (supermirror == sc.scalar(0)).data + try: + normalized.attrs['orso'] = sample.attrs['orso'] + normalized.attrs['orso'].value.reduction.corrections = list( + set( + sample.attrs['orso'].value.reduction.corrections + + supermirror.attrs['orso'].value.reduction.corrections + ) + ) + normalized.attrs[ + 'orso' + ].value.data_source.measurement.reference = supermirror.attrs[ + 'orso' + ].value.data_source.measurement.data_files + except KeyError: + orso.not_found_warning() + return normalized diff --git a/src/essreflectometry/amor/orso.py b/src/essreflectometry/amor/orso.py new file mode 100644 index 0000000..99bb687 --- /dev/null +++ b/src/essreflectometry/amor/orso.py @@ -0,0 +1,60 @@ +# SPDX-License-Identifier: BSD-3-Clause +# Copyright (c) 2023 Scipp contributors (https://github.com/scipp) +import platform +from datetime import datetime + +from orsopy import fileio + +from .. import __version__ + + +def make_orso( + owner: fileio.base.Person, + sample: fileio.data_source.Sample, + creator: fileio.base.Person, + reduction_script: str, +) -> fileio.orso.Orso: + """ + Generate the base Orso object for the Amor instrument. + Populate the Orso object for metadata storage. + + Parameters + ---------- + owner: + The owner of the data set, i.e. the main proposer of the measurement. + sample: + A description of the sample. + creator: + The creator of the reduced data, the person responsible for the + reduction process. + reduction_script: + The script or notebook used for reduction. + + Returns + ------- + : + Orso object with the default parameters for the Amor instrument. + """ + orso = fileio.orso.Orso.empty() + orso.data_source.experiment.probe = 'neutrons' + orso.data_source.experiment.facility = 'Paul Scherrer Institut' + orso.data_source.measurement.scheme = 'angle- and energy-dispersive' + orso.reduction.software = fileio.reduction.Software( + 'scipp-ess', __version__, platform.platform() + ) + orso.reduction.timestep = datetime.now() + orso.reduction.corrections = [] + orso.reduction.computer = platform.node() + orso.columns = [ + fileio.base.Column('Qz', '1/angstrom', 'wavevector transfer'), + fileio.base.Column('R', None, 'reflectivity'), + fileio.base.Column('sR', None, 'standard deivation of reflectivity'), + fileio.base.Column( + 'sQz', '1/angstrom', 'standard deviation of wavevector transfer resolution' + ), + ] + orso.data_source.owner = owner + orso.data_source.sample = sample + orso.reduction.creator = creator + orso.reduction.script = reduction_script + return orso diff --git a/src/essreflectometry/amor/resolution.py b/src/essreflectometry/amor/resolution.py new file mode 100644 index 0000000..7591a8e --- /dev/null +++ b/src/essreflectometry/amor/resolution.py @@ -0,0 +1,135 @@ +# SPDX-License-Identifier: BSD-3-Clause +# Copyright (c) 2023 Scipp contributors (https://github.com/scipp) +import scipp as sc + +from .tools import fwhm_to_std + + +def wavelength_resolution( + chopper_1_position: sc.Variable, + chopper_2_position: sc.Variable, + pixel_position: sc.Variable, +) -> sc.Variable: + """ + Find the wavelength resolution contribution as described in Section 4.3.3 of the + Amor publication (doi: 10.1016/j.nima.2016.03.007). + + Parameters + ---------- + chopper_1_position: + Position of first chopper (the one closer to the source). + chopper_2_position: + Position of second chopper (the one closer to the sample). + pixel_position: + Positions for detector pixels. + + Returns + ------- + : + The angular resolution variable, as standard deviation. + """ + distance_between_choppers = ( + chopper_2_position.data.fields.z - chopper_1_position.data.fields.z + ) + chopper_midpoint = (chopper_1_position.data + chopper_2_position.data) * sc.scalar( + 0.5 + ) + chopper_detector_distance = pixel_position.fields.z - chopper_midpoint.fields.z + return fwhm_to_std(distance_between_choppers / chopper_detector_distance) + + +def sample_size_resolution( + pixel_position: sc.Variable, sample_size: sc.Variable +) -> sc.Variable: + """ + The resolution from the projected sample size, where it may be bigger + than the detector pixel resolution as described in Section 4.3.3 of the Amor + publication (doi: 10.1016/j.nima.2016.03.007). + + Parameters + ---------- + pixel_position: + Positions for detector pixels. + sample_size: + Size of sample. + + Returns + ------- + : + Standard deviation of contribution from the sample size. + """ + return fwhm_to_std( + sc.to_unit(sample_size, 'm') + / sc.to_unit(pixel_position.fields.z, 'm', copy=False) + ) + + +def angular_resolution( + pixel_position: sc.Variable, + theta: sc.Variable, + detector_spatial_resolution: sc.Variable, +) -> sc.Variable: + """ + Determine the angular resolution as described in Section 4.3.3 of the Amor + publication (doi: 10.1016/j.nima.2016.03.007). + + Parameters + ---------- + pixel_position: + Positions for detector pixels. + theta: + Theta values for events. + detector_spatial_resolution: + FWHM of detector pixel resolution. + + Returns + ------- + : + Angular resolution standard deviation + """ + theta_unit = theta.bins.unit if theta.bins is not None else theta.unit + return ( + fwhm_to_std( + sc.to_unit( + sc.atan( + sc.to_unit(detector_spatial_resolution, 'm') + / sc.to_unit(pixel_position.fields.z, 'm', copy=False) + ), + theta_unit, + copy=False, + ) + ) + / theta + ) + + +def sigma_Q( + angular_resolution: sc.Variable, + wavelength_resolution: sc.Variable, + sample_size_resolution: sc.Variable, + q_bins: sc.Variable, +) -> sc.Variable: + """ + Combine all of the components of the resolution and add Q contribution. + + Parameters + ---------- + angular_resolution: + Angular resolution contribution. + wavelength_resolution: + Wavelength resolution contribution. + sample_size_resolution: + Sample size resolution contribution. + q_bins: + Q-bin values. + + Returns + ------- + : + Combined resolution function. + """ + return sc.sqrt( + angular_resolution**2 + + wavelength_resolution**2 + + sample_size_resolution**2 + ).max('detector_number') * sc.midpoints(q_bins) diff --git a/src/essreflectometry/amor/tools.py b/src/essreflectometry/amor/tools.py new file mode 100644 index 0000000..809daca --- /dev/null +++ b/src/essreflectometry/amor/tools.py @@ -0,0 +1,112 @@ +# SPDX-License-Identifier: BSD-3-Clause +# Copyright (c) 2023 Scipp contributors (https://github.com/scipp) +from typing import Union + +import numpy as np +import scipp as sc + +_STD_TO_FWHM = sc.scalar(2.0) * sc.sqrt(sc.scalar(2.0) * sc.log(sc.scalar(2.0))) + + +def fwhm_to_std(fwhm: sc.Variable) -> sc.Variable: + """ + Convert from full-width half maximum to standard deviation. + + Parameters + ---------- + fwhm: + Full-width half maximum. + + Returns + ------- + : + Standard deviation. + """ + # Enables the conversion from full width half + # maximum to standard deviation + return fwhm / _STD_TO_FWHM + + +def std_to_fwhm(std: sc.Variable) -> sc.Variable: + """ + Convert from standard deviation to full-width half maximum. + + Parameters + ---------- + std: + Standard deviation. + + Returns + ------- + : + Full-width half maximum. + """ + # Enables the conversion from full width half + # maximum to standard deviation + return std * _STD_TO_FWHM + + +def linlogspace( + dim: str, + edges: Union[list, np.ndarray], + scale: Union[list, str], + num: Union[list, int], + unit: str = None, +) -> sc.Variable: + """ + Generate a 1d array of bin edges with a mixture of linear and/or logarithmic + spacings. + + Examples: + + - Create linearly spaced edges (equivalent to `sc.linspace`): + linlogspace(dim='x', edges=[0.008, 0.08], scale='linear', num=50, unit='m') + - Create logarithmically spaced edges (equivalent to `sc.geomspace`): + linlogspace(dim='x', edges=[0.008, 0.08], scale='log', num=50, unit='m') + - Create edges with a linear and a logarithmic part: + linlogspace(dim='x', edges=[1, 3, 8], scale=['linear', 'log'], num=[16, 20]) + + Parameters + ---------- + dim: + The dimension of the output Variable. + edges: + The edges for the different parts of the mesh. + scale: + A string or list of strings specifying the scaling for the different + parts of the mesh. Possible values for the scaling are `"linear"` and `"log"`. + If a list is supplied, the length of the list must be one less than the length + of the `edges` parameter. + num: + An integer or a list of integers specifying the number of points to use + in each part of the mesh. If a list is supplied, the length of the list must be + one less than the length of the `edges` parameter. + unit: + The unit of the output Variable. + + Returns + ------- + : + Lin-log spaced Q-bin edges. + """ + if not isinstance(scale, list): + scale = [scale] + if not isinstance(num, list): + num = [num] + if len(scale) != len(edges) - 1: + raise ValueError( + "Sizes do not match. The length of edges should be one " + "greater than scale." + ) + + funcs = {"linear": sc.linspace, "log": sc.geomspace} + grids = [] + for i in range(len(edges) - 1): + # Skip the leading edge in the piece when concatenating + start = int(i > 0) + mesh = funcs[scale[i]]( + dim=dim, start=edges[i], stop=edges[i + 1], num=num[i] + start, unit=unit + ) + grids.append(mesh[dim, start:]) + + return sc.concat(grids, dim) diff --git a/src/essreflectometry/choppers/__init__.py b/src/essreflectometry/choppers/__init__.py new file mode 100644 index 0000000..e5b97f9 --- /dev/null +++ b/src/essreflectometry/choppers/__init__.py @@ -0,0 +1,16 @@ +# SPDX-License-Identifier: BSD-3-Clause +# Copyright (c) 2023 Scipp contributors (https://github.com/scipp) + +# flake8: noqa F401 + +from .make_chopper import make_chopper +from .utils import ( + angular_frequency, + cutout_angles_begin, + cutout_angles_center, + cutout_angles_end, + cutout_angles_width, + find_chopper_keys, + time_closed, + time_open, +) diff --git a/src/essreflectometry/choppers/make_chopper.py b/src/essreflectometry/choppers/make_chopper.py new file mode 100644 index 0000000..f2747dc --- /dev/null +++ b/src/essreflectometry/choppers/make_chopper.py @@ -0,0 +1,63 @@ +# SPDX-License-Identifier: BSD-3-Clause +# Copyright (c) 2023 Scipp contributors (https://github.com/scipp) +import scipp as sc + +from . import utils + + +def make_chopper( + frequency: sc.Variable, + position: sc.Variable, + phase: sc.Variable = None, + cutout_angles_center: sc.Variable = None, + cutout_angles_width: sc.Variable = None, + cutout_angles_begin: sc.Variable = None, + cutout_angles_end: sc.Variable = None, + kind: str = None, +) -> sc.Dataset: + """ + Create a Dataset that holds chopper parameters. + This ensures the Dataset is compatible with the other functions in the choppers + module. + Defining a chopper's cutout angles can either constructed from an array of cutout + centers and an array of angular widths, or two arrays containing the begin (leading) + and end (closing) angles of the cutout windows. + + :param frequency: The rotational frequency of the chopper. + :param position: The position vector of the chopper. + :param phase: The chopper phase. + :param cutout_angles_center: The centers of the chopper cutout angles. + :param cutout_angles_width: The angular widths of the chopper cutouts. + :param cutout_angles_begin: The starting/opening angles of the chopper cutouts. + :param cutout_angles_end: The ending/closing angles of the chopper cutouts. + :param kind: The chopper kind. Any string can be supplied, but WFM choppers should + be given 'wfm' as their kind. + """ + data = {"frequency": frequency, "position": position} + if phase is not None: + data["phase"] = phase + if cutout_angles_center is not None: + data["cutout_angles_center"] = cutout_angles_center + if cutout_angles_width is not None: + data["cutout_angles_width"] = cutout_angles_width + if cutout_angles_begin is not None: + data["cutout_angles_begin"] = cutout_angles_begin + if cutout_angles_end is not None: + data["cutout_angles_end"] = cutout_angles_end + if kind is not None: + data["kind"] = kind + chopper = sc.Dataset(data=data) + + # Sanitize input parameters + if (None not in [cutout_angles_begin, cutout_angles_end]) or ( + None not in [cutout_angles_center, cutout_angles_width] + ): + widths = utils.cutout_angles_width(chopper) + if (sc.min(widths) < sc.scalar(0.0, unit=widths.unit)).value: + raise ValueError("Negative window width found in chopper cutout angles.") + if not sc.allsorted(utils.cutout_angles_begin(chopper), dim=widths.dim): + raise ValueError("Chopper begin cutout angles are not monotonic.") + if not sc.allsorted(utils.cutout_angles_end(chopper), dim=widths.dim): + raise ValueError("Chopper end cutout angles are not monotonic.") + + return chopper diff --git a/src/essreflectometry/choppers/utils.py b/src/essreflectometry/choppers/utils.py new file mode 100644 index 0000000..483c20d --- /dev/null +++ b/src/essreflectometry/choppers/utils.py @@ -0,0 +1,140 @@ +# SPDX-License-Identifier: BSD-3-Clause +# Copyright (c) 2023 Scipp contributors (https://github.com/scipp) +import scipp as sc +from scipp.constants import pi + + +def cutout_angles_begin(chopper: sc.Dataset, unit="rad") -> sc.Variable: + """ + Get the starting/opening angles of the chopper cutouts. + + :param chopper: The Dataset containing the chopper parameters. + :param unit: Convert to this unit before returning. Default is `'rad'`. + """ + if "cutout_angles_begin" in chopper: + out = chopper["cutout_angles_begin"].data + elif all(x in chopper for x in ["cutout_angles_width", "cutout_angles_center"]): + out = ( + chopper["cutout_angles_center"].data + - 0.5 * chopper["cutout_angles_width"].data + ) + else: + raise KeyError( + "Chopper does not contain the information required to compute " + "the cutout_angles_begin." + ) + return sc.to_unit(out, unit, copy=False) + + +def cutout_angles_end(chopper: sc.Dataset, unit="rad") -> sc.Variable: + """ + Get the ending/closing angles of the chopper cutouts. + + :param chopper: The Dataset containing the chopper parameters. + :param unit: Convert to this unit before returning. Default is `'rad'`. + """ + if "cutout_angles_end" in chopper: + out = chopper["cutout_angles_end"].data + elif all(x in chopper for x in ["cutout_angles_width", "cutout_angles_center"]): + out = ( + chopper["cutout_angles_center"].data + + 0.5 * chopper["cutout_angles_width"].data + ) + else: + raise KeyError( + "Chopper does not contain the information required to compute " + "the cutout_angles_end." + ) + return sc.to_unit(out, unit, copy=False) + + +def cutout_angles_width(chopper: sc.Dataset, unit="rad") -> sc.Variable: + """ + Get the angular widths of the chopper cutouts. + + :param chopper: The Dataset containing the chopper parameters. + :param unit: Convert to this unit before returning. Default is `'rad'`. + """ + if "cutout_angles_width" in chopper: + out = chopper["cutout_angles_width"].data + elif all(x in chopper for x in ["cutout_angles_begin", "cutout_angles_end"]): + out = chopper["cutout_angles_end"].data - chopper["cutout_angles_begin"].data + else: + raise KeyError( + "Chopper does not contain the information required to compute " + "the cutout_angles_width." + ) + return sc.to_unit(out, unit, copy=False) + + +def cutout_angles_center(chopper: sc.Dataset, unit="rad") -> sc.Variable: + """ + Get the angular centers of the chopper cutouts. + + :param chopper: The Dataset containing the chopper parameters. + :param unit: Convert to this unit before returning. Default is `'rad'`. + """ + if "cutout_angles_center" in chopper: + out = chopper["cutout_angles_center"].data + elif all(x in chopper for x in ["cutout_angles_begin", "cutout_angles_end"]): + out = 0.5 * ( + chopper["cutout_angles_begin"].data + chopper["cutout_angles_end"].data + ) + else: + raise KeyError( + "Chopper does not contain the information required to compute " + "the cutout_angles_center." + ) + return sc.to_unit(out, unit, copy=False) + + +def angular_frequency(chopper: sc.Dataset) -> sc.Variable: + """ + Get the angular frequency of the chopper. + + :param chopper: The Dataset containing the chopper parameters. + """ + return (2.0 * sc.units.rad) * pi * chopper["frequency"].data + + +def time_open(chopper: sc.Dataset, unit: str = "us") -> sc.Variable: + """ + Get the times when a chopper window is open. + + :param chopper: The Dataset containing the chopper parameters. + :param unit: Convert to this unit before returning. Default is `'rad'`. + """ + return sc.to_unit( + (cutout_angles_begin(chopper) + sc.to_unit(chopper["phase"].data, "rad")) + / angular_frequency(chopper), + unit, + copy=False, + ) + + +def time_closed(chopper: sc.Dataset, unit: str = "us") -> sc.Variable: + """ + Get the times when a chopper window is closed. + + :param chopper: The Dataset containing the chopper parameters. + :param unit: Convert to this unit before returning. Default is `'rad'`. + """ + return sc.to_unit( + (cutout_angles_end(chopper) + sc.to_unit(chopper["phase"].data, "rad")) + / angular_frequency(chopper), + unit, + copy=False, + ) + + +def find_chopper_keys(da: sc.DataArray) -> list: + """ + Scan the coords of the data container and return a list of all entries starting with + `"chopper"`, ignoring case. + TODO: This is a very brittle mechanism. In Nexus, choppers are identified with + `NXdisk_chopper`. We could make use of these identifiers instead, but we currently + do not store the NX attributes when loading files in `scn.load_nexus()`. + + :param da: The DataArray containing the coordinates to scan. + """ + return [key for key in da.coords if key.lower().startswith("chopper")] diff --git a/src/essreflectometry/logging.py b/src/essreflectometry/logging.py new file mode 100644 index 0000000..5319687 --- /dev/null +++ b/src/essreflectometry/logging.py @@ -0,0 +1,366 @@ +# SPDX-License-Identifier: BSD-3-Clause +# Copyright (c) 2023 Scipp contributors (https://github.com/scipp) +# @author Jan-Lukas Wynen +""" +Utilities for logging. + +All code in the ess package should log auxiliary information +through this module. +Loggers should be obtained through :py:func:`ess.logging.get_logger` +and pass a name that reflects the current context. +E.g., in the loki package, pass ``'loki'`` (all lowercase). + +Logging can be configured using :py:func:`ess.logging.configure` or +:py:func:`ess.logging.configure_workflow`. +Use the latter at the beginning of workflow notebooks. +""" + +import functools +import inspect +import logging +import logging.config +from copy import copy +from os import PathLike +from typing import Any, Callable, List, Optional, Sequence, Union + +import scipp as sc +import scippneutron as scn +from scipp.utils import running_in_jupyter + + +def get_logger(subname: Optional[str] = None) -> logging.Logger: + """Return one of ess's loggers. + + Parameters + ---------- + subname: + Name of an instrument, technique, or workflow. + If given, return the logger with the given name + as a child of the ess logger. + Otherwise, return the general ess logger. + + Returns + ------- + : + The requested logger. + """ + name = 'scipp.ess' + ('.' + subname if subname else '') + return logging.getLogger(name) + + +def log_call( + *, instrument: str, message: str = None, level: Union[int, str] = logging.INFO +): + """ + Decorator that logs a message every time the function is called. + """ + level = logging.getLevelName(level) if isinstance(level, str) else level + + def deco(f: Callable): + @functools.wraps(f) + def impl(*args, **kwargs): + if message is not None: + get_logger(instrument).log(level, message) + else: + get_logger(instrument).log(level, 'Calling %s', _function_name(f)) + return f(*args, **kwargs) + + return impl + + return deco + + +class Formatter(logging.Formatter): + """ + Logging formatter that indents messages and optionally shows threading information. + """ + + def __init__(self, show_thread: bool, show_process: bool): + """ + Initialize the formatter. + + The formatting is mostly fixed. + Only printing of thread and processor names can be toggled using the + corresponding arguments. + Times are always printed in ISO 8601 format. + """ + fmt_proc = '%(processName)s' if show_process else '' + fmt_thread = '%(threadName)s' if show_thread else '' + if show_process: + fmt_proc_thread = fmt_proc + if show_thread: + fmt_proc_thread += ',' + fmt_thread + elif show_thread: + fmt_proc_thread = fmt_thread + else: + fmt_proc_thread = '' + fmt_pre = '[%(asctime)s] %(levelname)-8s ' + fmt_post = '<%(name)s> : %(message)s' + fmt = ( + fmt_pre + + ('{' + fmt_proc_thread + '} ' if fmt_proc_thread else '') + + fmt_post + ) + super().__init__(fmt, datefmt='%Y-%m-%dT%H:%M:%S%z') + + def format(self, record: logging.LogRecord) -> str: + record = copy(record) + record.msg = '\n ' + record.msg.replace('\n', '\n ') + return super().format(record) + + +def default_loggers_to_configure() -> List[logging.Logger]: + """ + Return a list of all loggers that get configured by ess by default. + """ + import pooch + + return [ + sc.get_logger(), + logging.getLogger('Mantid'), + pooch.get_logger(), + ] + + +def configure( + *, + filename: Optional[Union[str, PathLike]] = 'scipp.ess.log', + file_level: Union[str, int] = logging.INFO, + stream_level: Union[str, int] = logging.WARNING, + widget_level: Union[str, int] = logging.INFO, + show_thread: bool = False, + show_process: bool = False, + loggers: Optional[Sequence[Union[str, logging.Logger]]] = None, +): + """Set up logging for the ess package. + + This function is meant as a helper for application (or notebook) developers. + It configures the loggers of ess, scippneutron, scipp, and some + third party packages. + *Calling it from a library should be avoided* + because it can mess up a user's setup. + + Up to 3 handlers are configured: + + - *File Handler* Writes files to a file with given name. + Can be disabled using the `filename` argument. + - *Stream Handler* Writes to `sys.stderr`. + - *Widget Handler* Writes to a :py:class:`scipp.logging.LogWidget` + if Python is running in a Jupyter notebook. + + Parameters + ---------- + filename: + Name of the log file. Overwrites existing files. + Setting this to `None` disables logging to file. + file_level: + Log level for the file handler. + stream_level: + Log level for the stream handler. + widget_level: + Log level for the widget handler. + show_thread: + If `True`, log messages include the name of the thread + the message originates from. + show_process: + If `True`, log messages include the name of the process + the message originates from. + loggers: + Collection of loggers or names of loggers to configure. + If not given, uses :py:func:`default_loggers_to_configure`. + + See Also + -------- + ess.logging.configure_workflow: + Configure logging and do some additional setup for a reduction workflow. + """ + if configure.is_configured: + get_logger().warning( + 'Called `logging.configure` but logging is already configured' + ) + return + + handlers = _make_handlers( + filename, file_level, stream_level, widget_level, show_thread, show_process + ) + base_level = _base_level([file_level, stream_level, widget_level]) + loggers = { + logging.getLogger(logger) if isinstance(logger, str) else logger + for logger in (default_loggers_to_configure() if loggers is None else loggers) + } + for logger in loggers: + _configure_logger(logger, handlers, base_level) + if any(logger.name == 'Mantid' for logger in loggers): + _configure_mantid_logging('notice') + + configure.is_configured = True + + +configure.is_configured = False + + +def configure_workflow( + workflow_name: Optional[str] = None, *, display: Optional[bool] = None, **kwargs +) -> logging.Logger: + """Configure logging for a reduction workflow. + + Configures loggers, logs a greeting message, sets up a logger for a workflow, + and optionally creates and displays a log widget. + + Parameters + ---------- + workflow_name: + Used as the name of the returned logger. + display: + If `True`, show a :py:class:`scipp.logging.LogWidget` + in the outputs of the current cell. + Defaults to `True` in Jupyter and `False` otherwise. + kwargs: + Forwarded to :py:func:`ess.logging.configure`. + Refer to that function for details. + + Returns + ------- + : + A logger for use in the workflow. + + See Also + -------- + ess.logging.configure: + General logging setup. + """ + configure(**kwargs) + greet() + if (display is None and running_in_jupyter()) or display: + sc.display_logs() + return get_logger(workflow_name) + + +def greet(): + """Log a message showing the versions of important packages.""" + # Import here so we don't import from a partially built package. + from . import __version__ + + msg = f'''Software Versions: + ess: {__version__} (https://scipp.github.io/ess) + scippneutron: {scn.__version__} (https://scipp.github.io/scippneutron) + scipp: {sc.__version__} (https://scipp.github.io)''' + mantid_version = _mantid_version() + if mantid_version: + msg += f'\n Mantid: {mantid_version} (https://www.mantidproject.org)' + get_logger().info(msg) + + +_INSTRUMENTS = [ + 'amor', + 'beer', + 'bifrost', + 'cspec', + 'dream', + 'estia', + 'freia', + 'heimdal', + 'loki', + 'magic', + 'miracles', + 'nmx', + 'odin', + 'skadi', + 'trex', + 'v20', + 'vespa', +] + + +def _deduce_instrument_name(f: Any) -> Optional[str]: + # Assumes package name: ess.[.subpackage] + package = inspect.getmodule(f).__package__ + components = package.split('.', 2) + try: + if components[0] == 'ess': + candidate = components[1] + if candidate in _INSTRUMENTS: + return candidate + except IndexError: + pass + return None + + +def _function_name(f: Callable) -> str: + if hasattr(f, '__module__'): + return f'{f.__module__}.{f.__name__}' + return f.__name__ + + +def _make_stream_handler( + level: Union[str, int], show_thread: bool, show_process: bool +) -> logging.StreamHandler: + handler = logging.StreamHandler() + handler.setLevel(level) + handler.setFormatter(Formatter(show_thread, show_process)) + return handler + + +def _make_file_handler( + filename: Union[str, PathLike], + level: Union[str, int], + show_thread: bool, + show_process: bool, +) -> logging.FileHandler: + handler = logging.FileHandler(filename, mode='w') + handler.setLevel(level) + handler.setFormatter(Formatter(show_thread, show_process)) + return handler + + +def _make_handlers( + filename: Optional[Union[str, PathLike]], + file_level: Union[str, int], + stream_level: Union[str, int], + widget_level: Union[str, int], + show_thread: bool, + show_process: bool, +) -> List[logging.Handler]: + handlers = [_make_stream_handler(stream_level, show_thread, show_process)] + if filename is not None: + handlers.append( + _make_file_handler(filename, file_level, show_thread, show_process) + ) + if running_in_jupyter(): + handlers.append(sc.logging.make_widget_handler()) + return handlers + + +def _configure_logger( + logger: logging.Logger, handlers: List[logging.Handler], level: Union[str, int] +): + for handler in handlers: + logger.addHandler(handler) + logger.setLevel(level) + + +def _configure_mantid_logging(level: str): + try: + from mantid.utils.logging import log_to_python + + log_to_python(level) + except ImportError: + pass + + +def _base_level(levels: List[Union[str, int]]) -> int: + return min( + ( + logging.getLevelName(level) if isinstance(level, str) else level + for level in levels + ) + ) + + +def _mantid_version() -> Optional[str]: + try: + import mantid + + return mantid.__version__ + except ImportError: + return None diff --git a/src/essreflectometry/reflectometry/__init__.py b/src/essreflectometry/reflectometry/__init__.py new file mode 100644 index 0000000..fcb6518 --- /dev/null +++ b/src/essreflectometry/reflectometry/__init__.py @@ -0,0 +1,6 @@ +# SPDX-License-Identifier: BSD-3-Clause +# Copyright (c) 2023 Scipp contributors (https://github.com/scipp) + +# flake8: noqa: F401 + +from . import conversions, corrections, io diff --git a/src/essreflectometry/reflectometry/conversions.py b/src/essreflectometry/reflectometry/conversions.py new file mode 100644 index 0000000..bbabb27 --- /dev/null +++ b/src/essreflectometry/reflectometry/conversions.py @@ -0,0 +1,251 @@ +# SPDX-License-Identifier: BSD-3-Clause +# Copyright (c) 2023 Scipp contributors (https://github.com/scipp) +import scipp as sc +from scipp.constants import h, m_n, pi +from scippneutron._utils import elem_dtype, elem_unit +from scippneutron.conversion.graph import beamline, tof + +from . import orso + + +def theta( + gravity: sc.Variable, + wavelength: sc.Variable, + incident_beam: sc.Variable, + scattered_beam: sc.Variable, + sample_rotation: sc.Variable, +) -> sc.Variable: + """ + Compute the theta angle, including gravity correction, + This is similar to the theta calculation in SANS (see + https://docs.mantidproject.org/nightly/algorithms/Q1D-v2.html#q-unit-conversion), + but we ignore the horizontal `x` component. + See the schematic in Fig 5 of doi: 10.1016/j.nima.2016.03.007. + + Parameters + ---------- + gravity: + The three-dimensional vector describing gravity. + wavelength: + Wavelength values for the events. + incident_beam: + Vector for incident beam. + scatter_beam: + Vector for scattered beam. + sample_rotation: + Rotation of sample wrt to incoming beam. + + Returns + ------- + : + Theta values, accounting for gravity. + """ + grav = sc.norm(gravity) + L2 = sc.norm(scattered_beam) + y = sc.dot(scattered_beam, gravity) / grav + y_correction = sc.to_unit(wavelength, elem_unit(L2), copy=True) + y_correction *= y_correction + drop = L2**2 + drop *= grav * (m_n**2 / (2 * h**2)) + # Optimization when handling either the dense or the event coord of binned data: + # - For the event coord, both operands have same dims, and we can multiply in place + # - For the dense coord, we need to broadcast using non in-place operation + if set(drop.dims).issubset(set(y_correction.dims)): + y_correction *= drop + else: + y_correction = y_correction * drop + y_correction += y + out = sc.abs(y_correction, out=y_correction) + out /= L2 + out = sc.asin(out, out=out) + out -= sc.to_unit(sample_rotation, 'rad') + return out + + +def reflectometry_q(wavelength: sc.Variable, theta: sc.Variable) -> sc.Variable: + """ + Compute the Q vector from the theta angle computed as the difference + between gamma and omega. + Note that this is identical the 'normal' Q defined in scippneutron, except that + the `theta` angle is given as an input instead of `two_theta`. + + Parameters + ---------- + wavelength: + Wavelength values for the events. + theta: + Theta values, accounting for gravity. + + Returns + ------- + : + Q-values. + """ + dtype = elem_dtype(wavelength) + c = (4 * pi).astype(dtype) + return c * sc.sin(theta.astype(dtype, copy=False)) / wavelength + + +def specular_reflection() -> dict: + """ + Generate a coordinate transformation graph for specular reflection reflectometry. + + Returns + ------- + : + Specular reflectometry graph. + """ + graph = { + **beamline.beamline(scatter=True), + **tof.elastic_wavelength("tof"), + "theta": theta, + "Q": reflectometry_q, + } + return graph + + +def tof_to_wavelength( + data_array: sc.DataArray, wavelength_edges: sc.Variable = None, graph: dict = None +) -> sc.DataArray: + """ + Use :code:`transform_coords` to convert from ToF to wavelength, cutoff high and + low limits for wavelength, and add necessary ORSO metadata. + + Parameters + ---------- + data_array: + Data array to convert. + wavelength_edges: + The lower and upper limits for the wavelength. + If :code:`None`, no binning is performed. + graph: + Graph for :code:`transform_coords`. + + Returns + ------- + : + New data array with wavelength dimension. + """ + graph = graph if graph is not None else specular_reflection() + data_array_wav = data_array.transform_coords(["wavelength"], graph=graph) + if wavelength_edges is not None: + data_array_wav = data_array_wav.bin({wavelength_edges.dim: wavelength_edges}) + try: + from orsopy import fileio + + unit = data_array_wav.coords['wavelength'].unit + # This insures that when the unit is Å it is written as + # angstrom in the ORSO object. + if unit == 'angstrom': + unit = 'angstrom' + orso_measurement = data_array_wav.attrs['orso'].value.data_source.measurement + orso_measurement.instrument_settings.wavelength = fileio.base.ValueRange( + float(data_array_wav.coords['wavelength'].min().value), + float(data_array_wav.coords['wavelength'].max().value), + unit, + ) + except ImportError: + orso.not_found_warning() + return data_array_wav + + +def wavelength_to_theta( + data_array: sc.DataArray, theta_edges: sc.Variable = None, graph: dict = None +) -> sc.DataArray: + """ + Use :code:`transform_coords` to find the theta values for the events and + potentially add ORSO metadata. + + Parameters + ---------- + data_array: + Data array to convert. + theta_edges: + The lower and upper limits for the theta. If :code:`None`, no + binning is performed. + graph: + Graph for :code:`transform_coords`. + + Returns + ------- + : + New data array with theta coordinate. + """ + graph = graph if graph is not None else specular_reflection() + data_array_theta = data_array.transform_coords(['theta'], graph=graph) + if theta_edges is not None: + data_array_theta = data_array_theta.bin({theta_edges.dim: theta_edges}) + try: + from orsopy import fileio + + orso_measurement = data_array_theta.attrs['orso'].value.data_source.measurement + orso_measurement.instrument_settings.incident_angle = fileio.base.ValueRange( + float(data_array_theta.coords['theta'].min().value), + float(data_array_theta.coords['theta'].max().value), + data_array_theta.bins.coords['theta'].min().unit, + ) + import inspect + + # Determine if 'gravity' is in the graph and if to add the gravity correction + if any( + [ + 'gravity' in i.parameters.keys() + for i in map(inspect.signature, graph.values()) + ] + ): + data_array_theta.attrs['orso'].value.reduction.corrections += [ + 'gravity correction' + ] + except ImportError: + orso.not_found_warning() + return data_array_theta + + +def theta_to_q( + data_array: sc.DataArray, q_edges: sc.Variable = None, graph: dict = None +) -> sc.DataArray: + """ + Convert from theta to Q and if necessary bin in Q. + + Parameters + ---------- + data_array: + Data array to convert. + q_edges: + The lower and upper limits for the Q. If :code:`None`, no + binning is performed. + graph: + Graph for :code:`transform_coords`. + + Returns + ------- + : + New data array with theta coordinate. + """ + graph = graph if graph is not None else specular_reflection() + data_array_q = data_array.transform_coords(["Q"], graph=graph) + if q_edges is not None: + data_array_q = data_array_q.bin({q_edges.dim: q_edges}) + return data_array_q + + +def sum_bins(data_array: sc.DataArray): + """ + Sum the event bins and propagate the maximum resolution, where available. + + Parameters + ---------- + data_array: + Data array to be summed. + + Returns + ------- + : + Summed data array. + """ + data_array_summed = data_array.bins.sum() + if 'angular_resolution' in data_array.bins.coords: + data_array_summed.coords['angular_resolution'] = data_array.bins.coords[ + 'angular_resolution' + ].max('detector_number') + return data_array_summed diff --git a/src/essreflectometry/reflectometry/corrections.py b/src/essreflectometry/reflectometry/corrections.py new file mode 100644 index 0000000..a8833c5 --- /dev/null +++ b/src/essreflectometry/reflectometry/corrections.py @@ -0,0 +1,98 @@ +# SPDX-License-Identifier: BSD-3-Clause +# Copyright (c) 2023 Scipp contributors (https://github.com/scipp) +import numpy as np +import scipp as sc + +from ..amor.tools import fwhm_to_std +from . import orso + + +def footprint_correction(data_array: sc.DataArray) -> sc.DataArray: + """ + Perform the footprint correction on the data array that has a :code:`beam_size` and + binned :code:`theta` values. + + Parameters + ---------- + data_array: + Data array to perform footprint correction on. + + Returns + ------- + : + Footprint corrected data array. + """ + size_of_beam_on_sample = beam_on_sample( + data_array.coords['beam_size'], data_array.bins.coords['theta'] + ) + footprint_scale = sc.erf( + fwhm_to_std(data_array.coords['sample_size'] / size_of_beam_on_sample) + ) + data_array_fp_correction = data_array / footprint_scale.squeeze() + try: + data_array_fp_correction.attrs['orso'].value.reduction.corrections += [ + 'footprint correction' + ] + except KeyError: + orso.not_found_warning() + return data_array_fp_correction + + +def normalize_by_counts(data_array: sc.DataArray) -> sc.DataArray: + """ + Normalize the bin-summed data by the total number of counts. + If the data has variances, a check is performed to ensure that the counts in each + bin is much lower than the total counts. If this is not the case, an error is raised + because the normalization would introduce non-negligible correlations which are not + handled Scipp's basic error propagation. See Heybrock et al. (2023). + If the check passes, the input data is simply divided by the total number of counts, + ignoring the variances of the denominator. + + Parameters + ---------- + data_array: + Data array to be normalized. + + Returns + ------- + : + Normalized data array. + """ + # Dividing by ncounts fails because ncounts also has variances, and this introduces + # correlations. According to Heybrock et al. (2023), we can however safely drop the + # variances of ncounts if counts_in_bin / ncounts is small everywhere. + ncounts = sc.values(data_array.sum()) + norm = data_array / ncounts + if (data_array.variances is not None) and (norm.max().value > 0.1): + ind = np.argmax(data_array.values) + raise ValueError( + 'One or more bins contain a number of counts of the same order as the ' + 'total number of counts. It is not safe to drop the variances of the ' + 'denominator when normalizing by the total number of counts in this ' + f'regime. The maximum counts found is {data_array.values[ind]} at ' + f'index {ind}. The total number of counts is {ncounts.value}.' + ) + try: + norm.attrs['orso'].value.reduction.corrections += ['total counts'] + except KeyError: + orso.not_found_warning() + return norm + + +def beam_on_sample(beam_size: sc.Variable, theta: sc.Variable) -> sc.Variable: + """ + Size of the beam on the sample. + + Parameters + ---------- + beam_size: + Full width half maximum of the beam. + theta: + Angular of incidence with the sample. + + Returns + ------- + : + Size of the beam on the sample. + """ + return beam_size / sc.sin(theta) diff --git a/src/essreflectometry/reflectometry/io.py b/src/essreflectometry/reflectometry/io.py new file mode 100644 index 0000000..df2e982 --- /dev/null +++ b/src/essreflectometry/reflectometry/io.py @@ -0,0 +1,41 @@ +# SPDX-License-Identifier: BSD-3-Clause +# Copyright (c) 2023 Scipp contributors (https://github.com/scipp) + +from typing import Optional + +import numpy as np +import scipp as sc + + +def save_ort( + data_array: sc.DataArray, filename: str, dimension: Optional[str] = None +) -> None: + """ + Save a data array with the ORSO .ort file format. + + Parameters + ---------- + data_array: + Scipp-data array to save. + filename: + Filename. + dimension: + String for dimension to perform mean over. + """ + from orsopy import fileio + + if filename[:-4] == '.ort': + raise ValueError("The expected output file ending is .ort.") + if dimension is not None: + data_array = data_array.mean(dimension) + q = data_array.coords['Q'] + if data_array.coords.is_edges('Q'): + q = sc.midpoints(q) + R = data_array.data + sR = sc.stddevs(data_array.data) + sq = data_array.coords['sigma_Q'] + dataset = fileio.orso.OrsoDataset( + data_array.attrs['orso'].value, + np.array([q.values, R.values, sR.values, sq.values]).T, + ) + fileio.orso.save_orso([dataset], filename) diff --git a/src/essreflectometry/reflectometry/orso.py b/src/essreflectometry/reflectometry/orso.py new file mode 100644 index 0000000..a090676 --- /dev/null +++ b/src/essreflectometry/reflectometry/orso.py @@ -0,0 +1,14 @@ +# SPDX-License-Identifier: BSD-3-Clause +# Copyright (c) 2023 Scipp contributors (https://github.com/scipp) +import warnings + + +def not_found_warning(): + """ + A function to raise a orso specific error if necessary. + """ + warnings.warn( + "For metadata to be logged in the data array, " + "it is necessary to install the orsopy package.", + UserWarning, + ) From 0904f26a8f3a523e480bf9014f4bd77b76124469 Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Tue, 17 Oct 2023 13:16:10 +0200 Subject: [PATCH 2/7] deps: add dependencies --- pyproject.toml | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/pyproject.toml b/pyproject.toml index 9ac8576..e0819f8 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -27,6 +27,14 @@ classifiers = [ ] requires-python = ">=3.8" dependencies = [ + "dask", + "graphviz", + "plopp", + "pythreejs", + "orsopy", + "sciline>=23.9.1", + "scipp>=23.8.0", + "scippneutron>=23.9.0", ] dynamic = ["version"] @@ -42,6 +50,7 @@ addopts = "-ra -v" testpaths = "tests" filterwarnings = [ "error", + "ignore::DeprecationWarning", ] [tool.bandit] From f50aa6a9bd37f77f1fafbf5f94b74617002ac55a Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Tue, 17 Oct 2023 13:16:44 +0200 Subject: [PATCH 3/7] test: add tests from ess.reflectometry --- tests/amor/__init__.py | 2 + tests/amor/tools_test.py | 85 ++++++++++++++++++++++ tests/reflectometry/__init__.py | 2 + tests/reflectometry/corrections_test.py | 93 +++++++++++++++++++++++++ 4 files changed, 182 insertions(+) create mode 100644 tests/amor/__init__.py create mode 100644 tests/amor/tools_test.py create mode 100644 tests/reflectometry/__init__.py create mode 100644 tests/reflectometry/corrections_test.py diff --git a/tests/amor/__init__.py b/tests/amor/__init__.py new file mode 100644 index 0000000..df17698 --- /dev/null +++ b/tests/amor/__init__.py @@ -0,0 +1,2 @@ +# SPDX-License-Identifier: BSD-3-Clause +# Copyright (c) 2023 Scipp contributors (https://github.com/scipp) diff --git a/tests/amor/tools_test.py b/tests/amor/tools_test.py new file mode 100644 index 0000000..31eb753 --- /dev/null +++ b/tests/amor/tools_test.py @@ -0,0 +1,85 @@ +# SPDX-License-Identifier: BSD-3-Clause +# Copyright (c) 2023 Scipp contributors (https://github.com/scipp) +# @author Andrew R. McCluskey (arm61) +import pytest +import scipp as sc + +from essreflectometry.amor import tools + + +def test_linlogspace_linear(): + q_lin = tools.linlogspace( + dim='qz', edges=[0.008, 0.08], scale='linear', num=50, unit='1/angstrom' + ) + expected = sc.linspace(dim='qz', start=0.008, stop=0.08, num=50, unit='1/angstrom') + assert sc.allclose(q_lin, expected) + + +def test_linlogspace_linear_list_input(): + q_lin = tools.linlogspace( + dim='qz', edges=[0.008, 0.08], unit='1/angstrom', scale=['linear'], num=[50] + ) + expected = sc.linspace(dim='qz', start=0.008, stop=0.08, num=50, unit='1/angstrom') + assert sc.allclose(q_lin, expected) + + +def test_linlogspace_log(): + q_log = tools.linlogspace( + dim='qz', edges=[0.008, 0.08], unit='1/angstrom', scale='log', num=50 + ) + expected = sc.geomspace(dim='qz', start=0.008, stop=0.08, num=50, unit='1/angstrom') + assert sc.allclose(q_log, expected) + + +def test_linlogspace_linear_log(): + q_linlog = tools.linlogspace( + dim='qz', + edges=[0.008, 0.03, 0.08], + unit='1/angstrom', + scale=['linear', 'log'], + num=[16, 20], + ) + exp_lin = sc.linspace(dim='qz', start=0.008, stop=0.03, num=16, unit='1/angstrom') + exp_log = sc.geomspace(dim='qz', start=0.03, stop=0.08, num=21, unit='1/angstrom') + expected = sc.concat([exp_lin, exp_log['qz', 1:]], 'qz') + assert sc.allclose(q_linlog, expected) + + +def test_linlogspace_log_linear(): + q_loglin = tools.linlogspace( + dim='qz', + edges=[0.008, 0.03, 0.08], + unit='1/angstrom', + scale=['log', 'linear'], + num=[16, 20], + ) + exp_log = sc.geomspace(dim='qz', start=0.008, stop=0.03, num=16, unit='1/angstrom') + exp_lin = sc.linspace(dim='qz', start=0.03, stop=0.08, num=21, unit='1/angstrom') + expected = sc.concat([exp_log, exp_lin['qz', 1:]], 'qz') + assert sc.allclose(q_loglin, expected) + + +def test_linlogspace_linear_log_linear(): + q_linloglin = tools.linlogspace( + dim='qz', + edges=[0.008, 0.03, 0.08, 0.12], + unit='1/angstrom', + scale=['linear', 'log', 'linear'], + num=[16, 20, 10], + ) + exp_lin = sc.linspace(dim='qz', start=0.008, stop=0.03, num=16, unit='1/angstrom') + exp_log = sc.geomspace(dim='qz', start=0.03, stop=0.08, num=21, unit='1/angstrom') + exp_lin2 = sc.linspace(dim='qz', start=0.08, stop=0.12, num=11, unit='1/angstrom') + expected = sc.concat([exp_lin, exp_log['qz', 1:], exp_lin2['qz', 1:]], 'qz') + assert sc.allclose(q_linloglin, expected) + + +def test_linlogspace_bad_input(): + with pytest.raises(ValueError): + _ = tools.linlogspace( + dim='qz', + edges=[0.008, 0.03, 0.08, 0.12], + unit='1/angstrom', + scale=['linear', 'log'], + num=[16, 20], + ) diff --git a/tests/reflectometry/__init__.py b/tests/reflectometry/__init__.py new file mode 100644 index 0000000..df17698 --- /dev/null +++ b/tests/reflectometry/__init__.py @@ -0,0 +1,2 @@ +# SPDX-License-Identifier: BSD-3-Clause +# Copyright (c) 2023 Scipp contributors (https://github.com/scipp) diff --git a/tests/reflectometry/corrections_test.py b/tests/reflectometry/corrections_test.py new file mode 100644 index 0000000..98924e1 --- /dev/null +++ b/tests/reflectometry/corrections_test.py @@ -0,0 +1,93 @@ +# SPDX-License-Identifier: BSD-3-Clause +# Copyright (c) 2023 Scipp contributors (https://github.com/scipp) +import pytest +import scipp as sc +from orsopy import fileio + +from essreflectometry.reflectometry import corrections + + +def test_normalize_by_counts(): + """ + Tests the corrections.normalize_by_counts function without + a orsopy object present. + """ + N = 50 + values = [1.0] * N + data = sc.Variable( + dims=['position'], unit=sc.units.counts, values=values, variances=values + ) + array = sc.DataArray(data=data) + array_normalized = corrections.normalize_by_counts(array) + result = sc.DataArray( + data=sc.Variable( + dims=['position'], + unit=sc.units.dimensionless, + values=[1 / N] * N, + variances=[1 / (N * N)] * N, + ) + ) + assert sc.allclose(array_normalized.data, result.data) + + +def test_normalize_by_counts_fails_when_ncounts_is_too_small(): + N = 10 + values = [1.0] * N + values[3] = 20.0 + data = sc.Variable( + dims=['position'], unit=sc.units.counts, values=values, variances=values + ) + array = sc.DataArray(data=data) + with pytest.raises( + ValueError, + match='One or more bins contain a number of counts of the same ' + 'order as the total number of counts', + ): + _ = corrections.normalize_by_counts(array) + + +def test_normalize_by_counts_orso(): + """ + Tests the corrections.normalize_by_counts function + with a orsopy object present. + """ + N = 50 + values = [1.0] * N + data = sc.Variable( + dims=['position'], unit=sc.units.counts, values=values, variances=values + ) + array = sc.DataArray(data=data, attrs={'orso': sc.scalar(fileio.orso.Orso.empty())}) + array.attrs['orso'].value.reduction.corrections = [] + array_normalized = corrections.normalize_by_counts(array) + result = sc.DataArray( + data=sc.Variable( + dims=['position'], + unit=sc.units.dimensionless, + values=[1 / N] * N, + variances=[1 / (N * N)] * N, + ) + ) + assert sc.allclose(array_normalized.data, result.data) + assert 'total counts' in array.attrs['orso'].value.reduction.corrections + + +def test_beam_on_sample(): + """ + Tests the corrections.beam_on_sample function. + """ + beam_size = sc.scalar(1.0, unit=sc.units.mm) + theta = sc.scalar(0.1, unit=sc.units.rad) + expected_result = sc.scalar(10.01668613, unit=sc.units.mm) + assert sc.allclose(expected_result, corrections.beam_on_sample(beam_size, theta)) + + +def test_beam_on_sample_array(): + """ + Tests the corrections.beam_on_sample function with an array of theta. + """ + beam_size = sc.scalar(1.0, unit=sc.units.mm) + theta = sc.array(dims=['x'], values=[0.1, 0.5], unit=sc.units.rad) + expected_result = sc.array( + dims=['x'], values=[10.01668613, 2.085829643], unit=sc.units.mm + ) + assert sc.allclose(expected_result, corrections.beam_on_sample(beam_size, theta)) From d19fe355e0090b8369bcf84af8b6121deb29b11f Mon Sep 17 00:00:00 2001 From: jokasimr Date: Tue, 17 Oct 2023 10:58:46 +0000 Subject: [PATCH 4/7] Apply automatic formatting --- src/essreflectometry/amor/conversions.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/essreflectometry/amor/conversions.py b/src/essreflectometry/amor/conversions.py index adf16a6..8623faf 100644 --- a/src/essreflectometry/amor/conversions.py +++ b/src/essreflectometry/amor/conversions.py @@ -9,7 +9,7 @@ def incident_beam( *, source_chopper_1: sc.Variable, source_chopper_2: sc.Variable, - sample_position: sc.Variable + sample_position: sc.Variable, ) -> sc.Variable: """ Compute the incident beam vector from the source chopper position vector, From 2e7cfd67abc28d5d2b3f79ba52f98b50524790af Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Tue, 17 Oct 2023 13:24:34 +0200 Subject: [PATCH 5/7] fix: ignore warnings about missing orso data --- pyproject.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/pyproject.toml b/pyproject.toml index e0819f8..729a052 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -51,6 +51,7 @@ testpaths = "tests" filterwarnings = [ "error", "ignore::DeprecationWarning", + "ignore:For metadata to be logged in the data array, it is necessary to install the orsopy package.:UserWarning", ] [tool.bandit] From 85866883240f97a9bd68c55207f7f9684c1c9985 Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Tue, 17 Oct 2023 13:31:22 +0200 Subject: [PATCH 6/7] deps: min python3.10 --- .copier-answers.yml | 2 +- docs/developer/getting-started.md | 4 ++-- tox.ini | 2 +- 3 files changed, 4 insertions(+), 4 deletions(-) diff --git a/.copier-answers.yml b/.copier-answers.yml index d69633c..09abe54 100644 --- a/.copier-answers.yml +++ b/.copier-answers.yml @@ -4,7 +4,7 @@ _src_path: gh:scipp/copier_template description: Reflectometry data reduction for the European Spallation Source github_linux_image: ubuntu-20.04 max_python: '3.11' -min_python: '3.8' +min_python: '3.10' orgname: scipp projectname: essreflectometry year: 2023 diff --git a/docs/developer/getting-started.md b/docs/developer/getting-started.md index 14902df..70b797c 100644 --- a/docs/developer/getting-started.md +++ b/docs/developer/getting-started.md @@ -40,7 +40,7 @@ Alternatively, if you want a different workflow, take a look at ``tox.ini`` or ` Run the tests using ```sh -tox -e py38 +tox -e py310 ``` (or just `tox` if you want to run all environments). @@ -88,4 +88,4 @@ python -m sphinx -v -b doctest -d build/.doctrees docs build/html python -m sphinx -v -b linkcheck -d build/.doctrees docs build/html ``` ```` -````` \ No newline at end of file +````` diff --git a/tox.ini b/tox.ini index 1c39d14..0573b76 100644 --- a/tox.ini +++ b/tox.ini @@ -1,5 +1,5 @@ [tox] -envlist = py38 +envlist = py310 isolated_build = true [testenv] From 0b7c12bdcba3b84e8771b13657edcc2e0f4662ed Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Tue, 17 Oct 2023 14:26:28 +0200 Subject: [PATCH 7/7] test: more specific deprecation warning ignore --- pyproject.toml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index 729a052..27c49bf 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -50,8 +50,8 @@ addopts = "-ra -v" testpaths = "tests" filterwarnings = [ "error", - "ignore::DeprecationWarning", - "ignore:For metadata to be logged in the data array, it is necessary to install the orsopy package.:UserWarning", + 'ignore:\n.*Sentinel is not a public part of the traitlets API.*:DeprecationWarning', + "ignore:.*metadata to be logged in the data array, it is necessary to install the orsopy package.:UserWarning", ] [tool.bandit]