Merge pull request #160 from scipp/sample-rotation-offset

Sample rotation offset

Author: jokasimr

docs/user-guide/amor/amor-reduction.ipynb

@@ -100,7 +100,7 @@
    "source": [
     "workflow[Filename[ReferenceRun]] = amor.data.amor_run(614)\n",
     "# The sample rotation value in the file is slightly off, so we set it manually\n",
-    "workflow[SampleRotation[ReferenceRun]] = sc.scalar(0.65, unit='deg')\n",
+    "workflow[SampleRotationOffset[ReferenceRun]] = sc.scalar(0.05, unit='deg')\n",
     "\n",
     "reference_result = workflow.compute(ReducedReference)\n",
     "# Set the result back onto the pipeline to cache it\n",
@@ -167,29 +167,30 @@
     "    '608': {\n",
     "        # The sample rotation values in the files are slightly off, so we replace\n",
     "        # them with corrected values.\n",
-    "        SampleRotation[SampleRun]: sc.scalar(0.85, unit='deg'),\n",
+    "        SampleRotationOffset[SampleRun]: sc.scalar(0.05, unit='deg'),\n",
     "        Filename[SampleRun]: amor.data.amor_run(608),\n",
     "    },\n",
     "    '609': {\n",
-    "        SampleRotation[SampleRun]: sc.scalar(2.25, unit='deg'),\n",
+    "        SampleRotationOffset[SampleRun]: sc.scalar(0.05, unit='deg'),\n",
     "        Filename[SampleRun]: amor.data.amor_run(609),\n",
     "    },\n",
     "    '610': {\n",
-    "        SampleRotation[SampleRun]: sc.scalar(3.65, unit='deg'),\n",
+    "        SampleRotationOffset[SampleRun]: sc.scalar(0.05, unit='deg'),\n",
     "        Filename[SampleRun]: amor.data.amor_run(610),\n",
     "    },\n",
     "    '611': {\n",
-    "        SampleRotation[SampleRun]: sc.scalar(5.05, unit='deg'),\n",
+    "        SampleRotationOffset[SampleRun]: sc.scalar(0.05, unit='deg'),\n",
     "        Filename[SampleRun]: amor.data.amor_run(611),\n",
     "    },\n",
     "}\n",
     "\n",
     "\n",
     "reflectivity = {}\n",
     "for run_number, params in runs.items():\n",
-    "    workflow[Filename[SampleRun]] = params[Filename[SampleRun]]\n",
-    "    workflow[SampleRotation[SampleRun]] = params[SampleRotation[SampleRun]]\n",
-    "    reflectivity[run_number] = workflow.compute(ReflectivityOverQ).hist()\n",
+    "    wf = workflow.copy()\n",
+    "    for key, value in params.items():\n",
+    "        wf[key] = value\n",
+    "    reflectivity[run_number] = wf.compute(ReflectivityOverQ).hist()\n",
     "\n",
     "sc.plot(reflectivity, norm='log', vmin=1e-4)"
    ]
@@ -267,9 +268,10 @@
     "# Start by computing the `ReflectivityData` for each of the files\n",
     "diagnostics = {}\n",
     "for run_number, params in runs.items():\n",
-    "    workflow[Filename[SampleRun]] = params[Filename[SampleRun]]\n",
-    "    workflow[SampleRotation[SampleRun]] = params[SampleRotation[SampleRun]]\n",
-    "    diagnostics[run_number] = workflow.compute((ReflectivityOverZW, ThetaBins[SampleRun]))\n",
+    "    wf = workflow.copy()\n",
+    "    for key, value in params.items():\n",
+    "        wf[key] = value\n",
+    "    diagnostics[run_number] = wf.compute((ReflectivityOverZW, ThetaBins[SampleRun]))\n",
     "\n",
     "# Scale the results using the scale factors computed earlier\n",
     "for run_number, scale_factor in zip(reflectivity.keys(), scale_factors, strict=True):\n",
@@ -358,7 +360,7 @@
     "from ess.reflectometry.figures import wavelength_z_figure\n",
     "\n",
     "workflow[Filename[SampleRun]] = runs['608'][Filename[SampleRun]]\n",
