Rewrite the Amor workflow using Sciline

docs/examples/amor.ipynb

@@ -0,0 +1,144 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Divergent data reduction for Amor\n",
+ "\n",
+ "In this notebook, we will look at how to use the `essreflectometry` package with Sciline, for reflectometry data collected from the PSI instrument [Amor](https://www.psi.ch/en/sinq/amor) in [divergent beam mode](https://www.psi.ch/en/sinq/amor/selene).\n",
+ "\n",
+ "We will begin by importing the modules that are necessary for this notebook."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import scipp as sc\n",
+ "import sciline\n",
+ "from essreflectometry.amor import providers, default_parameters\n",
+ "from essreflectometry.reflectometry.types import *\n",
+ "from essreflectometry.amor.types import SampleRotation"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "params={\n",
+ " **default_parameters,\n",
+ " QBins: sc.geomspace(dim='Q', start=0.008, stop=0.075, num=200, unit='1/angstrom'),\n",
+ " SampleRotation[Sample]: sc.scalar(0.7989, unit='deg'),\n",
+ " Filename[Sample]: \"sample.nxs\",\n",
+ " SampleRotation[Reference]: sc.scalar(0.8389, unit='deg'),\n",
+ " Filename[Reference]: \"reference.nxs\",\n",
+ "}"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "pipeline = sciline.Pipeline(\n",
+ " providers,\n",
+ " params=params\n",
+ ")"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "pipeline.visualize((NormalizedIOfQ, QStd), graph_attr={'rankdir': 'LR'})"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Compute I over Q and the standard deviation of Q\n",
+ "ioq, qstd = pipeline.compute((NormalizedIOfQ, QStd)).values()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import matplotlib.pyplot as plt\n",
+ "\n",
+ "fig = plt.figure(figsize=(5, 7))\n",
+ "ax1 = fig.add_axes([0, 0.55, 1.0, 0.45])\n",
+ "ax2 = fig.add_axes([0, 0.0, 1.0, 0.45])\n",
+ "cax = fig.add_axes([1.05, 0.55, 0.03, 0.45])\n",
+ "fig1 = ioq.plot(norm='log', ax=ax1, cax=cax, grid=True)\n",
+ "fig2 = ioq.mean('detector_number').plot(norm='log', ax=ax2, grid=True)\n",
+ "fig1.canvas.xrange = fig2.canvas.xrange"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Make a $(\\lambda, \\theta)$ map\n",
+ "A good sanity check is to create a two-dimensional map of the counts in $\\lambda$ and $\\theta$ bins. To achieve this, we request the `ThetaData` from the pipeline. In the graph above we can see that `WavelengthData` is required to compute `ThetaData`, therefore it is also present in `ThetaData` so we don't need to require it separately."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from essreflectometry.reflectometry.types import ThetaData\n",
+ "pipeline.compute(ThetaData[Sample])\\\n",
+ " .bins.concat('detector_number')\\\n",
+ " .hist(\n",
+ " theta=sc.linspace(dim='theta', start=0.0, stop=1.2, num=165, unit='deg').to(unit='rad'),\n",
+ " wavelength=sc.linspace(dim='wavelength', start=0, stop=15.0, num=165, unit='angstrom'),\n",
+ " )\\\n",
+ " .plot()\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "This plot can be used to check if the value of the sample rotation angle $\\omega$ is correct. The bright triangles should be pointing back to the origin $\\lambda = \\theta = 0$."
+ ]
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "essreflectometry",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.10.12"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}

src/essreflectometry/amor/__init__.py

@@ -1,7 +1,41 @@
 # 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 itertools import chain
+
+import scipp as sc
+from scipp.constants import g
+
+from ..reflectometry import providers as reflectometry_providers
+from ..reflectometry.types import Run
+from . import beamline, calibrations, conversions, load, normalize, resolution, tools
+
+# from .beamline import instrument_view_components
 from .instrument_view import instrument_view
-from .load import load
+from .types import *
+
+providers = list(
+ chain(
+ reflectometry_providers,
+ load.providers,
+ calibrations.providers,
+ conversions.providers,
+ normalize.providers,
+ resolution.providers,
+ beamline.providers,
+ )
+)
+
+default_parameters = {
+ Supermirror[MValue]: sc.scalar(5, unit=sc.units.dimensionless),
+ Supermirror[CriticalEdge]: 0.022 * sc.Unit('1/angstrom'),
+ Supermirror[Alpha]: sc.scalar(0.25 / 0.088, unit=sc.units.angstrom),
+ BeamSize[Run]: 2.0 * sc.units.mm,
+ SampleSize[Run]: 10.0 * sc.units.mm,
+ DetectorSpatialResolution[Run]: 0.0025 * sc.units.m,
+ Gravity: sc.vector(value=[0, -1, 0]) * g,
+ ChopperFrequency[Run]: sc.scalar(20 / 3, unit='Hz'),
+ ChopperPhase[Run]: sc.scalar(-8.0, unit='deg'),
+ Chopper1Position[Run]: sc.vector(value=[0, 0, -15.5], unit='m'),
+ Chopper2Position[Run]: sc.vector(value=[0, 0, -14.5], unit='m'),
+}

