From 59621bf439c149cf062a0247180f57e98b35f085 Mon Sep 17 00:00:00 2001 From: Simon Heybrock Date: Thu, 19 Oct 2023 11:35:05 +0200 Subject: [PATCH] Rewrite the Amor workflow using Sciline --- docs/examples/amor.ipynb | 144 ++++++++++++++++ src/essreflectometry/amor/__init__.py | 40 ++++- src/essreflectometry/amor/beamline.py | 54 +++--- src/essreflectometry/amor/calibrations.py | 63 ++++--- src/essreflectometry/amor/conversions.py | 8 +- src/essreflectometry/amor/load.py | 84 +++++----- src/essreflectometry/amor/normalize.py | 44 ++--- src/essreflectometry/amor/resolution.py | 57 ++++++- src/essreflectometry/amor/types.py | 54 ++++++ .../reflectometry/__init__.py | 8 + .../reflectometry/conversions.py | 155 +++++++++--------- .../reflectometry/corrections.py | 61 +++++-- src/essreflectometry/reflectometry/types.py | 67 ++++++++ 13 files changed, 616 insertions(+), 223 deletions(-) create mode 100644 docs/examples/amor.ipynb create mode 100644 src/essreflectometry/amor/types.py create mode 100644 src/essreflectometry/reflectometry/types.py diff --git a/docs/examples/amor.ipynb b/docs/examples/amor.ipynb new file mode 100644 index 0000000..750a768 --- /dev/null +++ b/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 +} diff --git a/src/essreflectometry/amor/__init__.py b/src/essreflectometry/amor/__init__.py index a74676e..d4122e7 100644 --- a/src/essreflectometry/amor/__init__.py +++ b/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'), +} diff --git a/src/essreflectometry/amor/beamline.py b/src/essreflectometry/amor/beamline.py index 2599c79..133f807 100644 --- a/src/essreflectometry/amor/beamline.py +++ b/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] diff --git a/src/essreflectometry/amor/calibrations.py b/src/essreflectometry/amor/calibrations.py index fd3d149..38f95c4 100644 --- a/src/essreflectometry/amor/calibrations.py +++ b/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] diff --git a/src/essreflectometry/amor/conversions.py b/src/essreflectometry/amor/conversions.py index 8623faf..cf8e0e3 100644 --- a/src/essreflectometry/amor/conversions.py +++ b/src/essreflectometry/amor/conversions.py @@ -3,6 +3,7 @@ import scipp as sc from ..reflectometry.conversions import specular_reflection as spec_relf_graph +from ..reflectometry.types import SpecularReflectionCoordTransformGraph def incident_beam( @@ -22,10 +23,13 @@ def incident_beam( return sample_position - chopper_midpoint -def specular_reflection() -> dict: +def specular_reflection() -> SpecularReflectionCoordTransformGraph: """ Generate a coordinate transformation graph for Amor reflectometry. """ graph = spec_relf_graph() graph['incident_beam'] = incident_beam - return graph + return SpecularReflectionCoordTransformGraph(graph) + + +providers = [specular_reflection] diff --git a/src/essreflectometry/amor/load.py b/src/essreflectometry/amor/load.py index 42f27a5..4116f2a 100644 --- a/src/essreflectometry/amor/load.py +++ b/src/essreflectometry/amor/load.py @@ -1,14 +1,14 @@ # 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 +from typing import Union import scipp as sc import scippnexus as snx from ..logging import get_logger -from .beamline import make_beamline +from ..reflectometry.types import BeamlineParams, Filename, Raw, Run +from .data import get_path def _tof_correction(data: sc.DataArray, dim: str = 'tof') -> sc.DataArray: @@ -29,8 +29,9 @@ def _tof_correction(data: sc.DataArray, dim: str = 'tof') -> sc.DataArray: : ToF corrected data array. """ - if 'orso' in data.attrs: - data.attrs['orso'].value.reduction.corrections += ['chopper ToF correction'] + # TODO + # 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), @@ -83,19 +84,13 @@ def _load_nexus_entry(filename: Union[str, Path]) -> sc.DataGroup: return f['entry'][()] -def load( - filename: Union[str, Path], - orso: Optional[Any] = None, - beamline: Optional[dict] = None, -) -> sc.DataArray: +def load(filename: Filename[Run], beamline: BeamlineParams[Run]) -> Raw[Run]: """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. @@ -104,6 +99,7 @@ def load( : Data array object for Amor dataset. """ + filename = get_path(filename) get_logger('amor').info( "Loading '%s' as an Amor NeXus file", filename.filename if hasattr(filename, 'filename') else filename, @@ -124,37 +120,47 @@ def load( ) # 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) + # 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) + data = _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', + # Ad-hoc correction described in + # https://scipp.github.io/ess/instruments/amor/amor_reduction.html + data.coords['position'].fields.y += data.coords['position'].fields.z * sc.tan( + 2.0 * data.coords['sample_rotation'] - (0.955 * sc.units.deg) ) - orso.data_source.measurement.data_files = [filename] + return data + + +# TODO +# 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] + + +providers = [load] diff --git a/src/essreflectometry/amor/normalize.py b/src/essreflectometry/amor/normalize.py index ceead34..81b2330 100644 --- a/src/essreflectometry/amor/normalize.py +++ b/src/essreflectometry/amor/normalize.py @@ -2,12 +2,14 @@ # 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 Normalized, NormalizedIOfQ, Reference, Sample def normalize_by_supermirror( - sample: sc.DataArray, supermirror: sc.DataArray -) -> sc.DataArray: + sample: Normalized[Sample], + supermirror: Normalized[Reference], +) -> NormalizedIOfQ: """ Normalize the sample measurement by the (ideally calibrated) supermirror. @@ -26,19 +28,23 @@ def normalize_by_supermirror( """ 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 + # TODO + # 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 NormalizedIOfQ(normalized) + + +providers = [normalize_by_supermirror] diff --git a/src/essreflectometry/amor/resolution.py b/src/essreflectometry/amor/resolution.py index 7591a8e..909dfa4 100644 --- a/src/essreflectometry/amor/resolution.py +++ b/src/essreflectometry/amor/resolution.py @@ -2,10 +2,48 @@ # Copyright (c) 2023 Scipp contributors (https://github.com/scipp) import scipp as sc +from ..reflectometry.types import ( + AngularResolution, + QBins, + QData, + QStd, + Sample, + SampleSizeResolution, + WavelengthResolution, +) from .tools import fwhm_to_std -def wavelength_resolution( +def wavelength_resolution(da: QData[Sample]) -> WavelengthResolution: + return WavelengthResolution( + _wavelength_resolution( + chopper_1_position=da.coords['source_chopper_1'].value['position'], + chopper_2_position=da.coords['source_chopper_2'].value['position'], + pixel_position=da.coords['position'], + ) + ) + + +def angular_resolution(da: QData[Sample]) -> AngularResolution: + return AngularResolution( + _angular_resolution( + pixel_position=da.coords['position'], + theta=da.bins.coords['theta'], + detector_spatial_resolution=da.coords['detector_spatial_resolution'], + ) + ) + + +def sample_size_resolution(da: QData[Sample]) -> SampleSizeResolution: + return SampleSizeResolution( + _sample_size_resolution( + pixel_position=da.coords['position'], + sample_size=da.coords['sample_size'], + ) + ) + + +def _wavelength_resolution( chopper_1_position: sc.Variable, chopper_2_position: sc.Variable, pixel_position: sc.Variable, @@ -38,7 +76,7 @@ def wavelength_resolution( return fwhm_to_std(distance_between_choppers / chopper_detector_distance) -def sample_size_resolution( +def _sample_size_resolution( pixel_position: sc.Variable, sample_size: sc.Variable ) -> sc.Variable: """ @@ -64,7 +102,7 @@ def sample_size_resolution( ) -def angular_resolution( +def _angular_resolution( pixel_position: sc.Variable, theta: sc.Variable, detector_spatial_resolution: sc.Variable, @@ -104,11 +142,11 @@ def angular_resolution( def sigma_Q( - angular_resolution: sc.Variable, - wavelength_resolution: sc.Variable, - sample_size_resolution: sc.Variable, - q_bins: sc.Variable, -) -> sc.Variable: + angular_resolution: AngularResolution, + wavelength_resolution: WavelengthResolution, + sample_size_resolution: SampleSizeResolution, + q_bins: QBins, +) -> QStd: """ Combine all of the components of the resolution and add Q contribution. @@ -133,3 +171,6 @@ def sigma_Q( + wavelength_resolution**2 + sample_size_resolution**2 ).max('detector_number') * sc.midpoints(q_bins) + + +providers = [sigma_Q, angular_resolution, wavelength_resolution, sample_size_resolution] diff --git a/src/essreflectometry/amor/types.py b/src/essreflectometry/amor/types.py new file mode 100644 index 0000000..df9c6ac --- /dev/null +++ b/src/essreflectometry/amor/types.py @@ -0,0 +1,54 @@ +from typing import NewType, TypeVar + +import sciline +import scipp as sc + +from ..reflectometry.types import Run + +# TODO What do they mean? +# Supermirror parameters +MValue = NewType('MValue', str) +CriticalEdge = NewType('CriticalEdge', str) +Alpha = NewType('Alpha', str) +SupermirrorParameter = TypeVar('SupermirrorParameter', MValue, CriticalEdge, Alpha) + + +class Supermirror(sciline.Scope[SupermirrorParameter, sc.Variable], sc.Variable): + """Supermirror parameter scope.""" + + +class SampleRotation(sciline.Scope[Run, sc.Variable], sc.Variable): + """The rotation of the sample / the reference sample.""" + + +class BeamSize(sciline.Scope[Run, sc.Variable], sc.Variable): + """FWHM of the neutron beam.""" + + +class DetectorSpatialResolution(sciline.Scope[Run, sc.Variable], sc.Variable): + # TODO + """Don't know what this is.""" + + +class SampleSize(sciline.Scope[Run, sc.Variable], sc.Variable): + # TODO is this radius or total length? + """Size of the sample.""" + + +Gravity = NewType('Gravity', sc.Variable) + + +class ChopperFrequency(sciline.Scope[Run, sc.Variable], sc.Variable): + """Frequency of the choppers used in the run.""" + + +class ChopperPhase(sciline.Scope[Run, sc.Variable], sc.Variable): + """Phase of the choppers in the run.""" + + +class Chopper1Position(sciline.Scope[Run, sc.Variable], sc.Variable): + """Position of the first chopper relative the source of the beam.""" + + +class Chopper2Position(sciline.Scope[Run, sc.Variable], sc.Variable): + """Position of the second chopper relative to the source of the beam.""" diff --git a/src/essreflectometry/reflectometry/__init__.py b/src/essreflectometry/reflectometry/__init__.py index fcb6518..e5328ac 100644 --- a/src/essreflectometry/reflectometry/__init__.py +++ b/src/essreflectometry/reflectometry/__init__.py @@ -2,5 +2,13 @@ # Copyright (c) 2023 Scipp contributors (https://github.com/scipp) # flake8: noqa: F401 +from itertools import chain from . import conversions, corrections, io + +providers = list( + chain( + conversions.providers, + corrections.providers, + ) +) diff --git a/src/essreflectometry/reflectometry/conversions.py b/src/essreflectometry/reflectometry/conversions.py index bbabb27..c37cd04 100644 --- a/src/essreflectometry/reflectometry/conversions.py +++ b/src/essreflectometry/reflectometry/conversions.py @@ -5,7 +5,18 @@ from scippneutron._utils import elem_dtype, elem_unit from scippneutron.conversion.graph import beamline, tof -from . import orso +# from . import orso +from .types import ( + FootprintCorrected, + Histogrammed, + QBins, + QData, + Raw, + Run, + SpecularReflectionCoordTransformGraph, + ThetaData, + WavelengthData, +) def theta( @@ -86,7 +97,7 @@ def reflectometry_q(wavelength: sc.Variable, theta: sc.Variable) -> sc.Variable: return c * sc.sin(theta.astype(dtype, copy=False)) / wavelength -def specular_reflection() -> dict: +def specular_reflection() -> SpecularReflectionCoordTransformGraph: """ Generate a coordinate transformation graph for specular reflection reflectometry. @@ -101,12 +112,13 @@ def specular_reflection() -> dict: "theta": theta, "Q": reflectometry_q, } - return graph + return SpecularReflectionCoordTransformGraph(graph) def tof_to_wavelength( - data_array: sc.DataArray, wavelength_edges: sc.Variable = None, graph: dict = None -) -> sc.DataArray: + data_array: Raw[Run], + graph: SpecularReflectionCoordTransformGraph, +) -> WavelengthData[Run]: """ Use :code:`transform_coords` to convert from ToF to wavelength, cutoff high and low limits for wavelength, and add necessary ORSO metadata. @@ -115,9 +127,6 @@ def tof_to_wavelength( ---------- 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`. @@ -126,32 +135,31 @@ def tof_to_wavelength( : 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 + # TODO + # 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 WavelengthData[Run](data_array_wav) def wavelength_to_theta( - data_array: sc.DataArray, theta_edges: sc.Variable = None, graph: dict = None -) -> sc.DataArray: + data_array: WavelengthData[Run], + graph: SpecularReflectionCoordTransformGraph, +) -> ThetaData[Run]: """ Use :code:`transform_coords` to find the theta values for the events and potentially add ORSO metadata. @@ -160,9 +168,6 @@ def wavelength_to_theta( ---------- 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`. @@ -171,39 +176,39 @@ def wavelength_to_theta( : 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 + # TODO + # 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 ThetaData[Run](data_array_theta) def theta_to_q( - data_array: sc.DataArray, q_edges: sc.Variable = None, graph: dict = None -) -> sc.DataArray: + data_array: FootprintCorrected[Run], + q_edges: QBins, + graph: SpecularReflectionCoordTransformGraph, +) -> QData[Run]: """ Convert from theta to Q and if necessary bin in Q. @@ -212,8 +217,7 @@ def theta_to_q( data_array: Data array to convert. q_edges: - The lower and upper limits for the Q. If :code:`None`, no - binning is performed. + The lower and upper limits for the Q. graph: Graph for :code:`transform_coords`. @@ -222,14 +226,12 @@ def theta_to_q( : 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 + data_array_q = data_array_q.bin({q_edges.dim: q_edges}) + return QData[Run](data_array_q) -def sum_bins(data_array: sc.DataArray): +def sum_bins(data_array: QData[Run]) -> Histogrammed[Run]: """ Sum the event bins and propagate the maximum resolution, where available. @@ -243,9 +245,12 @@ def sum_bins(data_array: sc.DataArray): : 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 + return Histogrammed[Run](data_array.bins.sum()) + + +providers = [ + tof_to_wavelength, + wavelength_to_theta, + theta_to_q, + sum_bins, +] diff --git a/src/essreflectometry/reflectometry/corrections.py b/src/essreflectometry/reflectometry/corrections.py index a8833c5..e7f73b4 100644 --- a/src/essreflectometry/reflectometry/corrections.py +++ b/src/essreflectometry/reflectometry/corrections.py @@ -4,10 +4,21 @@ import scipp as sc from ..amor.tools import fwhm_to_std -from . import orso +# from . import orso +from .types import ( + CalibratedReference, + FootprintCorrected, + Histogrammed, + Normalized, + Reference, + Run, + Sample, + ThetaData, +) -def footprint_correction(data_array: sc.DataArray) -> sc.DataArray: + +def footprint_correction(data_array: ThetaData[Run]) -> FootprintCorrected[Run]: """ Perform the footprint correction on the data array that has a :code:`beam_size` and binned :code:`theta` values. @@ -29,16 +40,30 @@ def footprint_correction(data_array: sc.DataArray) -> sc.DataArray: 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 + # try: + # data_array_fp_correction.attrs['orso'].value.reduction.corrections += [ + # 'footprint correction' + # ] + # except KeyError: + # orso.not_found_warning() + return FootprintCorrected[Run](data_array_fp_correction) + + +def normalize_sample_by_counts( + data_array: Histogrammed[Sample], +) -> Normalized[Sample]: + return Normalized[Sample](normalize_by_counts(data_array)) + +def normalize_reference_by_counts( + data_array: CalibratedReference, +) -> Normalized[Reference]: + return Normalized[Reference](normalize_by_counts(data_array)) -def normalize_by_counts(data_array: sc.DataArray) -> sc.DataArray: + +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 @@ -72,10 +97,11 @@ def normalize_by_counts(data_array: sc.DataArray) -> sc.DataArray: 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() + # TODO + # try: + # norm.attrs['orso'].value.reduction.corrections += ['total counts'] + # except KeyError: + # orso.not_found_warning() return norm @@ -96,3 +122,10 @@ def beam_on_sample(beam_size: sc.Variable, theta: sc.Variable) -> sc.Variable: Size of the beam on the sample. """ return beam_size / sc.sin(theta) + + +providers = [ + footprint_correction, + normalize_sample_by_counts, + normalize_reference_by_counts, +] diff --git a/src/essreflectometry/reflectometry/types.py b/src/essreflectometry/reflectometry/types.py new file mode 100644 index 0000000..6cd24e4 --- /dev/null +++ b/src/essreflectometry/reflectometry/types.py @@ -0,0 +1,67 @@ +from typing import NewType, TypeVar + +import sciline +import scipp as sc + +Reference = NewType('Reference', str) +Sample = NewType('Sample', str) +Run = TypeVar('Run', Reference, Sample) + + +class Raw(sciline.Scope[Run, sc.DataArray], sc.DataArray): + """Raw data""" + + +class WavelengthData(sciline.Scope[Run, sc.DataArray], sc.DataArray): + """Raw data transformed to wavelength""" + + +class ThetaData(sciline.Scope[Run, sc.DataArray], sc.DataArray): + """Wavelength data transformed to theta""" + + +class QData(sciline.Scope[Run, sc.DataArray], sc.DataArray): + """Theta data transformed to momentum transfer""" + + +QStd = NewType('QStd', sc.Variable) + + +class FootprintCorrected(sciline.Scope[Run, sc.DataArray], sc.DataArray): + """Experiment data corrected by footprint on sample""" + + +WavelengthResolution = NewType('WavelengthResolution', sc.Variable) +AngularResolution = NewType('AngularResolution', sc.Variable) +SampleSizeResolution = NewType('SampleSizeResolution', sc.Variable) + + +class BeamlineParams(sciline.Scope[Run, dict], dict): + """Parameters describing the beamline""" + + +SpecularReflectionCoordTransformGraph = NewType( + 'SpecularReflectionCoordTransformGraph', dict +) + +CalibratedReference = NewType('CalibratedReference', sc.DataArray) + + +class Histogrammed(sciline.Scope[Run, sc.DataArray], sc.DataArray): + """Histogrammmed by Q and detector_number""" + + +class Normalized(sciline.Scope[Run, sc.DataArray], sc.DataArray): + """Normalized histogram""" + + +NormalizedIOfQ = NewType('NormalizedIOfQ', sc.DataArray) + + +''' Parameters for the workflow ''' + +QBins = NewType('QBins', sc.Variable) + + +class Filename(sciline.Scope[Run, str], str): + """Filename of the raw data"""