-    "workflow[SampleRotation[SampleRun]] = runs['608'][SampleRotation[SampleRun]]\n",
+    "workflow[SampleRotationOffset[SampleRun]] = runs['608'][SampleRotationOffset[SampleRun]]\n",
     "wavelength_z_figure(\n",
     "    workflow.compute(Sample),\n",
     "    wavelength_bins=workflow.compute(WavelengthBins),\n",

docs/user-guide/amor/compare-to-eos.ipynb

@@ -96,7 +96,7 @@
     "# Chopper phase value in the file is wrong, so we set it manually\n",
     "workflow[ChopperPhase[ReferenceRun]] = sc.scalar(-7.5, unit='deg')\n",
     "# The sample rotation value in the file is slightly off, so we set it manually\n",
-    "workflow[SampleRotation[ReferenceRun]] = sc.scalar(0.65, unit='deg')\n",
+    "workflow[SampleRotationOffset[ReferenceRun]] = sc.scalar(0.05, unit='deg')\n",
     "workflow[Filename[ReferenceRun]] = amor.data.amor_run(614)\n",
     "\n",
     "reference_result = workflow.compute(ReducedReference)\n",
@@ -123,9 +123,9 @@
     "results = sc.DataGroup({'ess': sc.DataGroup(), 'psi': sc.DataGroup()})\n",
     "\n",
     "# ESS results\n",
-    "for key, angle in mu.items():\n",
+    "for key in mu:\n",
     "    print(key, '... ', end='')\n",
-    "    workflow[SampleRotation[SampleRun]] = sc.scalar(angle + 0.05, unit='deg')\n",
+    "    workflow[SampleRotationOffset[SampleRun]] = sc.scalar(0.05, unit='deg')\n",
     "    workflow[Filename[SampleRun]] = amor.data.amor_run(key)\n",
     "    da = workflow.compute(ReflectivityOverQ).hist()\n",
     "    da.coords['Q'] = sc.midpoints(da.coords['Q'], dim='Q')\n",

src/ess/amor/__init__.py

@@ -14,6 +14,7 @@
     NeXusDetectorName,
     RunType,
     SamplePosition,
+    SampleRotationOffset,
 )
 from . import (
     conversions,
@@ -71,6 +72,7 @@ def default_parameters() -> dict:
             sc.scalar(0.75, unit='deg'),
         ),
         GravityToggle: True,
+        SampleRotationOffset[RunType]: sc.scalar(0.0, unit='deg'),
     }
 
 

src/ess/amor/load.py

@@ -11,6 +11,7 @@
     NeXusDetectorName,
     ProtonCurrent,
     RawDetectorData,
+    RawSampleRotation,
     RunType,
     SampleRotation,
     SampleSize,
@@ -100,7 +101,7 @@ def load_amor_ch_frequency(ch: RawChopper[RunType]) -> ChopperFrequency[RunType]
     raise ValueError("No unit was found for the chopper frequency")
 
 
-def load_amor_sample_rotation(fp: Filename[RunType]) -> SampleRotation[RunType]:
+def load_amor_sample_rotation(fp: Filename[RunType]) -> RawSampleRotation[RunType]:
     (mu,) = load_nx(fp, "NXentry/NXinstrument/master_parameters/mu")
     # Jochens Amor code reads the first value of this log
     # see https://github.com/jochenstahn/amor/blob/140e3192ddb7e7f28acee87e2acaee65ce1332aa/libeos/file_reader.py#L272  # noqa: E501

src/ess/estia/__init__.py

@@ -14,6 +14,7 @@
     NeXusDetectorName,
     RunType,
     SamplePosition,
+    SampleRotationOffset,
 )
 from . import conversions, load, maskings, normalization, orso, resolution, workflow
 from .types import (
@@ -57,6 +58,7 @@ def default_parameters() -> dict:
             sc.scalar(-0.75, unit='deg'),
             sc.scalar(0.75, unit='deg'),
         ),
+        SampleRotationOffset[RunType]: sc.scalar(0.0, unit='deg'),
     }
 
 

src/ess/estia/load.py

@@ -6,12 +6,14 @@
     Filename,
     RawDetectorData,
     RunType,