src/essreflectometry/amor/beamline.py

@@ -1,26 +1,37 @@
 # 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
+from ..reflectometry.types import BeamlineParams, Run
+from .types import (
+ BeamSize,
+ Chopper1Position,
+ Chopper2Position,
+ ChopperFrequency,
+ ChopperPhase,
+ DetectorSpatialResolution,
+ Gravity,
+ SampleRotation,
+ SampleSize,
+)
 
 
 @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:
+ sample_rotation: SampleRotation[Run],
+ beam_size: BeamSize[Run],
+ sample_size: SampleSize[Run],
+ detector_spatial_resolution: DetectorSpatialResolution[Run],
+ gravity: Gravity,
+ chopper_frequency: ChopperFrequency[Run],
+ chopper_phase: ChopperPhase[Run],
+ chopper_1_position: Chopper1Position[Run],
+ chopper_2_position: Chopper2Position[Run],
+) -> BeamlineParams[Run]:
  """
  Amor beamline components.
 
@@ -50,22 +61,6 @@ def make_beamline(
  :
  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,
@@ -91,7 +86,7 @@ def make_beamline(
  position=chopper_1_position,
  )
  )
- return beamline
+ return BeamlineParams(beamline)
 
 
 @log_call(instrument='amor', level='DEBUG')
@@ -131,3 +126,6 @@ def instrument_view_components(da: sc.DataArray) -> dict:
  'type': 'disk',
  },
  }
+
+
+providers = [make_beamline]

src/essreflectometry/amor/calibrations.py

@@ -2,15 +2,17 @@
 # Copyright (c) 2023 Scipp contributors (https://github.com/scipp)
 import scipp as sc
 
-from ..reflectometry import orso
+# from ..reflectometry import orso
+from ..reflectometry.types import CalibratedReference, Histogrammed, Reference
+from .types import Alpha, CriticalEdge, MValue, Supermirror
 
 
 def supermirror_calibration(
- data_array: sc.DataArray,
- m_value: sc.Variable = None,
- critical_edge: sc.Variable = None,
- alpha: sc.Variable = None,
-) -> sc.Variable:
+ data_array: Histogrammed[Reference],
+ m_value: Supermirror[MValue],
+ critical_edge: Supermirror[CriticalEdge],
+ alpha: Supermirror[Alpha],
+) -> CalibratedReference:
  """
  Calibrate supermirror measurements
 
@@ -19,39 +21,33 @@ def supermirror_calibration(
  data_array:
  Data array to get q-bins/values from.
  m_value:
- m-value for the supermirror. Defaults to 5.
+ m-value for the supermirror.
  critical_edge:
- Supermirror critical edge. Defaults to 0.022 1/angstrom.
+ Supermirror critical edge.
  alpha:
- Supermirror alpha value. Defaults to 0.25 / 0.088 angstrom.
-
+ Supermirror alpha value.
  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
+ # TODO
+ # try:
+ # data_array_cal.attrs['orso'].value.reduction.corrections += [
+ # 'supermirror calibration'
+ # ]
+ # except KeyError:
+ # orso.not_found_warning()
+ return Histogrammed[CalibratedReference](data_array_cal)
 
 
 def calibration_factor(
  data_array: sc.DataArray,
- m_value: sc.Variable = None,
- critical_edge: sc.Variable = None,
- alpha: sc.Variable = None,
+ m_value: sc.Variable,
+ critical_edge: sc.Variable,
+ alpha: sc.Variable,
 ) -> sc.Variable:
  """
  Return the calibration factor for the supermirror.
@@ -61,23 +57,17 @@ def calibration_factor(
  data_array:
  Data array to get q-bins/values from.
  m_value:
- m-value for the supermirror. Defaults to 5.
+ m-value for the supermirror.
  critical_edge:
- Supermirror critical edge. Defaults to 0.022 1/angstrom.
+ Supermirror critical edge.
  alpha:
- Supermirror alpha value. Defaults to 0.25 / 0.088 angstrom.
+ Supermirror alpha value.
 
  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)
@@ -87,3 +77,6 @@ def calibration_factor(
  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
+
+
+providers = [supermirror_calibration]