Merge pull request #137 from scipp/estia-orso

feat: add orso providers for estia

Author: jokasimr

src/ess/amor/orso.py

@@ -1,112 +1,16 @@
 # SPDX-License-Identifier: BSD-3-Clause
 # Copyright (c) 2024 Scipp contributors (https://github.com/scipp)
-"""ORSO utilities for Amor."""
+from ess.reflectometry.orso import OrsoCorrectionList
 
-import numpy as np
-import scipp as sc
-from orsopy.fileio import base as orso_base
-from orsopy.fileio import data_source as orso_data_source
-from orsopy.fileio.orso import Column, Orso, OrsoDataset
 
-from ..reflectometry.orso import (
-    OrsoDataSource,
-    OrsoInstrument,
-    OrsoIofQDataset,
-    OrsoReduction,
-)
-from ..reflectometry.types import ReflectivityOverQ
-
-
-def build_orso_instrument(events: ReflectivityOverQ) -> OrsoInstrument:
-    """Build ORSO instrument metadata from intermediate reduction results for Amor.
-
-    This assumes specular reflection and sets the incident angle equal to the computed
-    scattering angle.
-    """
-    return OrsoInstrument(
-        orso_data_source.InstrumentSettings(
-            wavelength=orso_base.ValueRange(*_limits_of_coord(events, "wavelength")),
-            incident_angle=orso_base.ValueRange(*_limits_of_coord(events, "theta")),
-            polarization=None,  # TODO how can we determine this from the inputs?
-        )
-    )
-
-
-def build_orso_iofq_dataset(
-    iofq: ReflectivityOverQ,
-    data_source: OrsoDataSource,
-    reduction: OrsoReduction,
-) -> OrsoIofQDataset:
-    """Build an ORSO dataset for reduced I-of-Q data and associated metadata."""
-    header = Orso(
-        data_source=data_source,
-        reduction=reduction,
-        columns=[
-            Column("Qz", "1/angstrom", "wavevector transfer"),
-            Column("R", None, "reflectivity"),
-            Column("sR", None, "standard deviation of reflectivity"),
-            Column(
-                "sQz",
-                "1/angstrom",
-                "standard deviation of wavevector transfer resolution",
-            ),
-        ],
+def orso_amor_corrections() -> OrsoCorrectionList:
+    return OrsoCorrectionList(
+        [
+            "chopper ToF correction",
+            "footprint correction",
+            "supermirror calibration",
+        ]
     )
-    iofq = iofq.hist()
-
-    qz = iofq.coords["Q"].to(unit="1/angstrom", copy=False)
-    if iofq.coords.is_edges("Q"):
-        qz = sc.midpoints(qz)
-    r = sc.values(iofq.data)
-    sr = sc.stddevs(iofq.data)
-    sqz = iofq.coords["Q_resolution"].to(unit="1/angstrom", copy=False)
-
-    data = np.column_stack(tuple(map(_extract_values_array, (qz, r, sr, sqz))))
-    data = data[np.isfinite(data).all(axis=-1)]
-    ds = OrsoIofQDataset(OrsoDataset(header, data))
-    ds.info.reduction.corrections = [
-        "chopper ToF correction",
-        "footprint correction",
-        "supermirror calibration",
-    ]
-    return ds
-
-
-def _extract_values_array(var: sc.Variable) -> np.ndarray:
-    if var.variances is not None:
-        raise sc.VariancesError(
-            "ORT columns must not have variances. "
-            "Store the uncertainties as standard deviations in a separate column."
-        )
-    if var.ndim != 1:
-        raise sc.DimensionError(f"ORT columns must be one-dimensional, got {var.sizes}")
-    return var.values
-
-
-def _limits_of_coord(data: sc.DataArray, name: str) -> tuple[float, float, str] | None:
-    if (coord := _get_coord(data, name)) is None:
-        return None
-    min_ = coord.min().value
-    max_ = coord.max().value
-    # Explicit conversions to float because orsopy does not like np.float* types.
-    return float(min_), float(max_), _ascii_unit(coord.unit)
-
-
-def _get_coord(data: sc.DataArray, name: str) -> sc.Variable | None:
-    if name in data.coords:
-        return sc.DataArray(data=data.coords[name], masks=data.masks)
-    if (data.bins is not None) and (name in data.bins.coords):
-        # Note that .bins.concat() applies the top-level masks
-        events = data.bins.concat().value
-        return sc.DataArray(data=events.coords[name], masks=events.masks)
-    return None
-
-
-def _ascii_unit(unit: sc.Unit) -> str:
-    unit = str(unit)
-    if unit == "Å":
-        return "angstrom"
-    return unit
 
 
-providers = (build_orso_instrument, build_orso_iofq_dataset)
+providers = (orso_amor_corrections,)

src/ess/estia/__init__.py

@@ -15,7 +15,7 @@
     RunType,
     SamplePosition,
 )