+    SampleRotationOffset,
 )
 from .mcstas import parse_events_ascii, parse_events_h5
 
 
 def load_mcstas_events(
     filename: Filename[RunType],
+    sample_rotation_offset: SampleRotationOffset[RunType],
 ) -> RawDetectorData[RunType]:
     """
     Load event data from a McStas run and reshape it
@@ -31,6 +33,9 @@ def load_mcstas_events(
     da.coords['detector_rotation'] = 2 * da.coords['sample_rotation'] + sc.scalar(
         1.65, unit='deg'
     )
+    da.coords['sample_rotation'] += sample_rotation_offset.to(
+        unit=da.coords['sample_rotation'].unit
+    )
 
     xbins = sc.linspace('x', -0.25, 0.25, 14 * 32 + 1)
     ybins = sc.linspace('y', -0.25, 0.25, 65)

src/ess/reflectometry/__init__.py

@@ -10,10 +10,11 @@
     __version__ = "0.0.0"
 
 
-from . import conversions, figures, normalization, orso
+from . import conversions, corrections, figures, normalization, orso
 from .load import load_reference, save_reference
 
 providers = (
+    *corrections.providers,
     *conversions.providers,
     *orso.providers,
     *normalization.providers,

src/ess/reflectometry/corrections.py

@@ -2,6 +2,7 @@
 import scipp as sc
 
 from .tools import fwhm_to_std
+from .types import RawSampleRotation, RunType, SampleRotation, SampleRotationOffset
 
 
 def footprint_on_sample(
@@ -45,3 +46,12 @@ def correct_by_footprint(da: sc.DataArray) -> sc.DataArray:
 def correct_by_proton_current(da: sc.DataArray) -> sc.DataArray:
     "Corrects the data by the proton current during the time of data collection"
     return da / da.bins.coords['proton_current']
+
+
+def correct_sample_rotation(
+    mu: RawSampleRotation[RunType], mu_offset: SampleRotationOffset[RunType]
+) -> SampleRotation[RunType]:
+    return mu + mu_offset.to(unit=mu.unit)
+
+
+providers = (correct_sample_rotation,)

src/ess/reflectometry/types.py

@@ -79,10 +79,19 @@ class Filename(sciline.Scope[RunType, str], str):
     """Filename of an event data nexus file."""
 
 
+class RawSampleRotation(sciline.Scope[RunType, sc.Variable], sc.Variable):
+    """The rotation of the sample registered in the NeXus file."""
+
+
 class SampleRotation(sciline.Scope[RunType, sc.Variable], sc.Variable):
     """The rotation of the sample relative to the center of the incoming beam."""
 
 
+class SampleRotationOffset(sciline.Scope[RunType, sc.Variable], sc.Variable):
+    """The difference between the true slope of the sample surface
+    and the sample rotation value in the file."""
+
+
 class DetectorRotation(sciline.Scope[RunType, sc.Variable], sc.Variable):
     """The rotation of the detector relative to the horizon"""
 

tests/amor/pipeline_test.py

@@ -17,10 +17,12 @@
     Filename,
     ProtonCurrent,
     QBins,
+    RawSampleRotation,
     ReducibleData,
     ReferenceRun,
     ReflectivityOverQ,
     SampleRotation,
+    SampleRotationOffset,
     SampleRun,
     SampleSize,
     WavelengthBins,
@@ -226,3 +228,18 @@ def test_proton_current(amor_pipeline: sciline.Pipeline):
     np.testing.assert_allclose(
         proton_current[np.searchsorted(timestamps, t) - 1], w_without / w_with
     )
+
+
+@pytest.mark.filterwarnings("ignore:Failed to convert .* into a transformation")
+@pytest.mark.filterwarnings("ignore:Invalid transformation, missing attribute")
+def test_sample_rotation_offset(amor_pipeline: sciline.Pipeline):
+    amor_pipeline[Filename[SampleRun]] = amor.data.amor_run(608)
+    amor_pipeline[SampleRotationOffset[SampleRun]] = sc.scalar(1.0, unit='deg')
+    mu, muoffset, muraw = amor_pipeline.compute(
+        (
+            SampleRotation[SampleRun],
+            SampleRotationOffset[SampleRun],
+            RawSampleRotation[SampleRun],
+        )
+    ).values()
+    assert mu == muoffset.to(unit=muraw.unit) + muraw