Merge pull request #4 from scipp/reflectometry-base

Add reflectometry related code from ess

Author: jokasimr

.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

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
 ```
 ````
-`````
+`````

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,8 @@ addopts = "-ra -v"
 testpaths = "tests"
 filterwarnings = [
   "error",
+  '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]

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

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',
+        },
+    }

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

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

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)

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
+    )