-from . import conversions, load, maskings, normalization, resolution, workflow
+from . import conversions, load, maskings, normalization, orso, resolution, workflow
 from .types import (
     AngularResolution,
     SampleSizeResolution,
@@ -35,6 +35,7 @@
     *maskings.providers,
     *workflow.providers,
     *normalization.providers,
+    *orso.providers,
 )
 """
 List of providers for setting up a Sciline pipeline.

src/ess/estia/orso.py

@@ -0,0 +1,16 @@
+# SPDX-License-Identifier: BSD-3-Clause
+# Copyright (c) 2025 Scipp contributors (https://github.com/scipp)
+from ess.reflectometry.orso import OrsoCorrectionList
+
+
+def orso_estia_corrections() -> OrsoCorrectionList:
+    return OrsoCorrectionList(
+        [
+            "chopper ToF correction",
+            "footprint correction",
+            "supermirror calibration",
+        ]
+    )
+
+
+providers = (orso_estia_corrections,)

src/ess/reflectometry/orso.py

@@ -11,14 +11,19 @@
 from datetime import datetime, timezone
 from typing import NewType
 
+import numpy as np
+import scipp as sc
 from dateutil.parser import parse as parse_datetime
 from orsopy.fileio import base as orso_base
 from orsopy.fileio import data_source, orso, reduction
+from orsopy.fileio.orso import Column, Orso, OrsoDataset
 
 from .load import load_nx
 from .types import (
     Filename,
+    ReducibleData,
     ReferenceRun,
+    ReflectivityOverQ,
     SampleRun,
 )
 
@@ -52,6 +57,8 @@
 OrsoSampleFilenames = NewType("OrsoSampleFilenames", list[orso_base.File])
 """Collection of filenames used to create the ORSO file"""
 
+OrsoCorrectionList = NewType("OrsoCorrectionList", list[str])
+
 
 def parse_orso_experiment(filename: Filename[SampleRun]) -> OrsoExperiment:
     """Parse ORSO experiment metadata from raw NeXus data."""
@@ -171,6 +178,95 @@ def build_orso_data_source(
     )
 
 
+def build_orso_instrument(events: ReducibleData[SampleRun]) -> OrsoInstrument:
+    """Build ORSO instrument metadata from intermediate reduction results.
+
+    This assumes specular reflection and sets the incident angle equal to the computed
+    scattering angle.
+    """
+    return OrsoInstrument(
+        data_source.InstrumentSettings(
+            wavelength=orso_base.ValueRange(*_limits_of_coord(events, "wavelength")),
+            incident_angle=orso_base.ValueRange(*_limits_of_coord(events, "theta")),
+            polarization=None,  # TODO how can we determine this from the inputs?
+        )
+    )
+
+
+def build_orso_iofq_dataset(
+    iofq: ReflectivityOverQ,
+    data_source: OrsoDataSource,
+    reduction: OrsoReduction,
+    corrections: OrsoCorrectionList,
+) -> OrsoIofQDataset:
+    """Build an ORSO dataset for reduced I-of-Q data and associated metadata."""
+    header = Orso(
+        data_source=data_source,
+        reduction=reduction,
+        columns=[
+            Column("Qz", "1/angstrom", "wavevector transfer"),
+            Column("R", None, "reflectivity"),
+            Column("sR", None, "standard deviation of reflectivity"),
+            Column(
+                "sQz",
+                "1/angstrom",
+                "standard deviation of wavevector transfer resolution",
+            ),
+        ],
+    )
+    iofq = iofq.hist()
+
+    qz = iofq.coords["Q"].to(unit="1/angstrom", copy=False)
+    if iofq.coords.is_edges("Q"):
+        qz = sc.midpoints(qz)
+    r = sc.values(iofq.data)
+    sr = sc.stddevs(iofq.data)
+    sqz = iofq.coords["Q_resolution"].to(unit="1/angstrom", copy=False)
+
+    data = np.column_stack(tuple(map(_extract_values_array, (qz, r, sr, sqz))))
+    data = data[np.isfinite(data).all(axis=-1)]
+    ds = OrsoIofQDataset(OrsoDataset(header, data))
+    ds.info.reduction.corrections = list(corrections)
+    return ds
+
+
+def _extract_values_array(var: sc.Variable) -> np.ndarray:
+    if var.variances is not None:
+        raise sc.VariancesError(
+            "ORT columns must not have variances. "
+            "Store the uncertainties as standard deviations in a separate column."
+        )
+    if var.ndim != 1:
+        raise sc.DimensionError(f"ORT columns must be one-dimensional, got {var.sizes}")
+    return var.values
+
+
+def _limits_of_coord(data: sc.DataArray, name: str) -> tuple[float, float, str] | None:
+    if (coord := _get_coord(data, name)) is None:
+        return None
+    min_ = coord.min().value
+    max_ = coord.max().value
+    # Explicit conversions to float because orsopy does not like np.float* types.
+    return float(min_), float(max_), _ascii_unit(coord.unit)
+
+
+def _get_coord(data: sc.DataArray, name: str) -> sc.Variable | None:
+    if name in data.coords:
+        return sc.DataArray(data=data.coords[name], masks=data.masks)
+    if (data.bins is not None) and (name in data.bins.coords):
+        # Note that .bins.concat() applies the top-level masks
+        events = data.bins.concat().value
+        return sc.DataArray(data=events.coords[name], masks=events.masks)
+    return None
+
+
+def _ascii_unit(unit: sc.Unit) -> str:
+    unit = str(unit)
+    if unit == "Å":
+        return "angstrom"
+    return unit
+
+
 providers = (
     build_orso_data_source,
     build_orso_measurement,
@@ -179,4 +275,6 @@ def build_orso_data_source(
     parse_orso_owner,
     parse_orso_sample,
     orso_data_files,
+    build_orso_instrument,
+    build_orso_iofq_dataset,
 )

tests/estia/mcstas_data_test.py

@@ -1,15 +1,20 @@
 # SPDX-License-Identifier: BSD-3-Clause
 # Copyright (c) 2024 Scipp contributors (https://github.com/scipp)
 # flake8: noqa: F403, F405
+from datetime import datetime
+from pathlib import Path
+from zoneinfo import ZoneInfo
 
 import numpy as np
 import pytest
 import sciline
 import scipp as sc
+from orsopy import fileio
 
 from ess.estia import EstiaWorkflow
 from ess.estia.data import estia_mcstas_reference_run, estia_mcstas_sample_run
 from ess.estia.load import load_mcstas_events
+from ess.reflectometry import orso
 from ess.reflectometry.types import (
     BeamDivergenceLimits,
     Filename,
@@ -35,6 +40,7 @@ def estia_mcstas_pipeline() -> sciline.Pipeline:
     wf[ZIndexLimits] = sc.scalar(0), sc.scalar(14 * 32)
     wf[BeamDivergenceLimits] = sc.scalar(-1.0, unit='deg'), sc.scalar(1.0, unit='deg')
     wf[WavelengthBins] = sc.geomspace('wavelength', 3.5, 12, 2001, unit='angstrom')
+    wf[QBins] = sc.geomspace('Q', 0.005, 0.1, 200, unit='1/angstrom')
     wf[ProtonCurrent[SampleRun]] = sc.DataArray(
         sc.array(dims=('time',), values=[]),
         coords={'time': sc.array(dims=('time',), values=[], unit='s')},
@@ -43,7 +49,30 @@ def estia_mcstas_pipeline() -> sciline.Pipeline:
         sc.array(dims=('time',), values=[]),
         coords={'time': sc.array(dims=('time',), values=[], unit='s')},
     )
-
+    wf[orso.OrsoCreator] = orso.OrsoCreator(
+        fileio.base.Person(
+            name="Max Mustermann",
+            affiliation="European Spallation Source ERIC",
+            contact="[email protected]",
+        )
+    )
+    wf[orso.OrsoExperiment] = orso.OrsoExperiment(
+        fileio.data_source.Experiment(
+            title='McStas run',
+            instrument='Estia',
+            facility='ESS',
+            start_date=datetime(2025, 3, 20, tzinfo=ZoneInfo("Europe/Stockholm")),
+            probe='neutron',
+        )
+    )
+    wf[orso.OrsoOwner] = orso.OrsoOwner(
+        fileio.base.Person(
+            name='John Doe',
+            contact='[email protected]',
+            affiliation='ESS',
+        )
+    )
+    wf[orso.OrsoSample] = orso.OrsoSample(fileio.data_source.Sample.empty())
     return wf
 
 
@@ -62,7 +91,6 @@ def test_mcstas_compute_reducible_data(estia_mcstas_pipeline: sciline.Pipeline):
 
 def test_can_compute_reflectivity_curve(estia_mcstas_pipeline: sciline.Pipeline):
     estia_mcstas_pipeline[Filename[SampleRun]] = estia_mcstas_sample_run(11)
-    estia_mcstas_pipeline[QBins] = sc.geomspace('Q', 0.005, 0.1, 200, unit='1/angstrom')
     r = estia_mcstas_pipeline.compute(ReflectivityOverQ)
     assert "Q" in r.coords
     assert "Q_resolution" in r.coords
@@ -80,3 +108,29 @@ def test_can_compute_reflectivity_curve(estia_mcstas_pipeline: sciline.Pipeline)
 
     assert max_q > sc.scalar(0.075, unit='1/angstrom')
     assert min_q < sc.scalar(0.007, unit='1/angstrom')
+
+
+def test_orso_pipeline(estia_mcstas_pipeline: sciline.Pipeline):
+    estia_mcstas_pipeline[Filename[SampleRun]] = estia_mcstas_sample_run(11)
+    res = estia_mcstas_pipeline.compute(orso.OrsoIofQDataset)
+    assert res.info.data_source.experiment.instrument == "Estia"
+    assert res.info.reduction.software.name == "ess.reflectometry"
+    assert res.info.reduction.corrections == [
+        "chopper ToF correction",
+        "footprint correction",
+        "supermirror calibration",
+    ]
+    assert res.data.ndim == 2
+    assert res.data.shape[1] == 4
+    assert np.all(res.data[:, 1] >= 0)
+    assert np.isfinite(res.data).all()
+
+
+def test_save_reduced_orso_file(
+    estia_mcstas_pipeline: sciline.Pipeline, output_folder: Path
+):
+    estia_mcstas_pipeline[Filename[SampleRun]] = estia_mcstas_sample_run(11)
+    res = estia_mcstas_pipeline.compute(orso.OrsoIofQDataset)
+    fileio.orso.save_orso(
+        datasets=[res], fname=output_folder / 'estia_reduced_iofq.ort'
+    )