From 7497ec30279a0a2a50f98c36500bb5e156e3b990 Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Fri, 11 Sep 2026 16:18:34 +0200 Subject: [PATCH 01/16] feat: freia choppers and more scaffolding --- .../freia/freia-mcstas-visualization.ipynb | 213 ++++++++++ .../freia/freia-wavelength-lookup-table.ipynb | 192 +++++++++ .../docs/user-guide/freia/index.md | 14 +- packages/essreflectometry/pyproject.toml | 1 + .../src/ess/freia/__init__.py | 12 +- .../essreflectometry/src/ess/freia/data.py | 16 +- .../essreflectometry/src/ess/freia/mcstas.py | 381 +++++++++++++++++- .../src/ess/freia/workflow.py | 42 +- .../tests/freia/mcstas_test.py | 226 +++++++++++ .../tests/freia/workflow_test.py | 22 - 10 files changed, 1068 insertions(+), 51 deletions(-) create mode 100644 packages/essreflectometry/docs/user-guide/freia/freia-mcstas-visualization.ipynb create mode 100644 packages/essreflectometry/docs/user-guide/freia/freia-wavelength-lookup-table.ipynb create mode 100644 packages/essreflectometry/tests/freia/mcstas_test.py diff --git a/packages/essreflectometry/docs/user-guide/freia/freia-mcstas-visualization.ipynb b/packages/essreflectometry/docs/user-guide/freia/freia-mcstas-visualization.ipynb new file mode 100644 index 000000000..8727d1cd7 --- /dev/null +++ b/packages/essreflectometry/docs/user-guide/freia/freia-mcstas-visualization.ipynb @@ -0,0 +1,213 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "0", + "metadata": {}, + "source": [ + "# FREIA McStas detector data\n", + "\n", + "This notebook loads the final FREIA detector with `mcstastox`, visualizes weighted events with Scipp and Plopp, and computes wavelengths with ESSreduce's analytical frame-unwrapping workflow. It is an initial detector-data inspection workflow; reflectivity normalization and conversion to Q are later steps.\n", + "\n", + "Use an environment containing the local `essreflectometry` and `essreduce` packages, Jupyter, and `ipympl` (the reflectometry documentation environment includes the plotting dependencies).\n", + "\n", + "The input must contain the **final `Multiblade` Mantid banana detector**, including an event list and pixel geometry. Enable its `mantid banana ... list all neutrons` output when running McStas. A `Multiblade_histogram` output alone has no event times and cannot be unwrapped. The loader reports an error if the detector events are absent.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1", + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib widget\n", + "import plopp as pp\n", + "import scipp as sc\n", + "\n", + "from ess import freia\n", + "from ess.freia import data\n", + "from ess.reduce.nexus.types import DiskChoppers, Filename, RawDetector\n", + "from ess.reduce.unwrap import (\n", + " ChopperFrameSequence,\n", + " DistanceResolution,\n", + " PulsePeriod,\n", + " TimeResolution,\n", + " WavelengthDetector,\n", + ")\n", + "from ess.reduce.unwrap.types import KeepEventTimeOffset\n", + "from ess.reflectometry.types import SampleRun\n" + ] + }, + { + "cell_type": "markdown", + "id": "2", + "metadata": {}, + "source": [ + "## Select a run\n", + "\n", + "Load the sample run using `ess.freia.data`.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3", + "metadata": {}, + "outputs": [], + "source": [ + "freia_mcstas = freia.FreiaMcStasWorkflow(wavelength_from='analytical')\n", + "freia_mcstas[Filename[SampleRun]] = data.freia_mcstas_sample_run()\n", + "freia_mcstas[KeepEventTimeOffset] = True\n", + "freia_mcstas[TimeResolution] = sc.scalar(20.0, unit='us')\n", + "freia_mcstas[DistanceResolution] = sc.scalar(0.01, unit='m')\n" + ] + }, + { + "cell_type": "markdown", + "id": "4", + "metadata": {}, + "source": [ + "## Inspect the WFM chopper cascade\n", + "\n", + "`freia.mcstas.wfm_choppers()` constructs the fixed **WFM** configuration used by this workflow. The configuration contains three bandwidth disks and five WFM disks with seven openings each. It is independent of the input filename.\n", + "\n", + "For another configuration, assign its chopper dictionary to `freia_mcstas[DiskChoppers[SampleRun]]`.\n", + "\n", + "The analytical model approximates the beam by a central ray. It projects chopper positions onto the global z axis and estimates detector flight paths as source-to-sample plus sample-to-pixel distance; it does not trace the curved guide or model the finite beam footprint on each disk.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5", + "metadata": {}, + "outputs": [], + "source": [ + "choppers = freia_mcstas.compute(DiskChoppers[SampleRun])\n", + "sc.DataGroup(choppers)\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6", + "metadata": {}, + "outputs": [], + "source": [ + "frames = freia_mcstas.compute(ChopperFrameSequence[SampleRun])\n", + "frames.draw()\n" + ] + }, + { + "cell_type": "markdown", + "id": "7", + "metadata": {}, + "source": [ + "## Load and unwrap the detector events\n", + "\n", + "The generic workflow combines the McStas input providers with ESSreduce's wavelength calculation. We request both outputs together so that the file's detector events are loaded once. Weights are retained, with per-event variances equal to the squared weights; empty pixels are also retained.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8", + "metadata": {}, + "outputs": [], + "source": [ + "results = freia_mcstas.compute((RawDetector[SampleRun], WavelengthDetector[SampleRun]))\n", + "raw = results[RawDetector[SampleRun]]\n", + "unwrapped = results[WavelengthDetector[SampleRun]]\n", + "raw\n" + ] + }, + { + "cell_type": "markdown", + "id": "9", + "metadata": {}, + "source": [ + "## Detector image and arrival times\n", + "\n", + "`longitude` and `height` are coordinates in the banana detector's local frame; the `position` vectors use the global instrument frame. All image intensities are sums of event weights.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "10", + "metadata": {}, + "outputs": [], + "source": [ + "detector_image = raw.hist(longitude=120, height=80, dim=raw.dims)\n", + "pp.plot(detector_image, norm='log', title='FREIA final detector')\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "11", + "metadata": {}, + "outputs": [], + "source": [ + "arrival_times = raw.hist(\n", + " event_time_offset=sc.linspace('event_time_offset', 0.0, freia_mcstas.compute(PulsePeriod).to(unit='s').value, 501, unit='s'),\n", + " dim=raw.dims,\n", + ")\n", + "arrival_times.coords['event_time_offset'] = arrival_times.coords['event_time_offset'].to(unit='ms')\n", + "pp.plot(arrival_times, title='Arrival time within the source period')\n" + ] + }, + { + "cell_type": "markdown", + "id": "12", + "metadata": {}, + "source": [ + "## Wavelengths from analytical frame unwrapping\n", + "\n", + "The chopper cascade and arrival times determine the wavelengths. Events outside the modeled transmission bands or above the workflow's relative wavelength uncertainty threshold have NaN wavelengths and do not contribute to the wavelength histograms. Inspect the assigned-event count when assessing the result.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "13", + "metadata": {}, + "outputs": [], + "source": [ + "event_wavelengths = unwrapped.bins.constituents['data'].coords['wavelength']\n", + "valid = sc.isfinite(event_wavelengths)\n", + "print(f'{sc.sum(valid).value:,} / {valid.size:,} events have an assigned wavelength')\n", + "\n", + "wavelength_bins = sc.linspace('wavelength', 1.0, 12.0, 441, unit='angstrom')\n", + "spectrum = unwrapped.hist(wavelength=wavelength_bins, dim=unwrapped.dims)\n", + "pp.plot(spectrum, title='FREIA wavelength spectrum')\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "14", + "metadata": {}, + "outputs": [], + "source": [ + "wavelength_image = unwrapped.hist(\n", + " wavelength=wavelength_bins, longitude=120, dim=unwrapped.dims\n", + ")\n", + "pp.plot(wavelength_image, norm='log', title='Wavelength across the detector')\n" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "name": "python" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/packages/essreflectometry/docs/user-guide/freia/freia-wavelength-lookup-table.ipynb b/packages/essreflectometry/docs/user-guide/freia/freia-wavelength-lookup-table.ipynb new file mode 100644 index 000000000..34a85e21b --- /dev/null +++ b/packages/essreflectometry/docs/user-guide/freia/freia-wavelength-lookup-table.ipynb @@ -0,0 +1,192 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "0", + "metadata": {}, + "source": [ + "# FREIA analytical wavelength lookup table\n", + "\n", + "Construct and visualize a wavelength lookup table using FREIA's fixed WFM chopper configuration and ESSreduce's generic analytical frame-unwrapping workflow.\n", + "\n", + "The calculation uses the final detector's pixel geometry and the source and sample positions. Detector events and histogram intensities are not needed to construct the lookup table.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1", + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib widget\n", + "from pathlib import Path\n", + "\n", + "import plopp as pp\n", + "import scipp as sc\n", + "from scippnexus import NXdetector\n", + "\n", + "from ess import freia\n", + "from ess.freia import data\n", + "from ess.reduce.nexus.types import DiskChoppers, Filename\n", + "from ess.reduce.unwrap import (\n", + " ChopperFrameSequence,\n", + " DetectorLtotal,\n", + " DistanceResolution,\n", + " LookupTable,\n", + " SourceBounds,\n", + " TimeResolution,\n", + ")\n", + "from ess.reflectometry.types import SampleRun\n" + ] + }, + { + "cell_type": "markdown", + "id": "2", + "metadata": {}, + "source": [ + "## Use the WFM configuration\n", + "\n", + "The workflow uses `freia.mcstas.wfm_choppers()` to construct the fixed **WFM** cascade. Detector geometry is read from the selected file, and ESSreduce computes the wavelength table analytically.\n", + "\n", + "This example uses 20 µs time resolution and 5 mm flight-path resolution. Source bounds are the generic ESS defaults, shown explicitly below. The central-ray model projects chopper positions onto the global z axis and approximates detector flight paths by source-to-sample plus sample-to-pixel distance.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3", + "metadata": {}, + "outputs": [], + "source": [ + "filename = data.freia_mcstas_sample_run()\n", + "\n", + "freia_mcstas = freia.FreiaMcStasWorkflow(wavelength_from='analytical')\n", + "freia_mcstas[Filename[SampleRun]] = filename\n", + "freia_mcstas[TimeResolution] = sc.scalar(20.0, unit='us')\n", + "freia_mcstas[DistanceResolution] = sc.scalar(0.005, unit='m')\n", + "freia_mcstas[SourceBounds] = SourceBounds(\n", + " time=(sc.scalar(0.0, unit='ms'), sc.scalar(5.0, unit='ms')),\n", + " wavelength=(sc.scalar(0.001, unit='angstrom'), sc.scalar(15.0, unit='angstrom')),\n", + ")\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4", + "metadata": {}, + "outputs": [], + "source": [ + "results = freia_mcstas.compute((\n", + " DiskChoppers[SampleRun],\n", + " ChopperFrameSequence[SampleRun],\n", + " DetectorLtotal[SampleRun],\n", + " LookupTable[SampleRun, NXdetector],\n", + "))\n", + "choppers = results[DiskChoppers[SampleRun]]\n", + "frames = results[ChopperFrameSequence[SampleRun]]\n", + "flight_paths = results[DetectorLtotal[SampleRun]]\n", + "lookup = results[LookupTable[SampleRun, NXdetector]]\n", + "\n", + "print(f'Input: {filename.name}')\n", + "print(f'{len(choppers)} disks; {len(frames.frames[-1].subframes)} transmitted subframes')\n", + "print(f'Detector flight paths: {flight_paths.min().value:.4f} to {flight_paths.max().value:.4f} m')\n", + "print(f'Table shape: {lookup.array.sizes}')\n" + ] + }, + { + "cell_type": "markdown", + "id": "5", + "metadata": {}, + "source": [ + "## Wavelength as a function of arrival time and flight path\n", + "\n", + "The heatmap shows the wavelength assigned to each arrival time and flight path. White regions have no modeled transmission. The lower panel shows a slice near the middle of the detector's flight-path range; shading spans the modeled lower and upper wavelength bounds. The table stores the squared half-width of this range as its variance; it describes wavelength ambiguity, not counting statistics.\n", + "\n", + "This is the unmasked lookup table: no additional relative-uncertainty cutoff has been applied. The lookup grid extends slightly beyond the detector's flight-path range to support interpolation.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6", + "metadata": {}, + "outputs": [], + "source": [ + "table = lookup.array.copy()\n", + "table.coords['event_time_offset'] = table.coords['event_time_offset'].to(unit='ms')\n", + "midpoint = (flight_paths.min() + flight_paths.max()) / 2\n", + "slice_index = int(abs(table.coords['distance'] - midpoint).values.argmin())\n", + "line = table['distance', slice_index].copy()\n", + "\n", + "heatmap = pp.plot(\n", + " sc.values(table),\n", + " title=f'FREIA wavelength lookup — {filename.name}',\n", + " xlabel='Arrival time within the source period [ms]',\n", + " ylabel='Source-to-detector flight path [m]',\n", + " cmap='viridis',\n", + ")\n", + "curve = pp.plot(\n", + " sc.values(line),\n", + " marker='',\n", + " linestyle='-',\n", + " linewidth=1.3,\n", + " title=f\"Slice at {line.coords['distance'].value:.4f} m\",\n", + " xlabel='Arrival time within the source period [ms]',\n", + " ylabel='Wavelength [Å]',\n", + ")\n", + "time = line.coords['event_time_offset'].values\n", + "half_width = sc.stddevs(line.data).values\n", + "\n", + "figure = pp.tiled(2, 1, figsize=(10, 7))\n", + "figure[0, 0] = heatmap\n", + "figure[1, 0] = curve\n", + "figure[0, 0].cax.set_ylabel('Wavelength [Å]')\n", + "figure[0, 0].ax.axhline(line.coords['distance'].value, color='black', linestyle='--', linewidth=0.8)\n", + "figure[1, 0].ax.fill_between(time, line.values - half_width, line.values + half_width, alpha=0.3)\n", + "figure[1, 0].ax.set_xlim(0, lookup.pulse_period.to(unit='ms').value)\n", + "figure[1, 0].ax.grid(alpha=0.2)\n", + "figure\n" + ] + }, + { + "cell_type": "markdown", + "id": "7", + "metadata": {}, + "source": [ + "## Save the lookup table and figure\n", + "\n", + "The HDF5 file retains the table, variances, pulse period, and interpolation resolutions. The PNG and PDF are standalone figures. Outputs are written to `build/freia` relative to the current working directory.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8", + "metadata": {}, + "outputs": [], + "source": [ + "output_dir = Path('build') / 'freia'\n", + "output_dir.mkdir(parents=True, exist_ok=True)\n", + "stem = f'{filename.stem}-wavelength-lookup-table'\n", + "lookup.save_hdf5(output_dir / f'{stem}.h5')\n", + "figure.save(str(output_dir / f'{stem}.png'), dpi=180, bbox_inches='tight')\n", + "figure.save(str(output_dir / f'{stem}.pdf'), bbox_inches='tight')\n", + "output_dir\n" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "name": "python" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/packages/essreflectometry/docs/user-guide/freia/index.md b/packages/essreflectometry/docs/user-guide/freia/index.md index 144f7c67c..e51482d56 100644 --- a/packages/essreflectometry/docs/user-guide/freia/index.md +++ b/packages/essreflectometry/docs/user-guide/freia/index.md @@ -1,3 +1,15 @@ # FREIA -FREIA-specific reduction guides will be added here once example workflows are available. +The initial FREIA workflow loads final-detector McStas events and uses the generic +ESSreduce analytical frame-unwrapping workflow to compute wavelengths. The notebook +guides below use the example download helpers in `ess.freia.data` and visualize the +detector with Scipp and Plopp. Choppers use the fixed **WFM** simulation configuration. +The wavelength lookup-table guide reads only detector geometry from the file, so it +can also be run on simulations without detector events. + +```{toctree} +:maxdepth: 1 + +freia-mcstas-visualization +freia-wavelength-lookup-table +``` diff --git a/packages/essreflectometry/pyproject.toml b/packages/essreflectometry/pyproject.toml index ba7b4b3b2..30880fa3b 100644 --- a/packages/essreflectometry/pyproject.toml +++ b/packages/essreflectometry/pyproject.toml @@ -31,6 +31,7 @@ requires-python = ">=3.12" dependencies = [ "dask>=2022.1.0", "graphviz>=0.20", + "mcstastox>=0.0.11", "python-dateutil>=2.9.0", "plopp>=26.5.0", "orsopy>=1.2", diff --git a/packages/essreflectometry/src/ess/freia/__init__.py b/packages/essreflectometry/src/ess/freia/__init__.py index 3860491e8..7e64e627d 100644 --- a/packages/essreflectometry/src/ess/freia/__init__.py +++ b/packages/essreflectometry/src/ess/freia/__init__.py @@ -3,7 +3,16 @@ import importlib.metadata from ..reflectometry import supermirror -from . import conversions, load, maskings, normalization, orso, resolution, workflow +from . import ( + conversions, + load, + maskings, + mcstas, + normalization, + orso, + resolution, + workflow, +) from .types import ( AngularResolution, SampleSizeResolution, @@ -44,6 +53,7 @@ "conversions", "load", "maskings", + "mcstas", "normalization", "orso", "resolution", diff --git a/packages/essreflectometry/src/ess/freia/data.py b/packages/essreflectometry/src/ess/freia/data.py index 7c92c113e..277709ad4 100644 --- a/packages/essreflectometry/src/ess/freia/data.py +++ b/packages/essreflectometry/src/ess/freia/data.py @@ -1,10 +1,24 @@ +# SPDX-License-Identifier: BSD-3-Clause # Copyright (c) 2025 Scipp contributors (https://github.com/scipp) from ess.reduce.data import make_registry +from ..reflectometry.types import Filename, ReferenceRun, SampleRun + _registry = make_registry( "ess/freia", version="1", files={}, ) -__all__ = [] + +def freia_mcstas_sample_run() -> Filename[SampleRun]: + """Return path to the McStas sample events file.""" + return Filename[SampleRun](_registry.get_path("mcstas-sample.h5")) + + +def freia_mcstas_reference_run() -> Filename[ReferenceRun]: + """Return path to the McStas reference events file.""" + return Filename[ReferenceRun](_registry.get_path("mcstas-reference.h5")) + + +__all__ = ["freia_mcstas_reference_run", "freia_mcstas_sample_run"] diff --git a/packages/essreflectometry/src/ess/freia/mcstas.py b/packages/essreflectometry/src/ess/freia/mcstas.py index 4c48502ca..edb8701c2 100644 --- a/packages/essreflectometry/src/ess/freia/mcstas.py +++ b/packages/essreflectometry/src/ess/freia/mcstas.py @@ -1,3 +1,380 @@ -# Copyright (c) 2025 Scipp contributors (https://github.com/scipp) +# SPDX-License-Identifier: BSD-3-Clause +# Copyright (c) 2026 Scipp contributors (https://github.com/scipp) +"""Adapters for FREIA McStas files. -providers = () +McStas conventions and the fixed WFM simulation configuration belong here. +The rest of the workflow uses the standard ESSreduce domain types. +""" + +import re +from pathlib import Path +from typing import TypedDict + +import mcstastox +import numpy as np +import scipp as sc +import scippnexus as snx +from ess.reduce.nexus.types import ( + DiskChoppers, + EmptyDetector, + Filename, + NeXusDetectorName, + Position, + RawDetector, + RunType, +) +from ess.reduce.unwrap import PulsePeriod +from scippneutron.chopper import DiskChopper + + +class _ChopperParameters(TypedDict): + frequency: float + position: list[float] + open: list[float] + close: list[float] + + +# WFM settings for FREIA_surface_test.instr with res=0.02, stopWBC=stopWFM=0. +# Frequencies are in Hz, global positions in m, and slit angles in degrees. +# Slit angles use ScippNeutron's convention with zero phase and beam position. +_WFM_PARAMETERS: dict[str, _ChopperParameters] = { + 'WBC1': { + 'frequency': 14.0, + 'position': [-0.025536188608227407, -0.2143487326684282, 6.313874089082927], + 'open': [265.75340934340136], + 'close': [347.41993036564133], + }, + 'PSC1': { + 'frequency': 56.0, + 'position': [-0.04279334439590373, -0.2487827747282733, 7.3252399717676475], + 'open': [ + 225.05060997711988, + 174.31325763665376, + 125.39199315800286, + 79.4636849757332, + 35.19478493987566, + 353.0332202961111, + 314.3279223973939, + ], + 'close': [ + 231.97085763665373, + 183.04959315800286, + 136.45713088537676, + 92.08261269250694, + 49.820457639881454, + 369.57011076447253, + 331.2358038278435, + ], + }, + 'PSC2': { + 'frequency': 56.0, + 'position': [-0.05036143156257555, -0.2612121196231304, 7.6903035721702], + 'open': [ + 218.75947574117998, + 165.5769221153046, + 114.32685543062891, + 66.18060316860294, + 19.799339992501075, + 335.62596717152013, + 295.0046293340228, + ], + 'close': [ + 225.67972340071387, + 174.31325763665373, + 125.39199315800283, + 78.79953088537671, + 34.42501269250694, + 352.16285763988157, + 311.9125107644725, + ], + }, + 'PSC3': { + 'frequency': 56.0, + 'position': [-0.06728900511986453, -0.2850923564964583, 8.391692529125612], + 'open': [ + 206.1409753631455, + 148.79069468752678, + 93.06600550591139, + 40.84355449212529, + 350.3878524351604, + 302.33399651829416, + 256.56765150000723, + ], + 'close': [ + 218.89724809831387, + 162.22530499276786, + 108.60576573022513, + 56.985424074290464, + 368.2383919327576, + 321.81186800372, + 277.5953889804347, + ], + }, + 'WBC2': { + 'frequency': 14.0, + 'position': [-0.08633018876666028, -0.3139402382145939, 9.238986706984262], + 'open': [227.07624477736135], + 'close': [339.9669648456413], + }, + 'PSC4': { + 'frequency': 42.0, + 'position': [-0.10516748527471297, -0.35051394569767463, 10.313196928583936], + 'open': [ + 218.21716131667478, + 166.56540575149816, + 115.57883431227881, + 68.09433093440703, + 21.893063791537656, + 337.8911186406503, + 297.67374389476004, + ], + 'close': [ + 238.93286862631388, + 186.29390878964273, + 136.42670903352186, + 88.51839823856325, + 43.20268716658999, + 360.04414290207035, + 317.96091970176843, + ], + }, + 'PSC5': { + 'frequency': 28.0, + 'position': [-0.11514600144160746, -0.5309902203230933, 15.613984633755468], + 'open': [ + 220.38174306267877, + 168.49262474947486, + 118.50890129461459, + 72.09339945600789, + 28.148411442764086, + 344.6532571208708, + 303.481198960639, + ], + 'close': [ + 251.37265421671387, + 197.68679894272015, + 146.514171366115, + 97.48011281665286, + 50.54011861520438, + 366.1021937591697, + 323.6450492619195, + ], + }, + 'WBC3': { + 'frequency': 14.0, + 'position': [-0.10215007498652816, -0.5905364416738501, 17.362922984132187], + 'open': [129.22970188564136], + 'close': [305.77460188564135], + }, +} + + +def wfm_choppers() -> DiskChoppers[RunType]: + """Construct the fixed WFM cascade used by the initial FREIA simulations. + + Returns three bandwidth disks and five pulse-shaping/frame-overlap disks, + each of the latter with seven openings. Settings are independent of the + input file; a fresh set of disks is returned on each call. + + Positions are in the simulation's global frame, with the source at the origin. + The analytical cascade projects them onto the z axis, approximating the + curved guide and finite beam width by a central ray. + """ + return DiskChoppers[RunType]( + { + name: DiskChopper( + frequency=sc.scalar(parameters['frequency'], unit='Hz'), + beam_position=sc.scalar(0.0, unit='deg'), + phase=sc.scalar(0.0, unit='deg'), + axle_position=sc.vector(parameters['position'], unit='m'), + slit_begin=sc.array( + dims=['cutout'], values=parameters['open'], unit='deg' + ), + slit_end=sc.array( + dims=['cutout'], values=parameters['close'], unit='deg' + ), + ) + for name, parameters in _WFM_PARAMETERS.items() + } + ) + + +def _text(value) -> str: + return value.decode() if isinstance(value, bytes) else str(value) + + +def _open_mcstas(filename: str | Path) -> mcstastox.Read: + filename = Path(filename) + return mcstastox.Read(filename.parent, filename.name) + + +def load_mcstas( + filename: str | Path, + detector_name: str = 'Multiblade', + *, + pulse_period: sc.Variable | None = None, +) -> sc.DataArray: + """Load weighted events and pixel geometry from the final FREIA detector. + + Only the selected component is read. The expected output is the + ``mantid banana`` event list, including ``p``, ``t``, ``id`` and pixel + geometry. Histogram-only and upstream debug monitors are not substitutes. + Empty detector pixels are retained, and weighted-event variances are ``p**2``. + Arrival times are split into ``event_time_zero`` and ``event_time_offset`` + using ``pulse_period``, which defaults to the ESS period of 1/14 s. + """ + if pulse_period is None: + pulse_period = sc.scalar(1 / 14, unit='s') + with _open_mcstas(filename) as data: + return _load_events(data, detector_name, pulse_period) + + +def _load_events( + data: mcstastox.Read, + detector_name: str, + pulse_period: sc.Variable, + geometry: sc.DataArray | None = None, +) -> sc.DataArray: + if detector_name not in data.get_components_with_ids(): + raise ValueError( + f'No Mantid detector events with pixel IDs found for {detector_name!r}. ' + 'Enable "mantid banana ... list all neutrons" in the simulation.' + ) + values = data.get_event_data( + variables=['p', 't', 'id'], component_name=detector_name, filter_zeros=True + ) + if geometry is None: + geometry = _detector_geometry(data, detector_name) + pixel_ids = geometry.coords['pixel_id'] + if not np.isin(values['id'], pixel_ids.values).all(): + raise ValueError('Detector events contain pixel IDs absent from the pixel map.') + time = sc.array(dims=['event'], values=values['t'], unit='s') + offset = time % pulse_period.to(unit=time.unit) + events = sc.DataArray( + sc.array( + dims=['event'], + values=values['p'], + variances=values['p'] ** 2, + unit='counts', + ), + coords={ + 'pixel_id': sc.array( + dims=['event'], values=values['id'], dtype='int64', unit=None + ), + 'event_time_offset': offset, + 'event_time_zero': sc.datetime(0, unit='ns') + + (time - offset).to(unit='ns', dtype='int64'), + }, + ).group(pixel_ids) + return events.assign_coords(geometry.coords) + + +def load_mcstas_provider( + filename: Filename[RunType], + detector_name: NeXusDetectorName, + geometry: EmptyDetector[RunType], + pulse_period: PulsePeriod, +) -> RawDetector[RunType]: + """Provide final-detector events, reusing the workflow's pixel geometry.""" + with _open_mcstas(filename) as data: + return RawDetector[RunType]( + _load_events(data, detector_name, pulse_period, geometry) + ) + + +def _component_position(filename, component): + with _open_mcstas(filename) as data: + return sc.vector(data.get_global_component_coordinates(component), unit='m') + + +def mcstas_source_position( + filename: Filename[RunType], +) -> Position[snx.NXsource, RunType]: + """Load the moderator position from the simulation geometry.""" + return Position[snx.NXsource, RunType](_component_position(filename, 'Source')) + + +def mcstas_sample_position( + filename: Filename[RunType], +) -> Position[snx.NXsample, RunType]: + """Load the sample position from the simulation geometry.""" + return Position[snx.NXsample, RunType](_component_position(filename, 'Arm_Sample')) + + +def _histogram_axis(histogram, axis, unit): + label = _text(histogram.attrs[f'{axis}label']) + name = re.sub(r'[^a-zA-Z]', '_', label) + # McStas encodes the units in the axis label, e.g., "y [m]". + axis_unit = label[label.index('[') + 1 : label.index(']')] + return ( + sc.array( + dims=['pixel'], values=histogram[name][:], unit=axis_unit, dtype='float64' + ) + .to(unit=unit) + .values + ) + + +def _detector_geometry(data, detector_name) -> sc.DataArray: + output = data.file_object.get_output_entry(detector_name) + if 'BINS' in output: + # Keep the file's ID order, including non-contiguous or permuted IDs. + local = data.get_component_local(detector_name) + position = data.get_component_global(detector_name) + pixel_ids = np.asarray( + data.file_object.get_pixels_entry(detector_name), dtype='int64' + ).ravel() + else: + # Histogram axes also describe geometry; no intensities are read. + geometry = data.file_object.get_geometry_dict(detector_name) + if geometry['shape'] != 'banana': + raise ValueError(f'Expected banana geometry for {detector_name!r}.') + histogram = data.file_object.get_info_entry(detector_name) + angle, height = np.meshgrid( + _histogram_axis(histogram, 'x', 'rad'), + _histogram_axis(histogram, 'y', 'm'), + ) + radius = geometry['radius'] + local = np.column_stack( + ( + (radius * np.sin(angle)).ravel(), + height.ravel(), + (radius * np.cos(angle)).ravel(), + ) + ) + pixel_ids = np.arange(len(local)) + position = data.transform(local, detector_name) + return sc.DataArray( + sc.zeros(dims=['pixel_id'], shape=[len(local)], unit='counts'), + coords={ + 'pixel_id': sc.array(dims=['pixel_id'], values=pixel_ids, unit=None), + 'position': sc.vectors(dims=['pixel_id'], values=position, unit='m'), + 'longitude': sc.array( + dims=['pixel_id'], + values=np.rad2deg(np.arctan2(local[:, 0], local[:, 2])), + unit='deg', + ), + 'height': sc.array(dims=['pixel_id'], values=local[:, 1], unit='m'), + }, + ) + + +def mcstas_detector_geometry( + filename: Filename[RunType], detector_name: NeXusDetectorName +) -> EmptyDetector[RunType]: + """Provide detector geometry to the generic flight-path calculation. + + Works with Mantid pixel maps or banana histogram axes without reading events + or intensities. The selected component's name and axes are read from metadata. + """ + with _open_mcstas(filename) as data: + return EmptyDetector[RunType](_detector_geometry(data, detector_name)) + + +providers = ( + wfm_choppers, + load_mcstas_provider, + mcstas_source_position, + mcstas_sample_position, + mcstas_detector_geometry, +) diff --git a/packages/essreflectometry/src/ess/freia/workflow.py b/packages/essreflectometry/src/ess/freia/workflow.py index 327a5f9d9..fae23f746 100644 --- a/packages/essreflectometry/src/ess/freia/workflow.py +++ b/packages/essreflectometry/src/ess/freia/workflow.py @@ -29,7 +29,7 @@ ) from .corrections import RunNormalization, insert_run_normalization -_general_providers = ( +providers = ( *reflectometry_providers, *conversions.providers, *corrections.providers, @@ -39,37 +39,23 @@ *load.providers, ) -mcstas_providers = ( - *_general_providers, - *mcstas.providers, -) -"""List of providers for setting up a Sciline pipeline for McStas data. - -This provides a default Freia workflow including providers for loadings files. -""" - -providers = (*_general_providers,) """List of providers for setting up a Sciline pipeline data. -This provides a default Freia workflow including providers for loadings files. +This provides a default Freia workflow including providers for loading files. """ def mcstas_default_parameters() -> dict: """Return default parameters for the McStas Freia workflow.""" - return { - DetectorSpatialResolution: 0.0025 * sc.units.m, - NeXusDetectorName: "detector", + return default_parameters() | { + NeXusDetectorName: "Multiblade", BeamDivergenceLimits: ( sc.scalar(-0.75, unit='deg'), sc.scalar(0.75, unit='deg'), ), - SampleRotationOffset[RunType]: sc.scalar(0.0, unit='deg'), - CorrectionsToApply: corrections.default_corrections, LookupTableRelativeErrorThreshold: { - "detector": 0.06, + "Multiblade": 0.06, }, - UncertaintyBroadcastMode: UncertaintyBroadcastMode.drop, } @@ -90,10 +76,17 @@ def default_parameters() -> dict: def FreiaMcStasWorkflow( *, run_norm: RunNormalization = RunNormalization.none, - wavelength_from: WavelengthLutMode = "file", + wavelength_from: WavelengthLutMode = "analytical", **kwargs, ) -> sciline.Pipeline: - """Workflow for reduction of McStas data for the Freia instrument. + """Workflow for loading and unwrapping FREIA McStas detector events. + + Builds on the generic NeXus/unwrapping workflow, replacing input providers + with adapters for the final ``Multiblade`` Mantid banana detector. Detector + geometry is read from the simulation file; choppers use the fixed WFM + configuration from :func:`ess.freia.mcstas.wfm_choppers`. This initial + workflow supports visualization through ``WavelengthDetector``; full + reflectivity reduction and normalization require further instrument providers. Parameters ---------- @@ -104,10 +97,11 @@ def FreiaMcStasWorkflow( 'analytical', 'simulation', and 'file'. See https://scipp.github.io/ess/reduce/user-guide/unwrap/lut-building-methods.html """ - workflow = beamline.LoadNeXusWorkflow(wavelength_from=wavelength_from, **kwargs) - for provider in mcstas_providers: + workflow = FreiaWorkflow( + run_norm=run_norm, wavelength_from=wavelength_from, **kwargs + ) + for provider in mcstas.providers: workflow.insert(provider) - insert_run_normalization(workflow, run_norm) for name, param in mcstas_default_parameters().items(): workflow[name] = param return workflow diff --git a/packages/essreflectometry/tests/freia/mcstas_test.py b/packages/essreflectometry/tests/freia/mcstas_test.py new file mode 100644 index 000000000..d7c2d01ec --- /dev/null +++ b/packages/essreflectometry/tests/freia/mcstas_test.py @@ -0,0 +1,226 @@ +# SPDX-License-Identifier: BSD-3-Clause +# Copyright (c) 2026 Scipp contributors (https://github.com/scipp) + +import math + +import h5py +import numpy as np +import pytest +import scipp as sc +from ess.reduce.nexus.types import ( + DiskChoppers, + Filename, + NeXusDetectorName, + RawDetector, +) +from ess.reduce.unwrap import ( + DetectorLtotal, + FrameUnwrapBackend, + LookupTable, + PulsePeriod, + SourceBounds, + WavelengthDetector, +) +from scipp.testing import assert_allclose, assert_identical +from scippneutron.chopper import DiskChopper +from scippnexus import NXdetector + +from ess.freia import FreiaMcStasWorkflow +from ess.freia.mcstas import load_mcstas +from ess.reflectometry.types import ReferenceRun, SampleRun + + +def _component(components, name, position): + group = components.create_group(f'{len(components):04d}_{name}') + group['Position'] = position + group['Rotation'] = np.eye(3) + return group + + +@pytest.fixture +def mcstas_file(tmp_path): + """Small on-disk McStas file read by the real mcstastox library.""" + filename = tmp_path / 'freia.h5' + with h5py.File(filename, 'w') as f: + entry = f.create_group('entry1') + entry.create_group('data') + simulation = entry.create_group('simulation') + simulation.attrs['program'] = np.bytes_('3.7.18, git') + simulation.create_group('Param') + instrument = entry.create_group('instrument') + components = instrument.create_group('components') + _component(components, 'Source', [0.0, 0.0, 0.0]) + _component(components, 'Arm_Sample', [0.0, 0.0, 20.0]) + # A debug monitor must never be included in the detector event sum. + debug = _component(components, 'Slit_event', [0.0, 0.0, 19.0]) + debug_output = debug.create_group('output/events') + debug_output.attrs['variables'] = np.bytes_('p t id') + debug_output['events'] = [[999.0, 0.001, 0.0]] + detector = _component(components, 'Multiblade', [0.0, 0.0, 20.0]) + # Rotate the banana into the vertical scattering plane. + detector['Rotation'][...] = [[0, 1, 0], [-1, 0, 0], [0, 0, 1]] + geometry = detector.create_group('Geometry') + geometry.attrs['Shape identifier'] = np.bytes_('4') + geometry.attrs['radius'] = np.bytes_('3') + bins = detector.create_group('output/BINS') + for key, value in { + 'xvar': 'th', + 'yvar': 'y', + 'xlabel': 'theta', + 'ylabel': 'height', + }.items(): + bins.attrs[key] = np.bytes_(value) + bins['theta'] = [1.0, 2.0] + bins['height'] = [-0.001, 0.001] + # IDs need not be contiguous or sorted in geometry order. + bins['pixels'] = [[12, 10], [99, 101]] + output = detector.create_group('output/detector_events') + output.attrs['variables'] = np.bytes_('p t id') + output['events'] = [ + [2.0, 0.025, 10.0], + [3.0, 0.025 + 1 / 14, 10.0], + [0.0, 0.025, 99.0], + [4.0, 0.025, 12.0], + ] + return filename + + +def test_loader_reads_only_final_detector_and_retains_empty_pixels(mcstas_file): + detector = load_mcstas(mcstas_file) + # Check pixel associations without requiring the loader to return a given order. + image = sc.sort(detector.bins.sum(), 'pixel_id') + np.testing.assert_array_equal(image.coords['pixel_id'].values, [10, 12, 99, 101]) + assert_identical( + image.data, + sc.array( + dims=['pixel_id'], + values=[5.0, 4.0, 0.0, 0.0], + variances=[13.0, 16.0, 0.0, 0.0], + unit='counts', + ), + ) + offsets = detector.bins.constituents['data'].coords['event_time_offset'] + assert_allclose( + offsets.to(unit='s'), + sc.full(sizes=offsets.sizes, value=0.025, unit='s'), + ) + # Pixel 12 is the first banana pixel, rotated into the global frame. + assert_allclose( + image.coords['position'][1], + sc.vector( + [0.001, 3 * math.sin(math.pi / 180), 20 + 3 * math.cos(math.pi / 180)], + unit='m', + ), + ) + + +def test_loader_rejects_missing_detector_without_using_debug_events(mcstas_file): + with h5py.File(mcstas_file, 'r+') as f: + del f['entry1/instrument/components/0003_Multiblade'] + with pytest.raises(ValueError, match='No Mantid detector events'): + load_mcstas(mcstas_file) + + +@pytest.mark.parametrize('pixel_id', [999, 10.5]) +def test_loader_rejects_events_missing_from_pixel_map(mcstas_file, pixel_id): + with h5py.File(mcstas_file, 'r+') as f: + events = f[ + 'entry1/instrument/components/0003_Multiblade/output/detector_events/events' + ] + events[0, 2] = pixel_id + with pytest.raises(ValueError, match='pixel IDs absent from the pixel map'): + load_mcstas(mcstas_file) + + +def test_lookup_table_uses_histogram_geometry_without_loading_events(mcstas_file): + with h5py.File(mcstas_file, 'r+') as f: + components = f['entry1/instrument/components'] + components.move('0003_Multiblade', '0003_Detector') + detector = components['0003_Detector'] + detector['Position'][...] = [0.0, -0.25, 20.0] + del detector['output'] + histogram = detector.create_group('output/histogram') + histogram.attrs['xlabel'] = np.bytes_('theta [deg]') + histogram.attrs['ylabel'] = np.bytes_('Height [cm]') + histogram['theta__deg_'] = [0.0, 15.0] + histogram['Height__cm_'] = [-25, 25] + # No intensities or event arrays are needed anywhere in the file. + del components['0002_Slit_event'] + + workflow = FreiaMcStasWorkflow() + workflow[Filename[SampleRun]] = str(mcstas_file) + workflow[NeXusDetectorName] = 'Detector' + results = workflow.compute( + (DetectorLtotal[SampleRun], LookupTable[SampleRun, NXdetector]) + ) + # The detector's vertical offset shortens the distance at larger angles. + expected = [ + 20 + math.sqrt(9.125 - 1.5 * math.sin(angle)) for angle in (0, math.pi / 12) + ] + assert_allclose( + results[DetectorLtotal[SampleRun]], + sc.array(dims=['pixel_id'], values=expected * 2, unit='m'), + ) + table = results[LookupTable[SampleRun, NXdetector]].array + assert sc.isfinite(table.data).any().value + # The histogram geometry must not make it possible to load fake events. + with pytest.raises(ValueError, match='No Mantid detector events'): + workflow.compute(RawDetector[SampleRun]) + + +@pytest.mark.parametrize('run', [SampleRun, ReferenceRun]) +def test_workflow_loads_and_unwraps_detector_with_generic_providers(mcstas_file, run): + workflow = FreiaMcStasWorkflow() + workflow[Filename[run]] = str(mcstas_file) + workflow[FrameUnwrapBackend] = FrameUnwrapBackend.scipy + # Override WFM with a single disk for an independently calculable wavelength. + workflow[DiskChoppers[run]] = { + 'test_chopper': DiskChopper( + frequency=sc.scalar(14.0, unit='Hz'), + beam_position=sc.scalar(0.0, unit='deg'), + phase=sc.scalar(0.0, unit='deg'), + axle_position=sc.vector([0.0, 0.0, 5.0], unit='m'), + slit_begin=sc.array(dims=['cutout'], values=[279.68], unit='deg'), + slit_end=sc.array(dims=['cutout'], values=[359.68], unit='deg'), + ), + } + workflow[SourceBounds] = SourceBounds( + time=(sc.scalar(0.9, unit='ms'), sc.scalar(1.1, unit='ms')), + wavelength=(sc.scalar(0.5, unit='angstrom'), sc.scalar(12.0, unit='angstrom')), + ) + result = workflow.compute((RawDetector[run], WavelengthDetector[run])) + raw = result[RawDetector[run]] + unwrapped = result[WavelengthDetector[run]] + assert_allclose(raw.bins.sum().data, unwrapped.bins.sum().data) + wavelength = unwrapped.bins.constituents['data'].coords['wavelength'] + # Arrival time 25 ms minus emission time 1 ms, flight path about 23 m. + expected = ( + sc.constants.h + / sc.constants.m_n + * sc.scalar(24.0, unit='ms') + / sc.scalar(23.0, unit='m') + ).to(unit='angstrom') + assert_allclose( + wavelength, + sc.full(sizes=wavelength.sizes, value=expected.value, unit=expected.unit), + rtol=sc.scalar(0.002), + ) + + +def test_loader_uses_workflow_pulse_period(mcstas_file): + workflow = FreiaMcStasWorkflow() + workflow[Filename[SampleRun]] = str(mcstas_file) + workflow[PulsePeriod] = sc.scalar(50.0, unit='ms') + events = workflow.compute(RawDetector[SampleRun]).bins.constituents['data'] + events = sc.sort(events, 'event_time_zero') + # The event at 25 ms + 1/14 s falls in the second 50 ms pulse period. + assert_identical( + events.coords['event_time_zero'].to(unit='ns'), + sc.datetimes(dims=['event'], values=[0, 0, 50_000_000], unit='ns'), + ) + assert_allclose( + events.coords['event_time_offset'].to(unit='s'), + sc.array( + dims=['event'], values=[0.025, 0.025, 0.025 + 1 / 14 - 0.05], unit='s' + ), + ) diff --git a/packages/essreflectometry/tests/freia/workflow_test.py b/packages/essreflectometry/tests/freia/workflow_test.py index 41bc7a887..fbf869f16 100644 --- a/packages/essreflectometry/tests/freia/workflow_test.py +++ b/packages/essreflectometry/tests/freia/workflow_test.py @@ -1,30 +1,8 @@ # SPDX-License-Identifier: BSD-3-Clause # Copyright (c) 2026 Scipp contributors (https://github.com/scipp) -import inspect - from ess.reduce import workflow as reduce_workflow from ess import freia -from ess.freia.corrections import RunNormalization -from ess.reflectometry.types import CorrectionsToApply, NeXusDetectorName - - -def test_freia_workflow_uses_freia_defaults(): - params = freia.workflow.default_parameters() - - assert params[NeXusDetectorName] == "multiblade_detector" - assert params[CorrectionsToApply] == freia.corrections.default_corrections - - -def test_freia_workflows_have_expected_run_normalization_defaults(): - assert ( - inspect.signature(freia.FreiaMcStasWorkflow).parameters["run_norm"].default - is RunNormalization.none - ) - assert ( - inspect.signature(freia.FreiaWorkflow).parameters["run_norm"].default - is RunNormalization.proton_charge - ) def test_freia_workflow_registers_run_normalization_variants(): From 38cadda57c58e8496980a7f1e85528aa54b14b71 Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Fri, 11 Sep 2026 16:26:30 +0200 Subject: [PATCH 02/16] docs: fix tone and content, remove unnecessary --- .../freia/freia-mcstas-visualization.ipynb | 16 +++----- .../freia/freia-wavelength-lookup-table.ipynb | 41 ++----------------- .../docs/user-guide/freia/index.md | 10 ++--- 3 files changed, 14 insertions(+), 53 deletions(-) diff --git a/packages/essreflectometry/docs/user-guide/freia/freia-mcstas-visualization.ipynb b/packages/essreflectometry/docs/user-guide/freia/freia-mcstas-visualization.ipynb index 8727d1cd7..b4cba128d 100644 --- a/packages/essreflectometry/docs/user-guide/freia/freia-mcstas-visualization.ipynb +++ b/packages/essreflectometry/docs/user-guide/freia/freia-mcstas-visualization.ipynb @@ -7,11 +7,7 @@ "source": [ "# FREIA McStas detector data\n", "\n", - "This notebook loads the final FREIA detector with `mcstastox`, visualizes weighted events with Scipp and Plopp, and computes wavelengths with ESSreduce's analytical frame-unwrapping workflow. It is an initial detector-data inspection workflow; reflectivity normalization and conversion to Q are later steps.\n", - "\n", - "Use an environment containing the local `essreflectometry` and `essreduce` packages, Jupyter, and `ipympl` (the reflectometry documentation environment includes the plotting dependencies).\n", - "\n", - "The input must contain the **final `Multiblade` Mantid banana detector**, including an event list and pixel geometry. Enable its `mantid banana ... list all neutrons` output when running McStas. A `Multiblade_histogram` output alone has no event times and cannot be unwrapped. The loader reports an error if the detector events are absent.\n" + "This notebook visualizes the final detector's image, arrival-time distribution, and wavelength spectrum for a FREIA simulation. Wavelengths are reconstructed using the WFM chopper settings.\n" ] }, { @@ -46,7 +42,7 @@ "source": [ "## Select a run\n", "\n", - "Load the sample run using `ess.freia.data`.\n" + "Select the sample simulation and set the resolution used to reconstruct wavelengths.\n" ] }, { @@ -60,7 +56,7 @@ "freia_mcstas[Filename[SampleRun]] = data.freia_mcstas_sample_run()\n", "freia_mcstas[KeepEventTimeOffset] = True\n", "freia_mcstas[TimeResolution] = sc.scalar(20.0, unit='us')\n", - "freia_mcstas[DistanceResolution] = sc.scalar(0.01, unit='m')\n" + "freia_mcstas[DistanceResolution] = sc.scalar(0.1, unit='m')\n" ] }, { @@ -70,7 +66,7 @@ "source": [ "## Inspect the WFM chopper cascade\n", "\n", - "`freia.mcstas.wfm_choppers()` constructs the fixed **WFM** configuration used by this workflow. The configuration contains three bandwidth disks and five WFM disks with seven openings each. It is independent of the input filename.\n", + "This example uses a fixed **WFM** configuration with three bandwidth disks and five WFM disks with seven openings each.\n", "\n", "For another configuration, assign its chopper dictionary to `freia_mcstas[DiskChoppers[SampleRun]]`.\n", "\n", @@ -106,7 +102,7 @@ "source": [ "## Load and unwrap the detector events\n", "\n", - "The generic workflow combines the McStas input providers with ESSreduce's wavelength calculation. We request both outputs together so that the file's detector events are loaded once. Weights are retained, with per-event variances equal to the squared weights; empty pixels are also retained.\n" + "Compute wavelengths from the event arrival times and chopper transmission bands.\n" ] }, { @@ -129,7 +125,7 @@ "source": [ "## Detector image and arrival times\n", "\n", - "`longitude` and `height` are coordinates in the banana detector's local frame; the `position` vectors use the global instrument frame. All image intensities are sums of event weights.\n" + "The image uses `longitude` and `height` in the detector's local frame. Intensities are sums of event weights.\n" ] }, { diff --git a/packages/essreflectometry/docs/user-guide/freia/freia-wavelength-lookup-table.ipynb b/packages/essreflectometry/docs/user-guide/freia/freia-wavelength-lookup-table.ipynb index 34a85e21b..7c4eb9c9a 100644 --- a/packages/essreflectometry/docs/user-guide/freia/freia-wavelength-lookup-table.ipynb +++ b/packages/essreflectometry/docs/user-guide/freia/freia-wavelength-lookup-table.ipynb @@ -7,7 +7,7 @@ "source": [ "# FREIA analytical wavelength lookup table\n", "\n", - "Construct and visualize a wavelength lookup table using FREIA's fixed WFM chopper configuration and ESSreduce's generic analytical frame-unwrapping workflow.\n", + "Explore which neutron wavelengths can reach the detector at each arrival time for the WFM chopper configuration.\n", "\n", "The calculation uses the final detector's pixel geometry and the source and sample positions. Detector events and histogram intensities are not needed to construct the lookup table.\n" ] @@ -20,8 +20,6 @@ "outputs": [], "source": [ "%matplotlib widget\n", - "from pathlib import Path\n", - "\n", "import plopp as pp\n", "import scipp as sc\n", "from scippnexus import NXdetector\n", @@ -34,7 +32,6 @@ " DetectorLtotal,\n", " DistanceResolution,\n", " LookupTable,\n", - " SourceBounds,\n", " TimeResolution,\n", ")\n", "from ess.reflectometry.types import SampleRun\n" @@ -47,9 +44,9 @@ "source": [ "## Use the WFM configuration\n", "\n", - "The workflow uses `freia.mcstas.wfm_choppers()` to construct the fixed **WFM** cascade. Detector geometry is read from the selected file, and ESSreduce computes the wavelength table analytically.\n", + "Use the fixed **WFM** chopper settings and read detector geometry from the sample simulation.\n", "\n", - "This example uses 20 µs time resolution and 5 mm flight-path resolution. Source bounds are the generic ESS defaults, shown explicitly below. The central-ray model projects chopper positions onto the global z axis and approximates detector flight paths by source-to-sample plus sample-to-pixel distance.\n" + "This example uses 20 µs time resolution and 10 cm flight-path resolution. The central-ray model projects chopper positions onto the global z axis and approximates detector flight paths by source-to-sample plus sample-to-pixel distance.\n" ] }, { @@ -64,11 +61,7 @@ "freia_mcstas = freia.FreiaMcStasWorkflow(wavelength_from='analytical')\n", "freia_mcstas[Filename[SampleRun]] = filename\n", "freia_mcstas[TimeResolution] = sc.scalar(20.0, unit='us')\n", - "freia_mcstas[DistanceResolution] = sc.scalar(0.005, unit='m')\n", - "freia_mcstas[SourceBounds] = SourceBounds(\n", - " time=(sc.scalar(0.0, unit='ms'), sc.scalar(5.0, unit='ms')),\n", - " wavelength=(sc.scalar(0.001, unit='angstrom'), sc.scalar(15.0, unit='angstrom')),\n", - ")\n" + "freia_mcstas[DistanceResolution] = sc.scalar(0.1, unit='m')\n" ] }, { @@ -149,32 +142,6 @@ "figure[1, 0].ax.grid(alpha=0.2)\n", "figure\n" ] - }, - { - "cell_type": "markdown", - "id": "7", - "metadata": {}, - "source": [ - "## Save the lookup table and figure\n", - "\n", - "The HDF5 file retains the table, variances, pulse period, and interpolation resolutions. The PNG and PDF are standalone figures. Outputs are written to `build/freia` relative to the current working directory.\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "8", - "metadata": {}, - "outputs": [], - "source": [ - "output_dir = Path('build') / 'freia'\n", - "output_dir.mkdir(parents=True, exist_ok=True)\n", - "stem = f'{filename.stem}-wavelength-lookup-table'\n", - "lookup.save_hdf5(output_dir / f'{stem}.h5')\n", - "figure.save(str(output_dir / f'{stem}.png'), dpi=180, bbox_inches='tight')\n", - "figure.save(str(output_dir / f'{stem}.pdf'), bbox_inches='tight')\n", - "output_dir\n" - ] } ], "metadata": { diff --git a/packages/essreflectometry/docs/user-guide/freia/index.md b/packages/essreflectometry/docs/user-guide/freia/index.md index e51482d56..f2b7dbb96 100644 --- a/packages/essreflectometry/docs/user-guide/freia/index.md +++ b/packages/essreflectometry/docs/user-guide/freia/index.md @@ -1,11 +1,9 @@ # FREIA -The initial FREIA workflow loads final-detector McStas events and uses the generic -ESSreduce analytical frame-unwrapping workflow to compute wavelengths. The notebook -guides below use the example download helpers in `ess.freia.data` and visualize the -detector with Scipp and Plopp. Choppers use the fixed **WFM** simulation configuration. -The wavelength lookup-table guide reads only detector geometry from the file, so it -can also be run on simulations without detector events. +Explore FREIA simulations through detector images, arrival-time distributions, and +wavelength spectra. The guides use fixed **WFM** chopper settings to reconstruct +wavelengths. The wavelength lookup-table guide requires only detector geometry and +the source and sample positions; detector events are not needed. ```{toctree} :maxdepth: 1 From d5a68508a88dfb0c98f7a24dc0fb7fa8922832e3 Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Mon, 14 Sep 2026 16:31:07 +0200 Subject: [PATCH 03/16] feat: normalizations by monitor, footprint and direct beam --- .../freia/freia-mcstas-visualization.ipynb | 40 ++- .../user-guide/freia/freia-reflectivity.ipynb | 266 ++++++++++++++++++ .../freia/freia-wavelength-lookup-table.ipynb | 12 +- .../docs/user-guide/freia/index.md | 7 +- .../src/ess/freia/beamline.py | 4 +- .../src/ess/freia/conversions.py | 209 +++++--------- .../src/ess/freia/corrections.py | 77 ++--- .../src/ess/freia/maskings.py | 75 ++--- .../essreflectometry/src/ess/freia/mcstas.py | 66 ++++- .../src/ess/freia/normalization.py | 71 ++++- .../essreflectometry/src/ess/freia/types.py | 22 ++ .../src/ess/freia/workflow.py | 25 +- .../tests/freia/mcstas_test.py | 38 ++- .../tests/freia/workflow_test.py | 200 +++++++++++++ 14 files changed, 847 insertions(+), 265 deletions(-) create mode 100644 packages/essreflectometry/docs/user-guide/freia/freia-reflectivity.ipynb diff --git a/packages/essreflectometry/docs/user-guide/freia/freia-mcstas-visualization.ipynb b/packages/essreflectometry/docs/user-guide/freia/freia-mcstas-visualization.ipynb index b4cba128d..c048916af 100644 --- a/packages/essreflectometry/docs/user-guide/freia/freia-mcstas-visualization.ipynb +++ b/packages/essreflectometry/docs/user-guide/freia/freia-mcstas-visualization.ipynb @@ -5,9 +5,9 @@ "id": "0", "metadata": {}, "source": [ - "# FREIA McStas detector data\n", + "# FREIA detector data\n", "\n", - "This notebook visualizes the final detector's image, arrival-time distribution, and wavelength spectrum for a FREIA simulation. Wavelengths are reconstructed using the WFM chopper settings.\n" + "Visualize detector images, arrival-time distributions, and wavelength spectra. Wavelengths are reconstructed using the WFM chopper settings.\n" ] }, { @@ -42,7 +42,7 @@ "source": [ "## Select a run\n", "\n", - "Select the sample simulation and set the resolution used to reconstruct wavelengths.\n" + "Select a run and set the resolution used to reconstruct wavelengths.\n" ] }, { @@ -66,11 +66,9 @@ "source": [ "## Inspect the WFM chopper cascade\n", "\n", - "This example uses a fixed **WFM** configuration with three bandwidth disks and five WFM disks with seven openings each.\n", + "Inspect the chopper timing used to reconstruct wavelengths. To use another configuration, assign its chopper dictionary to `freia_mcstas[DiskChoppers[SampleRun]]`.\n", "\n", - "For another configuration, assign its chopper dictionary to `freia_mcstas[DiskChoppers[SampleRun]]`.\n", - "\n", - "The analytical model approximates the beam by a central ray. It projects chopper positions onto the global z axis and estimates detector flight paths as source-to-sample plus sample-to-pixel distance; it does not trace the curved guide or model the finite beam footprint on each disk.\n" + "The analytical model approximates neutron flight paths.\n" ] }, { @@ -135,8 +133,8 @@ "metadata": {}, "outputs": [], "source": [ - "detector_image = raw.hist(longitude=120, height=80, dim=raw.dims)\n", - "pp.plot(detector_image, norm='log', title='FREIA final detector')\n" + "detector_image = raw.hist(longitude=120, height=64, dim=raw.dims)\n", + "pp.plot(detector_image, norm='log', title='FREIA final detector', vmin=1e-1)\n" ] }, { @@ -190,7 +188,18 @@ "wavelength_image = unwrapped.hist(\n", " wavelength=wavelength_bins, longitude=120, dim=unwrapped.dims\n", ")\n", - "pp.plot(wavelength_image, norm='log', title='Wavelength across the detector')\n" + "pp.plot(wavelength_image, norm='log', title='Wavelength across the detector', vmin=1e0)\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "15", + "metadata": {}, + "outputs": [], + "source": [ + "unwrapped.bins.coords['wavelength_error'] = (unwrapped.bins.coords['wavelength'] - unwrapped.bins.coords['wavelength_from_mcstas']) / unwrapped.bins.coords['wavelength_from_mcstas']\n", + "unwrapped.hist(wavelength_error=100, dim=unwrapped.dims).plot()" ] } ], @@ -201,7 +210,16 @@ "name": "python3" }, "language_info": { - "name": "python" + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.12.13" } }, "nbformat": 4, diff --git a/packages/essreflectometry/docs/user-guide/freia/freia-reflectivity.ipynb b/packages/essreflectometry/docs/user-guide/freia/freia-reflectivity.ipynb new file mode 100644 index 000000000..9decd8e4b --- /dev/null +++ b/packages/essreflectometry/docs/user-guide/freia/freia-reflectivity.ipynb @@ -0,0 +1,266 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "0", + "metadata": {}, + "source": [ + "# FREIA reflectivity reduction\n", + "\n", + "Reduce a reflected beam using a direct-beam measurement without a sample. The workflow converts events to specular Q, normalizes by an incident wavelength monitor, and divides the selected peaks to obtain R(Q).\n", + "\n", + "This example uses local McStas files. The sample and direct-beam runs must have matching slit and chopper settings." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1", + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib inline\n", + "import scipp as sc\n", + "\n", + "from ess.freia import FreiaMcStasWorkflow\n", + "from ess.freia.corrections import RunNormalization\n", + "from ess.freia.types import (\n", + " DetectorRegionOfInterest,\n", + " IncidentMonitor,\n", + " QDetector,\n", + " SampleIlluminatedFraction,\n", + " WavelengthMonitor,\n", + ")\n", + "from ess.reduce.nexus.types import NeXusName\n", + "from ess.reflectometry.types import (\n", + " BeamSize,\n", + " Filename,\n", + " QBins,\n", + " ReducibleData,\n", + " ReferenceRun,\n", + " ReflectivityOverQ,\n", + " SampleRun,\n", + " SampleSize,\n", + " WavelengthBins,\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "2", + "metadata": {}, + "source": [ + "## Select the runs\n", + "\n", + "Set the sample and no-sample filenames. These paths are relative to the notebook directory. Choose an incident wavelength monitor upstream of the sample, so its intensity is independent of reflectivity. If no monitor is available, use `RunNormalization.none`; the two runs must then already share an exposure and flux scale." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3", + "metadata": {}, + "outputs": [], + "source": [ + "workflow = FreiaMcStasWorkflow(run_norm=RunNormalization.monitor_histogram)\n", + "workflow[Filename[SampleRun]] = '../../../265149.h5'\n", + "workflow[Filename[ReferenceRun]] = '../../../265148.h5'\n", + "workflow[NeXusName[IncidentMonitor]] = 'GuideexitLambda'\n", + "\n", + "workflow[WavelengthBins] = sc.linspace('wavelength', 2.0, 10.0, 81, unit='angstrom')\n", + "workflow[QBins] = sc.geomspace('Q', 0.07, 0.4, 21, unit='1/angstrom')" + ] + }, + { + "cell_type": "markdown", + "id": "4", + "metadata": {}, + "source": [ + "## Inspect the reflected and direct peaks\n", + "\n", + "The angle θ is measured above the sample surface and includes gravity correction. Reflected beams have positive angles; the direct beams lie below the sample plane. The specular assumption gives $Q = 4\\pi\\sin(\\theta)/\\lambda$.\n", + "\n", + "Inspect both distributions before selecting corresponding peaks." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5", + "metadata": {}, + "outputs": [], + "source": [ + "detectors = workflow.compute((QDetector[SampleRun], QDetector[ReferenceRun]))\n", + "angle_bins = sc.linspace('theta', -4.0, 4.0, 161, unit='deg').to(unit='rad')\n", + "profiles = {}\n", + "for label, run in [('Sample', SampleRun), ('Direct beam', ReferenceRun)]:\n", + " detector = detectors[QDetector[run]]\n", + " profile = detector.hist(theta=angle_bins, dim=detector.dims)\n", + " profile.coords['theta'] = profile.coords['theta'].to(unit='deg')\n", + " profiles[label] = profile\n", + "sc.plot(profiles, norm='log', title='Reflected and direct beams', vmin=1e4, vmax=1e8)" + ] + }, + { + "cell_type": "markdown", + "id": "6", + "metadata": {}, + "source": [ + "## Select matching regions of interest\n", + "\n", + "Select the reflected peak and its corresponding direct peak separately. The example uses the beam near 3.5°. Additional bounds can be set on pixel coordinates such as `height` or `pixel_id`. Wavelength and Q bins are deliberately coarse because these runs have few events." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7", + "metadata": {}, + "outputs": [], + "source": [ + "workflow[DetectorRegionOfInterest[SampleRun]] = {\n", + " 'theta': (sc.scalar(3.2, unit='deg'), sc.scalar(3.9, unit='deg')),\n", + "}\n", + "workflow[DetectorRegionOfInterest[ReferenceRun]] = {\n", + " 'theta': (sc.scalar(-3.9, unit='deg'), sc.scalar(-3.2, unit='deg')),\n", + "}" + ] + }, + { + "cell_type": "markdown", + "id": "8", + "metadata": {}, + "source": [ + "## Inspect the wavelength normalization\n", + "\n", + "The incident monitor corrects wavelength-dependent flux differences between runs. Its wavelength range must cover the detector data." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9", + "metadata": {}, + "outputs": [], + "source": [ + "monitors = workflow.compute(\n", + " (WavelengthMonitor[SampleRun], WavelengthMonitor[ReferenceRun])\n", + ")\n", + "sc.plot(\n", + " {\n", + " 'Sample monitor': monitors[WavelengthMonitor[SampleRun]],\n", + " 'Direct-beam monitor': monitors[WavelengthMonitor[ReferenceRun]],\n", + " },\n", + " title='Incident wavelength spectra',\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "10", + "metadata": {}, + "outputs": [], + "source": [ + "selected = workflow.compute((ReducibleData[SampleRun], ReducibleData[ReferenceRun]))\n", + "spectra = {}\n", + "for label, run in [('Sample', SampleRun), ('Direct beam', ReferenceRun)]:\n", + " detector = selected[ReducibleData[run]]\n", + " spectra[label] = detector.hist(\n", + " wavelength=workflow.compute(WavelengthBins), dim=detector.dims\n", + " )\n", + "sc.plot(spectra, title='Selected peaks after monitor normalization')" + ] + }, + { + "cell_type": "markdown", + "id": "11", + "metadata": {}, + "source": [ + "## Configure footprint correction\n", + "\n", + "The Gaussian footprint model uses the sample length along the beam and the beam FWHM at the sample. It corrects only the reflected run. Different incident beams may require different widths.\n", + "\n", + "A beam width has not yet been established for these data, so this example leaves the footprint correction off. Replace `beam_size = None` with a measured width, for example `sc.scalar(width_in_mm, unit='mm')`, to enable it." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "12", + "metadata": {}, + "outputs": [], + "source": [ + "sample_size = sc.scalar(80.0, unit='mm')\n", + "beam_size = None\n", + "\n", + "reduction = workflow.copy()\n", + "if beam_size is None:\n", + " reduction[SampleIlluminatedFraction] = sc.scalar(1.0)\n", + "else:\n", + " reduction[SampleSize[SampleRun]] = sample_size\n", + " reduction[BeamSize[SampleRun]] = beam_size" + ] + }, + { + "cell_type": "markdown", + "id": "13", + "metadata": {}, + "source": [ + "## Compute reflectivity\n", + "\n", + "The workflow maps the direct beam to its corresponding specular Q, integrates the two selected peaks into matching Q bins, and divides their intensities. Both counting uncertainties propagate. Bins with no usable direct-beam intensity are masked." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "14", + "metadata": {}, + "outputs": [], + "source": [ + "reflectivity = reduction.compute(ReflectivityOverQ)\n", + "covered = (~reflectivity.masks['direct_beam']).sum().value\n", + "print(f'{covered} of {reflectivity.sizes[\"Q\"]} Q bins have direct-beam coverage.')\n", + "title = (\n", + " 'Reflectivity'\n", + " if beam_size is not None\n", + " else 'Reflectivity (footprint correction omitted)'\n", + ")\n", + "reflectivity.plot(norm='log', title=title)" + ] + }, + { + "cell_type": "markdown", + "id": "15", + "metadata": {}, + "source": [ + "The sparse statistics and omitted footprint correction limit the interpretation of this curve. Q resolution, background subtraction, finite-sample corrections, and complete ORSO export are not yet included." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "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.12.13" + }, + "nbsphinx": { + "execute": "never" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/packages/essreflectometry/docs/user-guide/freia/freia-wavelength-lookup-table.ipynb b/packages/essreflectometry/docs/user-guide/freia/freia-wavelength-lookup-table.ipynb index 7c4eb9c9a..133005254 100644 --- a/packages/essreflectometry/docs/user-guide/freia/freia-wavelength-lookup-table.ipynb +++ b/packages/essreflectometry/docs/user-guide/freia/freia-wavelength-lookup-table.ipynb @@ -7,9 +7,7 @@ "source": [ "# FREIA analytical wavelength lookup table\n", "\n", - "Explore which neutron wavelengths can reach the detector at each arrival time for the WFM chopper configuration.\n", - "\n", - "The calculation uses the final detector's pixel geometry and the source and sample positions. Detector events and histogram intensities are not needed to construct the lookup table.\n" + "Explore which wavelengths can reach the detector at each arrival time for the WFM chopper configuration. This calculation needs geometry and chopper settings; detector events are not required.\n" ] }, { @@ -44,9 +42,7 @@ "source": [ "## Use the WFM configuration\n", "\n", - "Use the fixed **WFM** chopper settings and read detector geometry from the sample simulation.\n", - "\n", - "This example uses 20 µs time resolution and 10 cm flight-path resolution. The central-ray model projects chopper positions onto the global z axis and approximates detector flight paths by source-to-sample plus sample-to-pixel distance.\n" + "Set the input file and lookup-table resolution. This example uses 20 µs time resolution and 10 cm flight-path resolution.\n" ] }, { @@ -95,9 +91,9 @@ "source": [ "## Wavelength as a function of arrival time and flight path\n", "\n", - "The heatmap shows the wavelength assigned to each arrival time and flight path. White regions have no modeled transmission. The lower panel shows a slice near the middle of the detector's flight-path range; shading spans the modeled lower and upper wavelength bounds. The table stores the squared half-width of this range as its variance; it describes wavelength ambiguity, not counting statistics.\n", + "The heatmap shows the wavelength assigned to each arrival time and flight path. White regions have no modeled transmission. The lower panel shows a slice near the middle of the detector's flight-path range; shading spans the modeled wavelength bounds.\n", "\n", - "This is the unmasked lookup table: no additional relative-uncertainty cutoff has been applied. The lookup grid extends slightly beyond the detector's flight-path range to support interpolation.\n" + "No additional relative-uncertainty cutoff has been applied.\n" ] }, { diff --git a/packages/essreflectometry/docs/user-guide/freia/index.md b/packages/essreflectometry/docs/user-guide/freia/index.md index f2b7dbb96..4844732c8 100644 --- a/packages/essreflectometry/docs/user-guide/freia/index.md +++ b/packages/essreflectometry/docs/user-guide/freia/index.md @@ -1,13 +1,12 @@ # FREIA -Explore FREIA simulations through detector images, arrival-time distributions, and -wavelength spectra. The guides use fixed **WFM** chopper settings to reconstruct -wavelengths. The wavelength lookup-table guide requires only detector geometry and -the source and sample positions; detector events are not needed. +Explore FREIA detector data, reconstruct wavelengths, and reduce reflectivity +using a direct-beam measurement. The examples use WFM chopper settings. ```{toctree} :maxdepth: 1 freia-mcstas-visualization freia-wavelength-lookup-table +freia-reflectivity ``` diff --git a/packages/essreflectometry/src/ess/freia/beamline.py b/packages/essreflectometry/src/ess/freia/beamline.py index e6ef59fa1..fc87b9c4e 100644 --- a/packages/essreflectometry/src/ess/freia/beamline.py +++ b/packages/essreflectometry/src/ess/freia/beamline.py @@ -4,6 +4,8 @@ from ess.reflectometry.types import ReferenceRun, SampleRun +from .types import IncidentMonitor + DETECTOR_BANK_SIZES = { "multiblade_detector": { "strip": 64, @@ -19,7 +21,7 @@ def LoadNeXusWorkflow(**kwargs) -> sciline.Pipeline: """ workflow = GenericUnwrapWorkflow( run_types=[SampleRun, ReferenceRun], - monitor_types=[], + monitor_types=[IncidentMonitor], **kwargs, ) workflow[DetectorBankSizes] = DETECTOR_BANK_SIZES diff --git a/packages/essreflectometry/src/ess/freia/conversions.py b/packages/essreflectometry/src/ess/freia/conversions.py index 35ef06657..036b277a7 100644 --- a/packages/essreflectometry/src/ess/freia/conversions.py +++ b/packages/essreflectometry/src/ess/freia/conversions.py @@ -1,171 +1,98 @@ # SPDX-License-Identifier: BSD-3-Clause -# Copyright (c) 2025 Scipp contributors (https://github.com/scipp) +# Copyright (c) 2026 Scipp contributors (https://github.com/scipp) import scipp as sc -from ess.reduce.nexus.types import DetectorBankSizes, Position -from scipp.constants import pi -from scippneutron._utils import elem_dtype -from scippneutron.conversion import graph +from ess.reduce.nexus.types import GravityVector, Position +from scippneutron.conversion import beamline, graph from scippnexus import NXsample, NXsource from ..reflectometry.conversions import reflectometry_q -from ..reflectometry.types import ( - CoordTransformationGraph, - DetectorRotation, - RunType, - SampleRotation, -) +from ..reflectometry.types import CoordTransformationGraph, RunType, WavelengthDetector +from .types import QDetector, SampleSurfaceNormal -def reflectometry_q_x( - wavelength: sc.Variable, theta: sc.Variable, sample_rotation: sc.Variable +def theta( + incident_beam: sc.Variable, + scattered_beam: sc.Variable, + wavelength: sc.Variable, + gravity: sc.Variable, + sample_surface_normal: sc.Variable, ) -> sc.Variable: - """ - Compute momentum transfer in off-specular direction. - - .. math:: - Q_x = \\frac{2 \\pi}{\\lambda} (cos(\\theta_o) - cos(\\theta_i)) - - Where :math:`\\theta_o` is the reflection angle (:func:`theta`) - and :math:`\\theta_i` is the incident angle on the sample surface. - - Note that here we assume the incident angle is equal to ``sample_rotation``. + """Signed, gravity-corrected exit angle above the sample plane. - Source: - `Frédéric Ott, "Off-specular data representations in neutron reflectivity" `_ + ScippNeutron reconstructs the outgoing direction at the sample. Project + that direction onto the sample normal to retain the sign and support a + tilted sample. Its reflectometry-specific scattering_angle_in_yz_plane + returns an unsigned angle, which cannot distinguish the direct beam. - Parameters - ---------- - wavelength: - Wavelength values for the events. - theta: - Angle of reflection for the events. - sample_rotation: - Angle of incidence. - - Returns - ------- - : - Qx-values. + The horizontal beam direction only defines a coordinate basis here; it + does not specify an incident angle. For Q and footprint we still assume + specular reflection, so the incidence angle equals this exit angle. """ - dtype = elem_dtype(wavelength) - c = (2 * pi).astype(dtype) - return ( - c - * ( - sc.cos(theta.astype(dtype, copy=False)) - - sc.cos(sample_rotation.to(unit=theta.unit, dtype=dtype)) - ) - / wavelength - ) - - -def theta( - divergence_angle: sc.Variable, - sample_rotation: sc.Variable, -): - ''' - Angle of reflection. - - Computes the angle between the scattering direction of - the neutron and the sample surface. - - Parameters - ------------ - divergence_angle: - Divergence angle of the scattered beam. - sample_rotation: - Rotation of the sample from to its zero position. - - Returns - ----------- - The reflection angle of the neutron. - ''' - return divergence_angle + sample_rotation.to( - unit=divergence_angle.unit, dtype='float64' + basis = beamline.beam_aligned_unit_vectors(incident_beam, gravity) + x, y, z = (basis[f'beam_aligned_unit_{axis}'] for axis in 'xyz') + # Use the horizontal reference axis: the source-to-sample line in FREIA + # is tilted and does not describe the incident direction at the sample. + angles = beamline.scattering_angles_with_gravity( + incident_beam=z * sc.scalar(1.0, unit='m'), + scattered_beam=scattered_beam, + wavelength=wavelength, + gravity=gravity, ) - - -def divergence_angle( - position: sc.Variable, - sample_position: sc.Variable, - detector_rotation: sc.Variable, -): - """ - Angle between the scattering ray and - the ray that travels parallel to the sample surface - when the sample rotation is zero. - - Parameters - ------------ - position: - Detector position where the neutron was detected. - sample_position: - Position of the sample. - detector_rotation: - Rotation of the detector from its zero position. - Returns - ---------- - The divergence angle of the scattered beam. - """ - p = position - sample_position.to(unit=position.unit) - return sc.atan2(y=p.fields.x, x=p.fields.z) - detector_rotation.to( - unit='rad', dtype='float64' + polar, azimuth = angles['two_theta'], angles['phi'] + normal = sample_surface_normal / sc.norm(sample_surface_normal) + return sc.asin( + sc.sin(polar) + * (sc.dot(normal, x) * sc.cos(azimuth) + sc.dot(normal, y) * sc.sin(azimuth)) + + sc.dot(normal, z) * sc.cos(polar) ) def coordinate_transformation_graph( source_position: Position[NXsource, RunType], sample_position: Position[NXsample, RunType], - sample_rotation: SampleRotation[RunType], - detector_rotation: DetectorRotation[RunType], - detector_bank_sizes: DetectorBankSizes, + sample_surface_normal: SampleSurfaceNormal[RunType], + gravity: GravityVector, ) -> CoordTransformationGraph[RunType]: - bank = detector_bank_sizes['multiblade_detector'] + """Build a specular conversion graph independent of detector pixel layout.""" + length = sc.norm(sample_surface_normal) + if ( + not sc.isfinite(length).value + or not (length > sc.scalar(0.0, unit=length.unit)).value + ): + raise ValueError('SampleSurfaceNormal must be a finite, nonzero vector.') return { **graph.beamline.beamline(scatter=True), - "theta": theta, - "divergence_angle": divergence_angle, - "Q": reflectometry_q, - "Qx": reflectometry_q_x, - 'sample_size': lambda: sc.scalar(20.0, unit='mm'), - 'blade': lambda: sc.arange('blade', bank['blade'] - 1, -1, -1), - 'wire': lambda: sc.arange('wire', bank['wire'] - 1, -1, -1), - 'strip': lambda: sc.arange('strip', bank['strip'] - 1, -1, -1), - 'z_index': lambda blade, wire: blade * wire, - 'sample_rotation': lambda: sample_rotation, - 'detector_rotation': lambda: detector_rotation, + 'theta': theta, + 'Q': reflectometry_q, 'source_position': lambda: source_position, 'sample_position': lambda: sample_position, + 'sample_surface_normal': lambda: sample_surface_normal, + 'gravity': lambda: gravity, } def add_coords( - da: sc.DataArray, - graph: dict, -) -> sc.DataArray: - "Adds scattering coordinates to the raw detector data." - return da.transform_coords( - ( - "wavelength", - "theta", - "divergence_angle", - "Q", - "Qx", - "L1", - "L2", - "blade", - "wire", - "strip", - "z_index", - "sample_rotation", - "detector_rotation", - "sample_size", - ), - graph, - rename_dims=False, - keep_intermediate=False, - keep_aliases=False, + da: WavelengthDetector[RunType], + graph: CoordTransformationGraph[RunType], +) -> QDetector[RunType]: + """Add Q without requiring an ROI, monitor, reference, or footprint inputs.""" + return QDetector[RunType]( + da.transform_coords( + ( + 'theta', + 'Q', + 'L1', + 'L2', + 'incident_beam', + 'sample_position', + 'sample_surface_normal', + ), + graph, + rename_dims=False, + keep_intermediate=False, + keep_aliases=False, + ) ) -providers = (coordinate_transformation_graph,) +providers = (coordinate_transformation_graph, add_coords) diff --git a/packages/essreflectometry/src/ess/freia/corrections.py b/packages/essreflectometry/src/ess/freia/corrections.py index a6f5c1238..d333589f1 100644 --- a/packages/essreflectometry/src/ess/freia/corrections.py +++ b/packages/essreflectometry/src/ess/freia/corrections.py @@ -7,20 +7,15 @@ from ..reflectometry import corrections as common_corrections from ..reflectometry.corrections import RunNormalization from ..reflectometry.types import ( - BeamDivergenceLimits, - CoordTransformationGraph, - CorrectionsToApply, + BeamSize, ReducibleData, RunType, RunUnnormalizedData, - WavelengthBins, - WavelengthDetector, - YIndexLimits, - ZIndexLimits, + Sample, + SampleRun, + SampleSize, ) -from .conversions import add_coords -from .maskings import add_masks -from .types import WavelengthMonitor +from .types import SampleIlluminatedFraction, WavelengthMonitor def normalize_by_monitor_histogram( @@ -89,33 +84,45 @@ def insert_run_normalization( ) -def add_coords_masks_and_apply_corrections( - da: WavelengthDetector[RunType], - ylim: YIndexLimits, - zlims: ZIndexLimits, - bdlim: BeamDivergenceLimits, - wbins: WavelengthBins, - graph: CoordTransformationGraph[RunType], - corrections_to_apply: CorrectionsToApply, -) -> RunUnnormalizedData[RunType]: - """ - Computes coordinates, masks and corrections that are - the same for the sample measurement and the reference measurement. - """ - da = add_coords(da, graph) - da = add_masks(da, ylim, zlims, bdlim, wbins) - - for correction in corrections_to_apply: - da = correction(da) - - return RunUnnormalizedData[RunType](da) +def sample_illuminated_fraction( + sample: ReducibleData[SampleRun], + beam_size: BeamSize[SampleRun], + sample_size: SampleSize[SampleRun], +) -> SampleIlluminatedFraction: + """Use Amor's Gaussian footprint model with the specular incidence angle. + Beam size is the FWHM at the sample; sample size is its length along the + beam. These must be supplied explicitly, since slit openings alone do not + determine the profile of the beam reaching the sample. + """ + for name, size in [('BeamSize', beam_size), ('SampleSize', sample_size)]: + if ( + not sc.isfinite(size).value + or not (size > sc.scalar(0.0, unit=size.unit)).value + ): + raise ValueError(f'{name} must be finite and positive.') + return SampleIlluminatedFraction( + common_corrections.footprint_on_sample( + sample.bins.coords['theta'], beam_size=beam_size, sample_size=sample_size + ) + ) -def correct_by_footprint(da: sc.DataArray) -> sc.DataArray: - """Corrects the data by the size of the footprint on the sample.""" - return da / sc.sin(da.coords['theta']) +def prepare_sample( + sample: ReducibleData[SampleRun], + illuminated_fraction: SampleIlluminatedFraction, +) -> Sample: + """Correct the reflected beam for footprint; the direct beam has no sample.""" + if illuminated_fraction.bins is None: + illuminated_fraction = sc.bins_like(sample, illuminated_fraction) + valid = sc.isfinite(illuminated_fraction) & (illuminated_fraction > sc.scalar(0.0)) + valid &= illuminated_fraction <= sc.scalar(1.0) + sample = sample.bins.assign_masks( + footprint=~valid, + non_reflected=sample.bins.coords['theta'] <= sc.scalar(0.0, unit='rad'), + ) + fraction = sc.where(valid, illuminated_fraction, sc.scalar(1.0)) + return Sample(sample / fraction) -default_corrections = {correct_by_footprint} -providers = (add_coords_masks_and_apply_corrections,) +providers = (sample_illuminated_fraction, prepare_sample) diff --git a/packages/essreflectometry/src/ess/freia/maskings.py b/packages/essreflectometry/src/ess/freia/maskings.py index 472e8bc93..8dea32f8d 100644 --- a/packages/essreflectometry/src/ess/freia/maskings.py +++ b/packages/essreflectometry/src/ess/freia/maskings.py @@ -1,54 +1,39 @@ # SPDX-License-Identifier: BSD-3-Clause -# Copyright (c) 2025 Scipp contributors (https://github.com/scipp) +# Copyright (c) 2026 Scipp contributors (https://github.com/scipp) import scipp as sc -from ..reflectometry.types import ( - BeamDivergenceLimits, - WavelengthBins, - YIndexLimits, - ZIndexLimits, -) +from ..reflectometry.types import RunType, RunUnnormalizedData, WavelengthBins +from .types import DetectorRegionOfInterest, QDetector -def _not_between(v, a, b): - return (v < a) | (v > b) - - -def add_masks( - da: sc.DataArray, - ylim: YIndexLimits, - zlims: ZIndexLimits, - bdlim: BeamDivergenceLimits, - wbins: WavelengthBins, -) -> sc.DataArray: - """ - Masks the data by ranges in the detector - coordinates ``z`` and ``y``, and by the divergence of the beam, - and by wavelength. - """ - da = da.assign_masks( - strip_range=_not_between(da.coords["strip"], *ylim), - z_range=_not_between(da.coords["z_index"], *zlims), - divergence_too_large=_not_between( - da.coords["divergence_angle"], - bdlim[0].to( - unit=da.coords["divergence_angle"].unit, - dtype='float64', +def select_events( + da: QDetector[RunType], + roi: DetectorRegionOfInterest[RunType], + wavelength_bins: WavelengthBins, +) -> RunUnnormalizedData[RunType]: + """Select a peak using pixel or event coordinates and wavelength.""" + masks = {} + event_masks = {} + for name, (low, high) in roi.items(): + is_event_coord = name in da.bins.coords + coord = da.bins.coords[name] if is_event_coord else da.coords[name] + low, high = low.to(unit=coord.unit), high.to(unit=coord.unit) + if not (low <= high).value: + raise ValueError(f'Reversed ROI bounds for {name!r}.') + (event_masks if is_event_coord else masks)[f'roi_{name}'] = ~( + (coord >= low) & (coord <= high) + ) + wavelength = da.bins.coords['wavelength'] + return RunUnnormalizedData[RunType]( + da.assign_masks(masks).bins.assign_masks( + event_masks, + wavelength=~( + sc.isfinite(wavelength) + & (wavelength >= wavelength_bins[0].to(unit=wavelength.unit)) + & (wavelength < wavelength_bins[-1].to(unit=wavelength.unit)) ), - bdlim[1].to( - unit=da.coords["divergence_angle"].unit, - dtype='float64', - ), - ), - ) - da = da.bins.assign_masks( - wavelength=_not_between( - da.bins.coords['wavelength'], - wbins[0], - wbins[-1], - ), + ) ) - return da -providers = () +providers = (select_events,) diff --git a/packages/essreflectometry/src/ess/freia/mcstas.py b/packages/essreflectometry/src/ess/freia/mcstas.py index edb8701c2..c35421d06 100644 --- a/packages/essreflectometry/src/ess/freia/mcstas.py +++ b/packages/essreflectometry/src/ess/freia/mcstas.py @@ -19,6 +19,7 @@ EmptyDetector, Filename, NeXusDetectorName, + NeXusName, Position, RawDetector, RunType, @@ -26,6 +27,8 @@ from ess.reduce.unwrap import PulsePeriod from scippneutron.chopper import DiskChopper +from .types import IncidentMonitor, SampleSurfaceNormal, WavelengthMonitor + class _ChopperParameters(TypedDict): frequency: float @@ -222,6 +225,7 @@ def load_mcstas( Empty detector pixels are retained, and weighted-event variances are ``p**2``. Arrival times are split into ``event_time_zero`` and ``event_time_offset`` using ``pulse_period``, which defaults to the ESS period of 1/14 s. + If present, ``L`` is kept as ``wavelength_from_mcstas`` in angstroms. """ if pulse_period is None: pulse_period = sc.scalar(1 / 14, unit='s') @@ -240,8 +244,11 @@ def _load_events( f'No Mantid detector events with pixel IDs found for {detector_name!r}. ' 'Enable "mantid banana ... list all neutrons" in the simulation.' ) + variables = ['p', 't', 'id'] + if 'L' in data.get_component_variables(detector_name): + variables.append('L') values = data.get_event_data( - variables=['p', 't', 'id'], component_name=detector_name, filter_zeros=True + variables=variables, component_name=detector_name, filter_zeros=True ) if geometry is None: geometry = _detector_geometry(data, detector_name) @@ -265,8 +272,12 @@ def _load_events( 'event_time_zero': sc.datetime(0, unit='ns') + (time - offset).to(unit='ns', dtype='int64'), }, - ).group(pixel_ids) - return events.assign_coords(geometry.coords) + ) + if 'L' in values: + events.coords['wavelength_from_mcstas'] = sc.array( + dims=['event'], values=values['L'], unit='angstrom' + ) + return events.group(pixel_ids).assign_coords(geometry.coords) def load_mcstas_provider( @@ -301,6 +312,53 @@ def mcstas_sample_position( return Position[snx.NXsample, RunType](_component_position(filename, 'Arm_Sample')) +def mcstas_sample_surface_normal( + filename: Filename[RunType], +) -> SampleSurfaceNormal[RunType]: + """Load the sample surface normal in global coordinates.""" + with _open_mcstas(filename) as data: + _, rotation = data.get_component_placement('Arm_Sample') + # mcstastox transforms row vectors from local to global with local @ rotation. + return SampleSurfaceNormal[RunType](sc.vector(rotation[1], unit='dimensionless')) + + +def load_mcstas_monitor(filename: str | Path, monitor_name: str) -> sc.DataArray: + """Read a selected one-dimensional McStas L_monitor histogram. + + Select an incident monitor upstream of the sample. In particular, + SampleLambda in FREIA_surface_test.instr is downstream of the sample. + McStas stores bin-integrated intensities and standard errors; xlimits + supplies the bounds of the uniformly spaced wavelength bins. + """ + with _open_mcstas(filename) as data: + histogram = data.file_object.get_info_entry(monitor_name) + if _text(histogram.attrs.get('xvar', '')) != 'L' or histogram['data'].ndim != 1: + raise ValueError(f'{monitor_name!r} is not a 1-D wavelength monitor.') + low, high = map(float, _text(histogram.attrs['xlimits']).split()) + intensity = histogram['data'][:] + return sc.DataArray( + sc.array( + dims=['wavelength'], + values=intensity, + variances=histogram['errors'][:] ** 2, + unit='counts', + ), + coords={ + 'wavelength': sc.linspace( + 'wavelength', low, high, len(intensity) + 1, unit='angstrom' + ), + }, + ) + + +def mcstas_wavelength_monitor( + filename: Filename[RunType], + monitor_name: NeXusName[IncidentMonitor], +) -> WavelengthMonitor[RunType]: + """Provide the explicitly selected incident wavelength histogram.""" + return WavelengthMonitor[RunType](load_mcstas_monitor(filename, monitor_name)) + + def _histogram_axis(histogram, axis, unit): label = _text(histogram.attrs[f'{axis}label']) name = re.sub(r'[^a-zA-Z]', '_', label) @@ -376,5 +434,7 @@ def mcstas_detector_geometry( load_mcstas_provider, mcstas_source_position, mcstas_sample_position, + mcstas_sample_surface_normal, + mcstas_wavelength_monitor, mcstas_detector_geometry, ) diff --git a/packages/essreflectometry/src/ess/freia/normalization.py b/packages/essreflectometry/src/ess/freia/normalization.py index 4c48502ca..39f8ce41f 100644 --- a/packages/essreflectometry/src/ess/freia/normalization.py +++ b/packages/essreflectometry/src/ess/freia/normalization.py @@ -1,3 +1,70 @@ -# Copyright (c) 2025 Scipp contributors (https://github.com/scipp) +# SPDX-License-Identifier: BSD-3-Clause +# Copyright (c) 2026 Scipp contributors (https://github.com/scipp) +import scipp as sc +from ess.reduce.nexus.types import GravityVector -providers = () +from ..reflectometry.conversions import reflectometry_q +from ..reflectometry.types import ( + QBins, + ReducibleData, + Reference, + ReferenceRun, + ReflectivityOverQ, + Sample, + SampleRun, +) +from .conversions import theta +from .types import SampleSurfaceNormal + + +def evaluate_direct_beam( + direct_beam: ReducibleData[ReferenceRun], + sample_surface_normal: SampleSurfaceNormal[SampleRun], + gravity: GravityVector, +) -> Reference: + """Map the direct-beam ROI to Q after reflection at the sample surface. + + A direct ray follows the incident direction. Specular reflection reverses + its normal component, so its exit angle is the negative of its angle above + the sample plane. Using the sample run's normal also supports a changed + sample orientation. No footprint or supermirror correction is applied. + """ + angle = -theta( + incident_beam=direct_beam.coords['incident_beam'], + scattered_beam=direct_beam.coords['position'] + - direct_beam.coords['sample_position'], + wavelength=direct_beam.bins.coords['wavelength'], + gravity=gravity, + sample_surface_normal=sample_surface_normal, + ) + reference = direct_beam.bins.assign_coords( + theta=angle, Q=reflectometry_q(direct_beam.bins.coords['wavelength'], angle) + ) + return Reference( + reference.bins.assign_masks(non_incident=angle <= sc.scalar(0.0, unit='rad')) + ) + + +def reduce_sample_over_q( + sample: Sample, + reference: Reference, + qbins: QBins, +) -> ReflectivityOverQ: + """Divide ROI intensities on a common Q grid, propagating both variances. + + Histogram before division so changing Q bin widths does not rescale R. + Empty, masked, or nonfinite direct-beam bins provide no normalization. + """ + numerator = sample.hist(Q=qbins, dim=sample.dims) + denominator = reference.hist(Q=qbins, dim=reference.dims) + valid = sc.isfinite(denominator.data) & ( + denominator.data > sc.scalar(0.0, unit=denominator.unit) + ) + # Keep unusable bins explicit instead of producing infinities at empty bins. + norm = sc.where( + valid, denominator.data, sc.scalar(float('nan'), unit=denominator.unit) + ) + return ReflectivityOverQ((numerator / norm).assign_masks(direct_beam=~valid)) + + +providers = (evaluate_direct_beam, reduce_sample_over_q) diff --git a/packages/essreflectometry/src/ess/freia/types.py b/packages/essreflectometry/src/ess/freia/types.py index 8b410217e..3b36b3b53 100644 --- a/packages/essreflectometry/src/ess/freia/types.py +++ b/packages/essreflectometry/src/ess/freia/types.py @@ -1,6 +1,7 @@ # Copyright (c) 2025 Scipp contributors (https://github.com/scipp) from typing import NewType +import sciline import scipp as sc from ess.reduce.unwrap.types import WavelengthMonitor as _WavelengthMonitor @@ -16,3 +17,24 @@ # generic alias at runtime, and it is subscripted again as # ``WavelengthMonitor[RunType]`` in providers. WavelengthMonitor = _WavelengthMonitor[RunType, IncidentMonitor] + + +class SampleSurfaceNormal(sciline.Scope[RunType, sc.Variable], sc.Variable): + """Normal pointing out of the reflecting surface, in global coordinates.""" + + +class QDetector(sciline.Scope[RunType, sc.DataArray], sc.DataArray): + """Detector events with specular Q and signed angle above the sample surface.""" + + +class DetectorRegionOfInterest(sciline.Scope[RunType, dict], dict): + """Pixel or event coordinates mapped to inclusive (lower, upper) bounds. + + Select corresponding reflected and direct peaks separately for SampleRun and + ReferenceRun, for example using ``theta`` and ``height``. + An empty dictionary explicitly selects the entire detector. + """ + + +SampleIlluminatedFraction = NewType('SampleIlluminatedFraction', sc.Variable) +"""Fraction of the incoming beam hitting the sample; set to 1 to skip footprint.""" diff --git a/packages/essreflectometry/src/ess/freia/workflow.py b/packages/essreflectometry/src/ess/freia/workflow.py index fae23f746..0af66d670 100644 --- a/packages/essreflectometry/src/ess/freia/workflow.py +++ b/packages/essreflectometry/src/ess/freia/workflow.py @@ -9,8 +9,6 @@ from ..reflectometry import providers as reflectometry_providers from ..reflectometry.types import ( - BeamDivergenceLimits, - CorrectionsToApply, DetectorSpatialResolution, LookupTableRelativeErrorThreshold, NeXusDetectorName, @@ -49,10 +47,6 @@ def mcstas_default_parameters() -> dict: """Return default parameters for the McStas Freia workflow.""" return default_parameters() | { NeXusDetectorName: "Multiblade", - BeamDivergenceLimits: ( - sc.scalar(-0.75, unit='deg'), - sc.scalar(0.75, unit='deg'), - ), LookupTableRelativeErrorThreshold: { "Multiblade": 0.06, }, @@ -64,7 +58,6 @@ def default_parameters() -> dict: return { NeXusDetectorName: "multiblade_detector", SampleRotationOffset[RunType]: sc.scalar(0.0, unit='deg'), - CorrectionsToApply: corrections.default_corrections, DetectorSpatialResolution: 0.0025 * sc.units.m, LookupTableRelativeErrorThreshold: { "multiblade_detector": float('inf'), @@ -79,14 +72,10 @@ def FreiaMcStasWorkflow( wavelength_from: WavelengthLutMode = "analytical", **kwargs, ) -> sciline.Pipeline: - """Workflow for loading and unwrapping FREIA McStas detector events. + """Workflow for reducing FREIA McStas events with a no-sample direct beam. - Builds on the generic NeXus/unwrapping workflow, replacing input providers - with adapters for the final ``Multiblade`` Mantid banana detector. Detector - geometry is read from the simulation file; choppers use the fixed WFM - configuration from :func:`ess.freia.mcstas.wfm_choppers`. This initial - workflow supports visualization through ``WavelengthDetector``; full - reflectivity reduction and normalization require further instrument providers. + Loads geometry and uses the default WFM chopper settings. Reduction inputs + and outputs are described in :func:`FreiaWorkflow`. Parameters ---------- @@ -115,6 +104,14 @@ def FreiaWorkflow( ) -> sciline.Pipeline: """Workflow for reduction of data for the Freia instrument. + ``QDetector`` provides specular Q with gravity correction. Reflectivity + requires separate sample/direct-beam ROIs, wavelength and Q bins, and beam + and sample sizes for the footprint correction. The reference run must be a + measurement without a sample, taken with matching slit and chopper settings. + + Monitor normalization requires an incident monitor selected through + ``NeXusName[IncidentMonitor]``, or supplied as ``WavelengthMonitor[RunType]``. + Parameters ---------- run_norm: diff --git a/packages/essreflectometry/tests/freia/mcstas_test.py b/packages/essreflectometry/tests/freia/mcstas_test.py index d7c2d01ec..c8754c063 100644 --- a/packages/essreflectometry/tests/freia/mcstas_test.py +++ b/packages/essreflectometry/tests/freia/mcstas_test.py @@ -11,6 +11,7 @@ DiskChoppers, Filename, NeXusDetectorName, + NeXusName, RawDetector, ) from ess.reduce.unwrap import ( @@ -27,6 +28,7 @@ from ess.freia import FreiaMcStasWorkflow from ess.freia.mcstas import load_mcstas +from ess.freia.types import IncidentMonitor, QDetector, WavelengthMonitor from ess.reflectometry.types import ReferenceRun, SampleRun @@ -188,10 +190,13 @@ def test_workflow_loads_and_unwraps_detector_with_generic_providers(mcstas_file, time=(sc.scalar(0.9, unit='ms'), sc.scalar(1.1, unit='ms')), wavelength=(sc.scalar(0.5, unit='angstrom'), sc.scalar(12.0, unit='angstrom')), ) - result = workflow.compute((RawDetector[run], WavelengthDetector[run])) + result = workflow.compute( + (RawDetector[run], WavelengthDetector[run], QDetector[run]) + ) raw = result[RawDetector[run]] unwrapped = result[WavelengthDetector[run]] assert_allclose(raw.bins.sum().data, unwrapped.bins.sum().data) + assert 'Q' in result[QDetector[run]].bins.coords wavelength = unwrapped.bins.constituents['data'].coords['wavelength'] # Arrival time 25 ms minus emission time 1 ms, flight path about 23 m. expected = ( @@ -224,3 +229,34 @@ def test_loader_uses_workflow_pulse_period(mcstas_file): dims=['event'], values=[0.025, 0.025, 0.025 + 1 / 14 - 0.05], unit='s' ), ) + + +def test_load_selected_wavelength_monitor_with_bin_edges_and_variances(mcstas_file): + with h5py.File(mcstas_file, 'r+') as f: + components = f['entry1/instrument/components'] + monitor = _component(components, 'IncidentLambda', [0.0, 0.0, 19.0]) + histogram = monitor.create_group('output/spectrum') + histogram.attrs['xvar'] = np.bytes_('L') + histogram.attrs['xlimits'] = np.bytes_('1 5') + histogram['data'] = [10.0, 20.0] + histogram['errors'] = [3.0, 4.0] + workflow = FreiaMcStasWorkflow() + workflow[Filename[SampleRun]] = str(mcstas_file) + workflow[NeXusName[IncidentMonitor]] = 'IncidentLambda' + result = workflow.compute(WavelengthMonitor[SampleRun]) + assert_identical( + result, + sc.DataArray( + sc.array( + dims=['wavelength'], + values=[10.0, 20.0], + variances=[9.0, 16.0], + unit='counts', + ), + coords={ + 'wavelength': sc.array( + dims=['wavelength'], values=[1.0, 3.0, 5.0], unit='angstrom' + ) + }, + ), + ) diff --git a/packages/essreflectometry/tests/freia/workflow_test.py b/packages/essreflectometry/tests/freia/workflow_test.py index fbf869f16..702e2ce9a 100644 --- a/packages/essreflectometry/tests/freia/workflow_test.py +++ b/packages/essreflectometry/tests/freia/workflow_test.py @@ -1,8 +1,35 @@ # SPDX-License-Identifier: BSD-3-Clause # Copyright (c) 2026 Scipp contributors (https://github.com/scipp) +# Small known intensities exercise the complete reduction graph without a LUT. +import numpy as np +import pytest +import scipp as sc from ess.reduce import workflow as reduce_workflow +from ess.reduce.nexus.types import Position +from scipp.testing import assert_allclose, assert_identical +from scippnexus import NXsample, NXsource from ess import freia +from ess.freia.corrections import RunNormalization +from ess.freia.types import ( + DetectorRegionOfInterest, + QDetector, + SampleIlluminatedFraction, + SampleSurfaceNormal, + WavelengthMonitor, +) +from ess.reflectometry.types import ( + BeamSize, + QBins, + Reference, + ReferenceRun, + ReflectivityOverQ, + RunUnnormalizedData, + SampleRun, + SampleSize, + WavelengthBins, + WavelengthDetector, +) def test_freia_workflow_registers_run_normalization_variants(): @@ -17,3 +44,176 @@ def test_freia_workflow_registers_run_normalization_variants(): freia.FreiaProtonChargeWorkflow, ): assert wf in reduce_workflow.workflow_registry + + +def _detector(angles, weights): + position = np.array([1.0, 2.0, 3.0]) + 3.0 * np.column_stack( + [np.zeros(len(angles)), np.sin(np.deg2rad(angles)), np.cos(np.deg2rad(angles))] + ) + return ( + sc.DataArray( + sc.array(dims=['event'], values=weights, variances=weights, unit='counts'), + coords={ + 'wavelength': sc.array( + dims=['event'], values=[2.0, 4.0, 2.0], unit='angstrom' + ), + 'pixel_id': sc.array(dims=['event'], values=[0, 0, 1], unit=None), + }, + ) + .group('pixel_id') + .assign_coords( + position=sc.vectors(dims=['pixel_id'], values=position[[0, 2]], unit='m') + ) + ) + + +def _workflow(run_norm=RunNormalization.none): + wf = freia.FreiaWorkflow(run_norm=run_norm) + for run, angles, weights in [ + (SampleRun, [1.0, 1.0, 2.0], [20.0, 40.0, 999.0]), + (ReferenceRun, [-1.0, -1.0, -2.0], [100.0, 100.0, 999.0]), + ]: + wf[WavelengthDetector[run]] = _detector(angles, weights) + wf[Position[NXsample, run]] = sc.vector([1.0, 2.0, 3.0], unit='m') + wf[Position[NXsource, run]] = sc.vector([1.0, 2.0, -17.0], unit='m') + wf[SampleSurfaceNormal[run]] = sc.vector([0.0, 1.0, 0.0]) + wf[DetectorRegionOfInterest[SampleRun]] = { + 'theta': (sc.scalar(0.5, unit='deg'), sc.scalar(1.5, unit='deg')) + } + wf[DetectorRegionOfInterest[ReferenceRun]] = { + 'theta': (sc.scalar(-1.5, unit='deg'), sc.scalar(-0.5, unit='deg')) + } + wf[WavelengthBins] = sc.array( + dims=['wavelength'], values=[1.0, 3.0, 5.0], unit='angstrom' + ) + wf[QBins] = sc.array(dims=['Q'], values=[0.04, 0.08, 0.13, 0.3], unit='1/angstrom') + return wf + + +def test_q_uses_sample_plane_and_translated_pixels_without_reduction_inputs(): + wf = _workflow() + result = wf.compute(QDetector[SampleRun]) + q = result.bins.constituents['data'].coords['Q'] + expected = 4 * np.pi * np.sin(_expected_exit_angles()) / [2.0, 4.0, 2.0] + assert_allclose(q, sc.array(dims=['event'], values=expected, unit='1/angstrom')) + assert_allclose( + result.coords['L2'], sc.full(sizes=result.sizes, value=3.0, unit='m') + ) + + +def test_q_follows_tilted_sample_surface(): + wf = _workflow() + tilt = np.deg2rad(0.5) + wf[SampleSurfaceNormal[SampleRun]] = sc.vector([0.0, np.cos(tilt), -np.sin(tilt)]) + result = wf.compute(QDetector[SampleRun]) + assert_allclose( + result.bins.constituents['data'].coords['theta'], + sc.array(dims=['event'], values=_expected_exit_angles() - tilt, unit='rad'), + ) + + +def _expected_exit_angles(sign=1.0): + """Independent calculation from displacement = velocity*time + gravity*time²/2.""" + angle = sign * np.deg2rad([1.0, 1.0, 2.0]) + speed = ( + ( + sc.constants.h + / sc.constants.m_n + / sc.array(dims=['event'], values=[2.0, 4.0, 2.0], unit='angstrom') + ) + .to(unit='m/s') + .values + ) + time = 3.0 / speed + y = 3.0 * np.sin(angle) + 0.5 * sc.constants.g.value * time**2 + return np.arctan2(y, 3.0 * np.cos(angle)) + + +def test_theta_roi_distinguishes_wavelengths_in_the_same_pixel(): + wf = _workflow() + wf[DetectorRegionOfInterest[SampleRun]] = { + 'theta': (sc.scalar(1.0005, unit='deg'), sc.scalar(1.001, unit='deg')), + } + selected = wf.compute(RunUnnormalizedData[SampleRun]) + assert selected.bins.sum().sum().value == 40.0 + + +def test_direct_beam_gravity_is_corrected_before_specular_mapping(): + wf = _workflow() + reference = wf.compute(Reference) + assert_allclose( + reference.bins.constituents['data'].coords['theta'], + sc.array(dims=['event'], values=-_expected_exit_angles(sign=-1.0), unit='rad'), + ) + + +@pytest.mark.parametrize( + ('run_norm', 'expected'), + [ + (RunNormalization.none, [0.8, 0.4]), + (RunNormalization.monitor_histogram, [0.2, 0.2]), + (RunNormalization.monitor_integrated, [0.8 / 3, 0.4 / 3]), + ], +) +def test_direct_beam_reduction_integrates_separate_rois(run_norm, expected): + wf = _workflow(run_norm) + wf[SampleIlluminatedFraction] = sc.scalar(0.5) + for run, values in [(SampleRun, [4.0, 8.0]), (ReferenceRun, [2.0, 2.0])]: + wf[WavelengthMonitor[run]] = sc.DataArray( + sc.array(dims=['wavelength'], values=values, unit='counts'), + coords={'wavelength': wf.compute(WavelengthBins)}, + ) + result = wf.compute(ReflectivityOverQ) + assert result.dims == ('Q',) + assert result.unit == sc.units.dimensionless + np.testing.assert_allclose(result.values[:2], expected) + assert_identical( + result.masks['direct_beam'], + sc.array(dims=['Q'], values=[False, False, True]), + ) + assert np.isnan(result.values[2]) + # The independent sample and direct-beam counting uncertainties both contribute. + np.testing.assert_allclose( + result.variances[:2], + np.array(expected) ** 2 * (1 / np.array([40, 20]) + 1 / 100), + ) + + +def test_reduction_is_ratio_of_integrals_when_q_bins_are_merged(): + wf = _workflow() + wf[SampleIlluminatedFraction] = sc.scalar(1.0) + wf[QBins] = sc.array(dims=['Q'], values=[0.04, 0.13], unit='1/angstrom') + result = wf.compute(ReflectivityOverQ) + np.testing.assert_allclose(result.values, [60.0 / 200.0]) + + +def test_zero_monitor_intensity_masks_corresponding_reflectivity_bin(): + wf = _workflow(RunNormalization.monitor_histogram) + wf[SampleIlluminatedFraction] = sc.scalar(1.0) + for run, values in [(SampleRun, [2.0, 2.0]), (ReferenceRun, [0.0, 2.0])]: + wf[WavelengthMonitor[run]] = sc.DataArray( + sc.array(dims=['wavelength'], values=values, unit='counts'), + coords={'wavelength': wf.compute(WavelengthBins)}, + ) + result = wf.compute(ReflectivityOverQ) + assert_identical( + result.masks['direct_beam'], + sc.array(dims=['Q'], values=[False, True, True]), + ) + + +def test_footprint_uses_amor_model_and_requires_only_sample_sizes(): + wf = _workflow() + wf[SampleSize[SampleRun]] = sc.scalar(10.0, unit='mm') + wf[BeamSize[SampleRun]] = sc.scalar(10.0 * np.sin(np.deg2rad(1.0)), unit='mm') + result = wf.compute(ReflectivityOverQ) + # The two wavelengths in the same pixel have different corrected footprints. + projected_size_over_beam = np.sin(_expected_exit_angles()[:2]) / np.sin( + np.deg2rad(1.0) + ) + fraction = sc.erf( + sc.array( + dims=['event'], values=projected_size_over_beam / np.sqrt(8.0 * np.log(2.0)) + ) + ).values + np.testing.assert_allclose(result.values[:2], np.array([0.4, 0.2]) / fraction[::-1]) From 1dfbd6917cce420db12848d897ccd3b37509963b Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Tue, 15 Sep 2026 17:49:07 +0200 Subject: [PATCH 04/16] fix: clean up --- .../src/ess/freia/conversions.py | 41 ++-- .../src/ess/freia/corrections.py | 5 +- .../src/ess/freia/maskings.py | 6 +- .../essreflectometry/src/ess/freia/mcstas.py | 35 ++- .../src/ess/freia/normalization.py | 35 +-- .../essreflectometry/src/ess/freia/types.py | 4 +- .../src/ess/freia/workflow.py | 1 + .../tests/freia/conversions_test.py | 48 ++++ .../tests/freia/mcstas_test.py | 194 ++++------------ .../tests/freia/workflow_test.py | 208 ++++-------------- 10 files changed, 184 insertions(+), 393 deletions(-) create mode 100644 packages/essreflectometry/tests/freia/conversions_test.py diff --git a/packages/essreflectometry/src/ess/freia/conversions.py b/packages/essreflectometry/src/ess/freia/conversions.py index 036b277a7..4a6a04239 100644 --- a/packages/essreflectometry/src/ess/freia/conversions.py +++ b/packages/essreflectometry/src/ess/freia/conversions.py @@ -2,7 +2,7 @@ # Copyright (c) 2026 Scipp contributors (https://github.com/scipp) import scipp as sc from ess.reduce.nexus.types import GravityVector, Position -from scippneutron.conversion import beamline, graph +from scippneutron.conversion import graph, tof from scippnexus import NXsample, NXsource from ..reflectometry.conversions import reflectometry_q @@ -11,40 +11,25 @@ def theta( - incident_beam: sc.Variable, scattered_beam: sc.Variable, wavelength: sc.Variable, gravity: sc.Variable, sample_surface_normal: sc.Variable, ) -> sc.Variable: - """Signed, gravity-corrected exit angle above the sample plane. + """Signed, gravity-corrected angle above the sample plane. - ScippNeutron reconstructs the outgoing direction at the sample. Project - that direction onto the sample normal to retain the sign and support a - tilted sample. Its reflectometry-specific scattering_angle_in_yz_plane - returns an unsigned angle, which cannot distinguish the direct beam. - - The horizontal beam direction only defines a coordinate basis here; it - does not specify an incident angle. For Q and footprint we still assume - specular reflection, so the incidence angle equals this exit angle. + Approximate the flight time using the straight sample-to-detector distance, + as in ScippNeutron's gravity correction. Positive angles point toward the + sample surface normal. """ - basis = beamline.beam_aligned_unit_vectors(incident_beam, gravity) - x, y, z = (basis[f'beam_aligned_unit_{axis}'] for axis in 'xyz') - # Use the horizontal reference axis: the source-to-sample line in FREIA - # is tilted and does not describe the incident direction at the sample. - angles = beamline.scattering_angles_with_gravity( - incident_beam=z * sc.scalar(1.0, unit='m'), - scattered_beam=scattered_beam, - wavelength=wavelength, - gravity=gravity, + flight_time = tof.tof_from_wavelength( + wavelength=wavelength, Ltotal=sc.norm(scattered_beam) + ).to(unit='s') + outgoing_beam = scattered_beam - (0.5 * gravity * flight_time**2).to( + unit=scattered_beam.unit ) - polar, azimuth = angles['two_theta'], angles['phi'] normal = sample_surface_normal / sc.norm(sample_surface_normal) - return sc.asin( - sc.sin(polar) - * (sc.dot(normal, x) * sc.cos(azimuth) + sc.dot(normal, y) * sc.sin(azimuth)) - + sc.dot(normal, z) * sc.cos(polar) - ) + return sc.asin(sc.dot(outgoing_beam, normal) / sc.norm(outgoing_beam)) def coordinate_transformation_graph( @@ -53,7 +38,7 @@ def coordinate_transformation_graph( sample_surface_normal: SampleSurfaceNormal[RunType], gravity: GravityVector, ) -> CoordTransformationGraph[RunType]: - """Build a specular conversion graph independent of detector pixel layout.""" + """Build a specular conversion graph.""" length = sc.norm(sample_surface_normal) if ( not sc.isfinite(length).value @@ -75,7 +60,7 @@ def add_coords( da: WavelengthDetector[RunType], graph: CoordTransformationGraph[RunType], ) -> QDetector[RunType]: - """Add Q without requiring an ROI, monitor, reference, or footprint inputs.""" + """Add specular Q and the gravity-corrected angle to detector events.""" return QDetector[RunType]( da.transform_coords( ( diff --git a/packages/essreflectometry/src/ess/freia/corrections.py b/packages/essreflectometry/src/ess/freia/corrections.py index d333589f1..0cf9885a1 100644 --- a/packages/essreflectometry/src/ess/freia/corrections.py +++ b/packages/essreflectometry/src/ess/freia/corrections.py @@ -117,10 +117,7 @@ def prepare_sample( illuminated_fraction = sc.bins_like(sample, illuminated_fraction) valid = sc.isfinite(illuminated_fraction) & (illuminated_fraction > sc.scalar(0.0)) valid &= illuminated_fraction <= sc.scalar(1.0) - sample = sample.bins.assign_masks( - footprint=~valid, - non_reflected=sample.bins.coords['theta'] <= sc.scalar(0.0, unit='rad'), - ) + sample = sample.bins.assign_masks(footprint=~valid) fraction = sc.where(valid, illuminated_fraction, sc.scalar(1.0)) return Sample(sample / fraction) diff --git a/packages/essreflectometry/src/ess/freia/maskings.py b/packages/essreflectometry/src/ess/freia/maskings.py index 8dea32f8d..84eec42a3 100644 --- a/packages/essreflectometry/src/ess/freia/maskings.py +++ b/packages/essreflectometry/src/ess/freia/maskings.py @@ -6,12 +6,12 @@ from .types import DetectorRegionOfInterest, QDetector -def select_events( +def add_masks( da: QDetector[RunType], roi: DetectorRegionOfInterest[RunType], wavelength_bins: WavelengthBins, ) -> RunUnnormalizedData[RunType]: - """Select a peak using pixel or event coordinates and wavelength.""" + """Mask events outside the ROI and wavelength range.""" masks = {} event_masks = {} for name, (low, high) in roi.items(): @@ -36,4 +36,4 @@ def select_events( ) -providers = (select_events,) +providers = (add_masks,) diff --git a/packages/essreflectometry/src/ess/freia/mcstas.py b/packages/essreflectometry/src/ess/freia/mcstas.py index c35421d06..17624dd6e 100644 --- a/packages/essreflectometry/src/ess/freia/mcstas.py +++ b/packages/essreflectometry/src/ess/freia/mcstas.py @@ -1,10 +1,6 @@ # SPDX-License-Identifier: BSD-3-Clause # Copyright (c) 2026 Scipp contributors (https://github.com/scipp) -"""Adapters for FREIA McStas files. - -McStas conventions and the fixed WFM simulation configuration belong here. -The rest of the workflow uses the standard ESSreduce domain types. -""" +"""Adapters for FREIA McStas files.""" import re from pathlib import Path @@ -171,6 +167,10 @@ class _ChopperParameters(TypedDict): }, } +# Non-WFM settings in run 265080: stopWBC=0, stopWFM=1. The bandwidth +# choppers use the same settings as in WFM mode; the PSC disks are stopped. +_NON_WFM_PARAMETERS = {name: _WFM_PARAMETERS[name] for name in ('WBC1', 'WBC2', 'WBC3')} + def wfm_choppers() -> DiskChoppers[RunType]: """Construct the fixed WFM cascade used by the initial FREIA simulations. @@ -183,6 +183,18 @@ def wfm_choppers() -> DiskChoppers[RunType]: The analytical cascade projects them onto the z axis, approximating the curved guide and finite beam width by a central ray. """ + return _make_choppers(_WFM_PARAMETERS) + + +def non_wfm_choppers() -> DiskChoppers[RunType]: + """Construct the non-WFM cascade with only the three bandwidth choppers. + + Insert this provider into the workflow to replace the default WFM cascade. + """ + return _make_choppers(_NON_WFM_PARAMETERS) + + +def _make_choppers(settings: dict[str, _ChopperParameters]) -> DiskChoppers[RunType]: return DiskChoppers[RunType]( { name: DiskChopper( @@ -197,7 +209,7 @@ def wfm_choppers() -> DiskChoppers[RunType]: dims=['cutout'], values=parameters['close'], unit='deg' ), ) - for name, parameters in _WFM_PARAMETERS.items() + for name, parameters in settings.items() } ) @@ -217,11 +229,11 @@ def load_mcstas( *, pulse_period: sc.Variable | None = None, ) -> sc.DataArray: - """Load weighted events and pixel geometry from the final FREIA detector. + """Load weighted detector events and pixel geometry. Only the selected component is read. The expected output is the ``mantid banana`` event list, including ``p``, ``t``, ``id`` and pixel - geometry. Histogram-only and upstream debug monitors are not substitutes. + geometry. Empty detector pixels are retained, and weighted-event variances are ``p**2``. Arrival times are split into ``event_time_zero`` and ``event_time_offset`` using ``pulse_period``, which defaults to the ESS period of 1/14 s. @@ -286,7 +298,7 @@ def load_mcstas_provider( geometry: EmptyDetector[RunType], pulse_period: PulsePeriod, ) -> RawDetector[RunType]: - """Provide final-detector events, reusing the workflow's pixel geometry.""" + """Load detector events using the workflow's pixel geometry.""" with _open_mcstas(filename) as data: return RawDetector[RunType]( _load_events(data, detector_name, pulse_period, geometry) @@ -325,8 +337,6 @@ def mcstas_sample_surface_normal( def load_mcstas_monitor(filename: str | Path, monitor_name: str) -> sc.DataArray: """Read a selected one-dimensional McStas L_monitor histogram. - Select an incident monitor upstream of the sample. In particular, - SampleLambda in FREIA_surface_test.instr is downstream of the sample. McStas stores bin-integrated intensities and standard errors; xlimits supplies the bounds of the uniformly spaced wavelength bins. """ @@ -383,7 +393,8 @@ def _detector_geometry(data, detector_name) -> sc.DataArray: data.file_object.get_pixels_entry(detector_name), dtype='int64' ).ravel() else: - # Histogram axes also describe geometry; no intensities are read. + # McStasToX's pixel reader requires BINS. Older histogram-only files + # need their axes read separately. geometry = data.file_object.get_geometry_dict(detector_name) if geometry['shape'] != 'banana': raise ValueError(f'Expected banana geometry for {detector_name!r}.') diff --git a/packages/essreflectometry/src/ess/freia/normalization.py b/packages/essreflectometry/src/ess/freia/normalization.py index 39f8ce41f..9aeb35326 100644 --- a/packages/essreflectometry/src/ess/freia/normalization.py +++ b/packages/essreflectometry/src/ess/freia/normalization.py @@ -1,7 +1,6 @@ # SPDX-License-Identifier: BSD-3-Clause # Copyright (c) 2026 Scipp contributors (https://github.com/scipp) import scipp as sc -from ess.reduce.nexus.types import GravityVector from ..reflectometry.conversions import reflectometry_q from ..reflectometry.types import ( @@ -11,37 +10,17 @@ ReferenceRun, ReflectivityOverQ, Sample, - SampleRun, ) -from .conversions import theta -from .types import SampleSurfaceNormal def evaluate_direct_beam( direct_beam: ReducibleData[ReferenceRun], - sample_surface_normal: SampleSurfaceNormal[SampleRun], - gravity: GravityVector, ) -> Reference: - """Map the direct-beam ROI to Q after reflection at the sample surface. - - A direct ray follows the incident direction. Specular reflection reverses - its normal component, so its exit angle is the negative of its angle above - the sample plane. Using the sample run's normal also supports a changed - sample orientation. No footprint or supermirror correction is applied. - """ - angle = -theta( - incident_beam=direct_beam.coords['incident_beam'], - scattered_beam=direct_beam.coords['position'] - - direct_beam.coords['sample_position'], - wavelength=direct_beam.bins.coords['wavelength'], - gravity=gravity, - sample_surface_normal=sample_surface_normal, - ) - reference = direct_beam.bins.assign_coords( - theta=angle, Q=reflectometry_q(direct_beam.bins.coords['wavelength'], angle) - ) + """Compute reference Q using the direct beam's incidence angle.""" + theta = -direct_beam.bins.coords['theta'] + wavelength = direct_beam.bins.coords['wavelength'] return Reference( - reference.bins.assign_masks(non_incident=angle <= sc.scalar(0.0, unit='rad')) + direct_beam.bins.assign_coords(Q=reflectometry_q(wavelength, theta)) ) @@ -60,11 +39,9 @@ def reduce_sample_over_q( valid = sc.isfinite(denominator.data) & ( denominator.data > sc.scalar(0.0, unit=denominator.unit) ) - # Keep unusable bins explicit instead of producing infinities at empty bins. - norm = sc.where( - valid, denominator.data, sc.scalar(float('nan'), unit=denominator.unit) + return ReflectivityOverQ( + (numerator / denominator.data).assign_masks(direct_beam=~valid) ) - return ReflectivityOverQ((numerator / norm).assign_masks(direct_beam=~valid)) providers = (evaluate_direct_beam, reduce_sample_over_q) diff --git a/packages/essreflectometry/src/ess/freia/types.py b/packages/essreflectometry/src/ess/freia/types.py index 3b36b3b53..c5e7a401a 100644 --- a/packages/essreflectometry/src/ess/freia/types.py +++ b/packages/essreflectometry/src/ess/freia/types.py @@ -24,14 +24,14 @@ class SampleSurfaceNormal(sciline.Scope[RunType, sc.Variable], sc.Variable): class QDetector(sciline.Scope[RunType, sc.DataArray], sc.DataArray): - """Detector events with specular Q and signed angle above the sample surface.""" + """Detector events with specular Q and signed angle to the sample surface.""" class DetectorRegionOfInterest(sciline.Scope[RunType, dict], dict): """Pixel or event coordinates mapped to inclusive (lower, upper) bounds. Select corresponding reflected and direct peaks separately for SampleRun and - ReferenceRun, for example using ``theta`` and ``height``. + ReferenceRun, for example using signed ``theta`` and ``height``. An empty dictionary explicitly selects the entire detector. """ diff --git a/packages/essreflectometry/src/ess/freia/workflow.py b/packages/essreflectometry/src/ess/freia/workflow.py index 0af66d670..c2562252b 100644 --- a/packages/essreflectometry/src/ess/freia/workflow.py +++ b/packages/essreflectometry/src/ess/freia/workflow.py @@ -108,6 +108,7 @@ def FreiaWorkflow( requires separate sample/direct-beam ROIs, wavelength and Q bins, and beam and sample sizes for the footprint correction. The reference run must be a measurement without a sample, taken with matching slit and chopper settings. + Set its ``SampleSurfaceNormal`` to the sample run's orientation. Monitor normalization requires an incident monitor selected through ``NeXusName[IncidentMonitor]``, or supplied as ``WavelengthMonitor[RunType]``. diff --git a/packages/essreflectometry/tests/freia/conversions_test.py b/packages/essreflectometry/tests/freia/conversions_test.py new file mode 100644 index 000000000..e61d0909c --- /dev/null +++ b/packages/essreflectometry/tests/freia/conversions_test.py @@ -0,0 +1,48 @@ +# SPDX-License-Identifier: BSD-3-Clause +# Copyright (c) 2026 Scipp contributors (https://github.com/scipp) +import numpy as np +import scipp as sc +from scipp.testing import assert_allclose + +from ess.freia.conversions import theta + + +def test_theta_relative_to_sample_surface(): + result = theta( + scattered_beam=sc.vectors( + dims=['event'], values=[[0, 1, 1], [0, -1, 1]], unit='m' + ), + wavelength=sc.scalar(4.0, unit='angstrom'), + gravity=sc.vector([0.0, 0.0, 0.0], unit='m/s^2'), + sample_surface_normal=sc.vector([0.0, np.sqrt(3) / 2, -0.5]), + ) + + assert_allclose( + result, + sc.array(dims=['event'], values=[15.0, -75.0], unit='deg').to(unit='rad'), + ) + + +def test_theta_with_gravity(): + # Horizontal rays at 1000 and 500 m/s fall 0.49 and 1.96 mm over 10 m. + wavelength = ( + sc.constants.h + / sc.constants.m_n + / sc.array(dims=['event'], values=[1000.0, 500.0], unit='m/s') + ) + result = theta( + scattered_beam=sc.vectors( + dims=['event'], + values=[[0, -0.0004903325, 10], [0, -0.00196133, 10]], + unit='m', + ), + wavelength=wavelength.to(unit='angstrom'), + gravity=sc.vector([0.0, -9.80665, 0.0], unit='m/s^2'), + sample_surface_normal=sc.vector([0.0, 1.0, 0.0]), + ) + + assert_allclose( + result, + sc.zeros(dims=['event'], shape=[2], unit='rad'), + atol=sc.scalar(1e-10, unit='rad'), + ) diff --git a/packages/essreflectometry/tests/freia/mcstas_test.py b/packages/essreflectometry/tests/freia/mcstas_test.py index c8754c063..d83073b74 100644 --- a/packages/essreflectometry/tests/freia/mcstas_test.py +++ b/packages/essreflectometry/tests/freia/mcstas_test.py @@ -1,35 +1,18 @@ # SPDX-License-Identifier: BSD-3-Clause # Copyright (c) 2026 Scipp contributors (https://github.com/scipp) - -import math +from pathlib import Path import h5py import numpy as np import pytest import scipp as sc -from ess.reduce.nexus.types import ( - DiskChoppers, - Filename, - NeXusDetectorName, - NeXusName, - RawDetector, -) -from ess.reduce.unwrap import ( - DetectorLtotal, - FrameUnwrapBackend, - LookupTable, - PulsePeriod, - SourceBounds, - WavelengthDetector, -) +from ess.reduce.nexus.types import Filename, NeXusName, RawDetector from scipp.testing import assert_allclose, assert_identical -from scippneutron.chopper import DiskChopper -from scippnexus import NXdetector from ess.freia import FreiaMcStasWorkflow -from ess.freia.mcstas import load_mcstas +from ess.freia.mcstas import mcstas_detector_geometry from ess.freia.types import IncidentMonitor, QDetector, WavelengthMonitor -from ess.reflectometry.types import ReferenceRun, SampleRun +from ess.reflectometry.types import SampleRun def _component(components, name, position): @@ -40,8 +23,8 @@ def _component(components, name, position): @pytest.fixture -def mcstas_file(tmp_path): - """Small on-disk McStas file read by the real mcstastox library.""" +def mcstas_file(tmp_path: Path) -> Path: + """Minimal file for the McStasToX reader; FREIA has no published event fixture.""" filename = tmp_path / 'freia.h5' with h5py.File(filename, 'w') as f: entry = f.create_group('entry1') @@ -49,17 +32,10 @@ def mcstas_file(tmp_path): simulation = entry.create_group('simulation') simulation.attrs['program'] = np.bytes_('3.7.18, git') simulation.create_group('Param') - instrument = entry.create_group('instrument') - components = instrument.create_group('components') + components = entry.create_group('instrument/components') _component(components, 'Source', [0.0, 0.0, 0.0]) _component(components, 'Arm_Sample', [0.0, 0.0, 20.0]) - # A debug monitor must never be included in the detector event sum. - debug = _component(components, 'Slit_event', [0.0, 0.0, 19.0]) - debug_output = debug.create_group('output/events') - debug_output.attrs['variables'] = np.bytes_('p t id') - debug_output['events'] = [[999.0, 0.001, 0.0]] detector = _component(components, 'Multiblade', [0.0, 0.0, 20.0]) - # Rotate the banana into the vertical scattering plane. detector['Rotation'][...] = [[0, 1, 0], [-1, 0, 0], [0, 0, 1]] geometry = detector.create_group('Geometry') geometry.attrs['Shape identifier'] = np.bytes_('4') @@ -74,7 +50,6 @@ def mcstas_file(tmp_path): bins.attrs[key] = np.bytes_(value) bins['theta'] = [1.0, 2.0] bins['height'] = [-0.001, 0.001] - # IDs need not be contiguous or sorted in geometry order. bins['pixels'] = [[12, 10], [99, 101]] output = detector.create_group('output/detector_events') output.attrs['variables'] = np.bytes_('p t id') @@ -87,11 +62,18 @@ def mcstas_file(tmp_path): return filename -def test_loader_reads_only_final_detector_and_retains_empty_pixels(mcstas_file): - detector = load_mcstas(mcstas_file) - # Check pixel associations without requiring the loader to return a given order. - image = sc.sort(detector.bins.sum(), 'pixel_id') - np.testing.assert_array_equal(image.coords['pixel_id'].values, [10, 12, 99, 101]) +def test_load_detector_and_compute_q(mcstas_file): + workflow = FreiaMcStasWorkflow() + workflow[Filename[SampleRun]] = mcstas_file + + result = workflow.compute((RawDetector[SampleRun], QDetector[SampleRun])) + + raw = result[RawDetector[SampleRun]] + image = sc.sort(raw.bins.sum(), 'pixel_id') + assert_identical( + image.coords['pixel_id'], + sc.array(dims=['pixel_id'], values=[10, 12, 99, 101], unit=None), + ) assert_identical( image.data, sc.array( @@ -101,149 +83,53 @@ def test_loader_reads_only_final_detector_and_retains_empty_pixels(mcstas_file): unit='counts', ), ) - offsets = detector.bins.constituents['data'].coords['event_time_offset'] - assert_allclose( - offsets.to(unit='s'), - sc.full(sizes=offsets.sizes, value=0.025, unit='s'), - ) - # Pixel 12 is the first banana pixel, rotated into the global frame. + events = raw.bins.constituents['data'] assert_allclose( - image.coords['position'][1], - sc.vector( - [0.001, 3 * math.sin(math.pi / 180), 20 + 3 * math.cos(math.pi / 180)], - unit='m', - ), + events.coords['event_time_offset'], + sc.full(sizes=events.sizes, value=0.025, unit='s'), ) + q = result[QDetector[SampleRun]].bins.constituents['data'].coords['Q'] + assert sc.isfinite(q).all().value -def test_loader_rejects_missing_detector_without_using_debug_events(mcstas_file): +def test_load_histogram_geometry(mcstas_file): with h5py.File(mcstas_file, 'r+') as f: - del f['entry1/instrument/components/0003_Multiblade'] - with pytest.raises(ValueError, match='No Mantid detector events'): - load_mcstas(mcstas_file) - - -@pytest.mark.parametrize('pixel_id', [999, 10.5]) -def test_loader_rejects_events_missing_from_pixel_map(mcstas_file, pixel_id): - with h5py.File(mcstas_file, 'r+') as f: - events = f[ - 'entry1/instrument/components/0003_Multiblade/output/detector_events/events' - ] - events[0, 2] = pixel_id - with pytest.raises(ValueError, match='pixel IDs absent from the pixel map'): - load_mcstas(mcstas_file) - - -def test_lookup_table_uses_histogram_geometry_without_loading_events(mcstas_file): - with h5py.File(mcstas_file, 'r+') as f: - components = f['entry1/instrument/components'] - components.move('0003_Multiblade', '0003_Detector') - detector = components['0003_Detector'] - detector['Position'][...] = [0.0, -0.25, 20.0] + detector = f['entry1/instrument/components/0002_Multiblade'] del detector['output'] histogram = detector.create_group('output/histogram') histogram.attrs['xlabel'] = np.bytes_('theta [deg]') histogram.attrs['ylabel'] = np.bytes_('Height [cm]') - histogram['theta__deg_'] = [0.0, 15.0] + histogram['theta__deg_'] = [0.0, 90.0] histogram['Height__cm_'] = [-25, 25] - # No intensities or event arrays are needed anywhere in the file. - del components['0002_Slit_event'] - - workflow = FreiaMcStasWorkflow() - workflow[Filename[SampleRun]] = str(mcstas_file) - workflow[NeXusDetectorName] = 'Detector' - results = workflow.compute( - (DetectorLtotal[SampleRun], LookupTable[SampleRun, NXdetector]) - ) - # The detector's vertical offset shortens the distance at larger angles. - expected = [ - 20 + math.sqrt(9.125 - 1.5 * math.sin(angle)) for angle in (0, math.pi / 12) - ] - assert_allclose( - results[DetectorLtotal[SampleRun]], - sc.array(dims=['pixel_id'], values=expected * 2, unit='m'), - ) - table = results[LookupTable[SampleRun, NXdetector]].array - assert sc.isfinite(table.data).any().value - # The histogram geometry must not make it possible to load fake events. - with pytest.raises(ValueError, match='No Mantid detector events'): - workflow.compute(RawDetector[SampleRun]) + detector = mcstas_detector_geometry(mcstas_file, 'Multiblade') -@pytest.mark.parametrize('run', [SampleRun, ReferenceRun]) -def test_workflow_loads_and_unwraps_detector_with_generic_providers(mcstas_file, run): - workflow = FreiaMcStasWorkflow() - workflow[Filename[run]] = str(mcstas_file) - workflow[FrameUnwrapBackend] = FrameUnwrapBackend.scipy - # Override WFM with a single disk for an independently calculable wavelength. - workflow[DiskChoppers[run]] = { - 'test_chopper': DiskChopper( - frequency=sc.scalar(14.0, unit='Hz'), - beam_position=sc.scalar(0.0, unit='deg'), - phase=sc.scalar(0.0, unit='deg'), - axle_position=sc.vector([0.0, 0.0, 5.0], unit='m'), - slit_begin=sc.array(dims=['cutout'], values=[279.68], unit='deg'), - slit_end=sc.array(dims=['cutout'], values=[359.68], unit='deg'), - ), - } - workflow[SourceBounds] = SourceBounds( - time=(sc.scalar(0.9, unit='ms'), sc.scalar(1.1, unit='ms')), - wavelength=(sc.scalar(0.5, unit='angstrom'), sc.scalar(12.0, unit='angstrom')), - ) - result = workflow.compute( - (RawDetector[run], WavelengthDetector[run], QDetector[run]) - ) - raw = result[RawDetector[run]] - unwrapped = result[WavelengthDetector[run]] - assert_allclose(raw.bins.sum().data, unwrapped.bins.sum().data) - assert 'Q' in result[QDetector[run]].bins.coords - wavelength = unwrapped.bins.constituents['data'].coords['wavelength'] - # Arrival time 25 ms minus emission time 1 ms, flight path about 23 m. - expected = ( - sc.constants.h - / sc.constants.m_n - * sc.scalar(24.0, unit='ms') - / sc.scalar(23.0, unit='m') - ).to(unit='angstrom') assert_allclose( - wavelength, - sc.full(sizes=wavelength.sizes, value=expected.value, unit=expected.unit), - rtol=sc.scalar(0.002), - ) - - -def test_loader_uses_workflow_pulse_period(mcstas_file): - workflow = FreiaMcStasWorkflow() - workflow[Filename[SampleRun]] = str(mcstas_file) - workflow[PulsePeriod] = sc.scalar(50.0, unit='ms') - events = workflow.compute(RawDetector[SampleRun]).bins.constituents['data'] - events = sc.sort(events, 'event_time_zero') - # The event at 25 ms + 1/14 s falls in the second 50 ms pulse period. - assert_identical( - events.coords['event_time_zero'].to(unit='ns'), - sc.datetimes(dims=['event'], values=[0, 0, 50_000_000], unit='ns'), - ) - assert_allclose( - events.coords['event_time_offset'].to(unit='s'), - sc.array( - dims=['event'], values=[0.025, 0.025, 0.025 + 1 / 14 - 0.05], unit='s' + detector.coords['position'], + sc.vectors( + dims=['pixel_id'], + values=[[0.25, 0, 23], [0.25, 3, 20], [-0.25, 0, 23], [-0.25, 3, 20]], + unit='m', ), ) -def test_load_selected_wavelength_monitor_with_bin_edges_and_variances(mcstas_file): +def test_load_monitor(mcstas_file): with h5py.File(mcstas_file, 'r+') as f: - components = f['entry1/instrument/components'] - monitor = _component(components, 'IncidentLambda', [0.0, 0.0, 19.0]) + monitor = _component( + f['entry1/instrument/components'], 'IncidentLambda', [0, 0, 19] + ) histogram = monitor.create_group('output/spectrum') histogram.attrs['xvar'] = np.bytes_('L') histogram.attrs['xlimits'] = np.bytes_('1 5') histogram['data'] = [10.0, 20.0] histogram['errors'] = [3.0, 4.0] workflow = FreiaMcStasWorkflow() - workflow[Filename[SampleRun]] = str(mcstas_file) + workflow[Filename[SampleRun]] = mcstas_file workflow[NeXusName[IncidentMonitor]] = 'IncidentLambda' + result = workflow.compute(WavelengthMonitor[SampleRun]) + assert_identical( result, sc.DataArray( diff --git a/packages/essreflectometry/tests/freia/workflow_test.py b/packages/essreflectometry/tests/freia/workflow_test.py index 702e2ce9a..319ab8ce5 100644 --- a/packages/essreflectometry/tests/freia/workflow_test.py +++ b/packages/essreflectometry/tests/freia/workflow_test.py @@ -1,30 +1,25 @@ # SPDX-License-Identifier: BSD-3-Clause # Copyright (c) 2026 Scipp contributors (https://github.com/scipp) -# Small known intensities exercise the complete reduction graph without a LUT. import numpy as np -import pytest import scipp as sc -from ess.reduce import workflow as reduce_workflow -from ess.reduce.nexus.types import Position +from ess.reduce.nexus.types import GravityVector, Position from scipp.testing import assert_allclose, assert_identical from scippnexus import NXsample, NXsource -from ess import freia +from ess.freia import FreiaWorkflow from ess.freia.corrections import RunNormalization from ess.freia.types import ( DetectorRegionOfInterest, - QDetector, SampleIlluminatedFraction, SampleSurfaceNormal, WavelengthMonitor, ) +from ess.reflectometry.corrections import footprint_on_sample from ess.reflectometry.types import ( BeamSize, QBins, - Reference, ReferenceRun, ReflectivityOverQ, - RunUnnormalizedData, SampleRun, SampleSize, WavelengthBins, @@ -32,26 +27,14 @@ ) -def test_freia_workflow_registers_run_normalization_variants(): - for wf in ( - freia.FreiaMcStasUnnormalizedWorkflow, - freia.FreiaMcStasMonitorHistogramWorkflow, - freia.FreiaMcStasMonitorIntegratedWorkflow, - freia.FreiaMcStasProtonChargeWorkflow, - freia.FreiaUnnormalizedWorkflow, - freia.FreiaMonitorHistogramWorkflow, - freia.FreiaMonitorIntegratedWorkflow, - freia.FreiaProtonChargeWorkflow, - ): - assert wf in reduce_workflow.workflow_registry - - -def _detector(angles, weights): - position = np.array([1.0, 2.0, 3.0]) + 3.0 * np.column_stack( - [np.zeros(len(angles)), np.sin(np.deg2rad(angles)), np.cos(np.deg2rad(angles))] - ) - return ( - sc.DataArray( +def _make_workflow(run_norm): + wf = FreiaWorkflow(run_norm=run_norm) + wf[GravityVector] = sc.vector([0.0, 0.0, 0.0], unit='m/s^2') + for run, sign, weights in [ + (SampleRun, 1, [20.0, 40.0, 999.0]), + (ReferenceRun, -1, [100.0, 100.0, 999.0]), + ]: + events = sc.DataArray( sc.array(dims=['event'], values=weights, variances=weights, unit='counts'), coords={ 'wavelength': sc.array( @@ -59,161 +42,64 @@ def _detector(angles, weights): ), 'pixel_id': sc.array(dims=['event'], values=[0, 0, 1], unit=None), }, + ).group('pixel_id') + # The first pixel is at 30 degrees; the second is outside the ROI. + events.coords['position'] = sc.vectors( + dims=['pixel_id'], + values=[[0, sign * 0.5, np.sqrt(3) / 2], [0, sign * np.sqrt(3) / 2, 0.5]], + unit='m', ) - .group('pixel_id') - .assign_coords( - position=sc.vectors(dims=['pixel_id'], values=position[[0, 2]], unit='m') - ) - ) - - -def _workflow(run_norm=RunNormalization.none): - wf = freia.FreiaWorkflow(run_norm=run_norm) - for run, angles, weights in [ - (SampleRun, [1.0, 1.0, 2.0], [20.0, 40.0, 999.0]), - (ReferenceRun, [-1.0, -1.0, -2.0], [100.0, 100.0, 999.0]), - ]: - wf[WavelengthDetector[run]] = _detector(angles, weights) - wf[Position[NXsample, run]] = sc.vector([1.0, 2.0, 3.0], unit='m') - wf[Position[NXsource, run]] = sc.vector([1.0, 2.0, -17.0], unit='m') + wf[WavelengthDetector[run]] = events + wf[Position[NXsample, run]] = sc.vector([0.0, 0.0, 0.0], unit='m') + wf[Position[NXsource, run]] = sc.vector([0.0, 0.0, -20.0], unit='m') wf[SampleSurfaceNormal[run]] = sc.vector([0.0, 1.0, 0.0]) - wf[DetectorRegionOfInterest[SampleRun]] = { - 'theta': (sc.scalar(0.5, unit='deg'), sc.scalar(1.5, unit='deg')) - } - wf[DetectorRegionOfInterest[ReferenceRun]] = { - 'theta': (sc.scalar(-1.5, unit='deg'), sc.scalar(-0.5, unit='deg')) - } + low, high = sorted([sign * 20.0, sign * 40.0]) + wf[DetectorRegionOfInterest[run]] = { + 'theta': (sc.scalar(low, unit='deg'), sc.scalar(high, unit='deg')) + } wf[WavelengthBins] = sc.array( dims=['wavelength'], values=[1.0, 3.0, 5.0], unit='angstrom' ) - wf[QBins] = sc.array(dims=['Q'], values=[0.04, 0.08, 0.13, 0.3], unit='1/angstrom') + wf[QBins] = sc.array(dims=['Q'], values=[1.0, 2.0, 4.0, 6.0], unit='1/angstrom') return wf -def test_q_uses_sample_plane_and_translated_pixels_without_reduction_inputs(): - wf = _workflow() - result = wf.compute(QDetector[SampleRun]) - q = result.bins.constituents['data'].coords['Q'] - expected = 4 * np.pi * np.sin(_expected_exit_angles()) / [2.0, 4.0, 2.0] - assert_allclose(q, sc.array(dims=['event'], values=expected, unit='1/angstrom')) - assert_allclose( - result.coords['L2'], sc.full(sizes=result.sizes, value=3.0, unit='m') - ) - - -def test_q_follows_tilted_sample_surface(): - wf = _workflow() - tilt = np.deg2rad(0.5) - wf[SampleSurfaceNormal[SampleRun]] = sc.vector([0.0, np.cos(tilt), -np.sin(tilt)]) - result = wf.compute(QDetector[SampleRun]) - assert_allclose( - result.bins.constituents['data'].coords['theta'], - sc.array(dims=['event'], values=_expected_exit_angles() - tilt, unit='rad'), - ) - - -def _expected_exit_angles(sign=1.0): - """Independent calculation from displacement = velocity*time + gravity*time²/2.""" - angle = sign * np.deg2rad([1.0, 1.0, 2.0]) - speed = ( - ( - sc.constants.h - / sc.constants.m_n - / sc.array(dims=['event'], values=[2.0, 4.0, 2.0], unit='angstrom') - ) - .to(unit='m/s') - .values - ) - time = 3.0 / speed - y = 3.0 * np.sin(angle) + 0.5 * sc.constants.g.value * time**2 - return np.arctan2(y, 3.0 * np.cos(angle)) - - -def test_theta_roi_distinguishes_wavelengths_in_the_same_pixel(): - wf = _workflow() - wf[DetectorRegionOfInterest[SampleRun]] = { - 'theta': (sc.scalar(1.0005, unit='deg'), sc.scalar(1.001, unit='deg')), - } - selected = wf.compute(RunUnnormalizedData[SampleRun]) - assert selected.bins.sum().sum().value == 40.0 - - -def test_direct_beam_gravity_is_corrected_before_specular_mapping(): - wf = _workflow() - reference = wf.compute(Reference) - assert_allclose( - reference.bins.constituents['data'].coords['theta'], - sc.array(dims=['event'], values=-_expected_exit_angles(sign=-1.0), unit='rad'), - ) - - -@pytest.mark.parametrize( - ('run_norm', 'expected'), - [ - (RunNormalization.none, [0.8, 0.4]), - (RunNormalization.monitor_histogram, [0.2, 0.2]), - (RunNormalization.monitor_integrated, [0.8 / 3, 0.4 / 3]), - ], -) -def test_direct_beam_reduction_integrates_separate_rois(run_norm, expected): - wf = _workflow(run_norm) - wf[SampleIlluminatedFraction] = sc.scalar(0.5) +def test_reduce_reflectivity_with_monitor_and_footprint(): + wf = _make_workflow(RunNormalization.monitor_histogram) + sample_size = sc.scalar(10.0, unit='mm') + beam_size = sc.scalar(5.0, unit='mm') + wf[SampleSize[SampleRun]] = sample_size + wf[BeamSize[SampleRun]] = beam_size for run, values in [(SampleRun, [4.0, 8.0]), (ReferenceRun, [2.0, 2.0])]: wf[WavelengthMonitor[run]] = sc.DataArray( sc.array(dims=['wavelength'], values=values, unit='counts'), coords={'wavelength': wf.compute(WavelengthBins)}, ) + result = wf.compute(ReflectivityOverQ) - assert result.dims == ('Q',) - assert result.unit == sc.units.dimensionless - np.testing.assert_allclose(result.values[:2], expected) + + fraction = footprint_on_sample(sc.scalar(30.0, unit='deg'), beam_size, sample_size) + # Q bins contain wavelengths 4 and 2 angstrom, respectively. + expected = ( + sc.array(dims=['Q'], values=[0.1, 0.1], variances=[0.00035, 0.0006]) / fraction + ) + assert_allclose(result['Q', :2].data, expected) assert_identical( result.masks['direct_beam'], sc.array(dims=['Q'], values=[False, False, True]), ) - assert np.isnan(result.values[2]) - # The independent sample and direct-beam counting uncertainties both contribute. - np.testing.assert_allclose( - result.variances[:2], - np.array(expected) ** 2 * (1 / np.array([40, 20]) + 1 / 100), - ) -def test_reduction_is_ratio_of_integrals_when_q_bins_are_merged(): - wf = _workflow() +def test_rebinning_integrates_before_dividing(): + wf = _make_workflow(RunNormalization.none) wf[SampleIlluminatedFraction] = sc.scalar(1.0) - wf[QBins] = sc.array(dims=['Q'], values=[0.04, 0.13], unit='1/angstrom') - result = wf.compute(ReflectivityOverQ) - np.testing.assert_allclose(result.values, [60.0 / 200.0]) + wf[QBins] = sc.array(dims=['Q'], values=[1.0, 4.0], unit='1/angstrom') - -def test_zero_monitor_intensity_masks_corresponding_reflectivity_bin(): - wf = _workflow(RunNormalization.monitor_histogram) - wf[SampleIlluminatedFraction] = sc.scalar(1.0) - for run, values in [(SampleRun, [2.0, 2.0]), (ReferenceRun, [0.0, 2.0])]: - wf[WavelengthMonitor[run]] = sc.DataArray( - sc.array(dims=['wavelength'], values=values, unit='counts'), - coords={'wavelength': wf.compute(WavelengthBins)}, - ) result = wf.compute(ReflectivityOverQ) - assert_identical( - result.masks['direct_beam'], - sc.array(dims=['Q'], values=[False, True, True]), - ) - -def test_footprint_uses_amor_model_and_requires_only_sample_sizes(): - wf = _workflow() - wf[SampleSize[SampleRun]] = sc.scalar(10.0, unit='mm') - wf[BeamSize[SampleRun]] = sc.scalar(10.0 * np.sin(np.deg2rad(1.0)), unit='mm') - result = wf.compute(ReflectivityOverQ) - # The two wavelengths in the same pixel have different corrected footprints. - projected_size_over_beam = np.sin(_expected_exit_angles()[:2]) / np.sin( - np.deg2rad(1.0) - ) - fraction = sc.erf( + assert_allclose( + result.data, sc.array( - dims=['event'], values=projected_size_over_beam / np.sqrt(8.0 * np.log(2.0)) - ) - ).values - np.testing.assert_allclose(result.values[:2], np.array([0.4, 0.2]) / fraction[::-1]) + dims=['Q'], values=[60.0 / 200.0], variances=[0.3**2 * (1 / 60 + 1 / 200)] + ), + ) From 1e38d751e6603f3eebbc0f456f4851169e98faf3 Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Tue, 15 Sep 2026 17:58:30 +0200 Subject: [PATCH 05/16] docs: cleanup --- .../freia/freia-mcstas-visualization.ipynb | 61 ++++++++----------- .../user-guide/freia/freia-reflectivity.ipynb | 26 ++++---- .../freia/freia-wavelength-lookup-table.ipynb | 56 ++++++++++------- .../docs/user-guide/freia/index.md | 6 +- 4 files changed, 80 insertions(+), 69 deletions(-) diff --git a/packages/essreflectometry/docs/user-guide/freia/freia-mcstas-visualization.ipynb b/packages/essreflectometry/docs/user-guide/freia/freia-mcstas-visualization.ipynb index c048916af..9932bbc57 100644 --- a/packages/essreflectometry/docs/user-guide/freia/freia-mcstas-visualization.ipynb +++ b/packages/essreflectometry/docs/user-guide/freia/freia-mcstas-visualization.ipynb @@ -7,7 +7,7 @@ "source": [ "# FREIA detector data\n", "\n", - "Visualize detector images, arrival-time distributions, and wavelength spectra. Wavelengths are reconstructed using the WFM chopper settings.\n" + "Visualize detector images, arrival-time distributions, and wavelength spectra. Wavelengths are reconstructed using the WFM chopper settings." ] }, { @@ -22,7 +22,6 @@ "import scipp as sc\n", "\n", "from ess import freia\n", - "from ess.freia import data\n", "from ess.reduce.nexus.types import DiskChoppers, Filename, RawDetector\n", "from ess.reduce.unwrap import (\n", " ChopperFrameSequence,\n", @@ -32,7 +31,7 @@ " WavelengthDetector,\n", ")\n", "from ess.reduce.unwrap.types import KeepEventTimeOffset\n", - "from ess.reflectometry.types import SampleRun\n" + "from ess.reflectometry.types import SampleRun" ] }, { @@ -42,7 +41,7 @@ "source": [ "## Select a run\n", "\n", - "Select a run and set the resolution used to reconstruct wavelengths.\n" + "Set the path to a local event file, relative to this notebook, and the resolution used to reconstruct wavelengths." ] }, { @@ -53,10 +52,10 @@ "outputs": [], "source": [ "freia_mcstas = freia.FreiaMcStasWorkflow(wavelength_from='analytical')\n", - "freia_mcstas[Filename[SampleRun]] = data.freia_mcstas_sample_run()\n", + "freia_mcstas[Filename[SampleRun]] = '../../../265224.h5'\n", "freia_mcstas[KeepEventTimeOffset] = True\n", "freia_mcstas[TimeResolution] = sc.scalar(20.0, unit='us')\n", - "freia_mcstas[DistanceResolution] = sc.scalar(0.1, unit='m')\n" + "freia_mcstas[DistanceResolution] = sc.scalar(0.1, unit='m')" ] }, { @@ -66,9 +65,7 @@ "source": [ "## Inspect the WFM chopper cascade\n", "\n", - "Inspect the chopper timing used to reconstruct wavelengths. To use another configuration, assign its chopper dictionary to `freia_mcstas[DiskChoppers[SampleRun]]`.\n", - "\n", - "The analytical model approximates neutron flight paths.\n" + "Inspect the chopper timing used to reconstruct wavelengths. For a non-WFM run, call `freia_mcstas.insert(freia.mcstas.non_wfm_choppers)` before computing wavelengths. Other configurations can be assigned to `freia_mcstas[DiskChoppers[SampleRun]]`." ] }, { @@ -79,7 +76,7 @@ "outputs": [], "source": [ "choppers = freia_mcstas.compute(DiskChoppers[SampleRun])\n", - "sc.DataGroup(choppers)\n" + "sc.DataGroup(choppers)" ] }, { @@ -90,7 +87,7 @@ "outputs": [], "source": [ "frames = freia_mcstas.compute(ChopperFrameSequence[SampleRun])\n", - "frames.draw()\n" + "frames.draw()" ] }, { @@ -100,7 +97,7 @@ "source": [ "## Load and unwrap the detector events\n", "\n", - "Compute wavelengths from the event arrival times and chopper transmission bands.\n" + "Compute wavelengths from the event arrival times and chopper transmission bands." ] }, { @@ -113,7 +110,7 @@ "results = freia_mcstas.compute((RawDetector[SampleRun], WavelengthDetector[SampleRun]))\n", "raw = results[RawDetector[SampleRun]]\n", "unwrapped = results[WavelengthDetector[SampleRun]]\n", - "raw\n" + "raw" ] }, { @@ -123,7 +120,7 @@ "source": [ "## Detector image and arrival times\n", "\n", - "The image uses `longitude` and `height` in the detector's local frame. Intensities are sums of event weights.\n" + "The image uses `longitude` and `height` in the detector's local frame. Intensities are sums of event weights." ] }, { @@ -134,7 +131,7 @@ "outputs": [], "source": [ "detector_image = raw.hist(longitude=120, height=64, dim=raw.dims)\n", - "pp.plot(detector_image, norm='log', title='FREIA final detector', vmin=1e-1)\n" + "pp.plot(detector_image, norm='log', title='FREIA detector', vmin=1e-1)" ] }, { @@ -144,12 +141,17 @@ "metadata": {}, "outputs": [], "source": [ + "pulse_period = freia_mcstas.compute(PulsePeriod).to(unit='s').value\n", "arrival_times = raw.hist(\n", - " event_time_offset=sc.linspace('event_time_offset', 0.0, freia_mcstas.compute(PulsePeriod).to(unit='s').value, 501, unit='s'),\n", + " event_time_offset=sc.linspace(\n", + " 'event_time_offset', 0.0, pulse_period, 501, unit='s'\n", + " ),\n", " dim=raw.dims,\n", ")\n", - "arrival_times.coords['event_time_offset'] = arrival_times.coords['event_time_offset'].to(unit='ms')\n", - "pp.plot(arrival_times, title='Arrival time within the source period')\n" + "arrival_times.coords['event_time_offset'] = arrival_times.coords[\n", + " 'event_time_offset'\n", + "].to(unit='ms')\n", + "pp.plot(arrival_times, title='Arrival time within the source period')" ] }, { @@ -159,7 +161,7 @@ "source": [ "## Wavelengths from analytical frame unwrapping\n", "\n", - "The chopper cascade and arrival times determine the wavelengths. Events outside the modeled transmission bands or above the workflow's relative wavelength uncertainty threshold have NaN wavelengths and do not contribute to the wavelength histograms. Inspect the assigned-event count when assessing the result.\n" + "The chopper cascade and arrival times determine the wavelengths. Events outside the modeled transmission bands or above the workflow's relative wavelength uncertainty threshold have NaN wavelengths and do not contribute to the wavelength histograms. Inspect the assigned-event count when assessing the result." ] }, { @@ -175,7 +177,7 @@ "\n", "wavelength_bins = sc.linspace('wavelength', 1.0, 12.0, 441, unit='angstrom')\n", "spectrum = unwrapped.hist(wavelength=wavelength_bins, dim=unwrapped.dims)\n", - "pp.plot(spectrum, title='FREIA wavelength spectrum')\n" + "pp.plot(spectrum, title='FREIA wavelength spectrum')" ] }, { @@ -188,18 +190,7 @@ "wavelength_image = unwrapped.hist(\n", " wavelength=wavelength_bins, longitude=120, dim=unwrapped.dims\n", ")\n", - "pp.plot(wavelength_image, norm='log', title='Wavelength across the detector', vmin=1e0)\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "15", - "metadata": {}, - "outputs": [], - "source": [ - "unwrapped.bins.coords['wavelength_error'] = (unwrapped.bins.coords['wavelength'] - unwrapped.bins.coords['wavelength_from_mcstas']) / unwrapped.bins.coords['wavelength_from_mcstas']\n", - "unwrapped.hist(wavelength_error=100, dim=unwrapped.dims).plot()" + "pp.plot(wavelength_image, norm='log', title='Wavelength across the detector', vmin=1e0)" ] } ], @@ -218,8 +209,10 @@ "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.12.13" + "pygments_lexer": "ipython3" + }, + "nbsphinx": { + "execute": "never" } }, "nbformat": 4, diff --git a/packages/essreflectometry/docs/user-guide/freia/freia-reflectivity.ipynb b/packages/essreflectometry/docs/user-guide/freia/freia-reflectivity.ipynb index 9decd8e4b..e4cca52fc 100644 --- a/packages/essreflectometry/docs/user-guide/freia/freia-reflectivity.ipynb +++ b/packages/essreflectometry/docs/user-guide/freia/freia-reflectivity.ipynb @@ -7,9 +7,9 @@ "source": [ "# FREIA reflectivity reduction\n", "\n", - "Reduce a reflected beam using a direct-beam measurement without a sample. The workflow converts events to specular Q, normalizes by an incident wavelength monitor, and divides the selected peaks to obtain R(Q).\n", + "Compute reflectivity from a sample measurement and a direct-beam measurement taken without a sample. The workflow converts events to specular Q, normalizes by an incident wavelength monitor, and divides the selected peaks to obtain R(Q).\n", "\n", - "This example uses local McStas files. The sample and direct-beam runs must have matching slit and chopper settings." + "This example uses local McStas files with WFM chopper settings. The sample and direct-beam runs must have matching slit and chopper settings." ] }, { @@ -31,7 +31,7 @@ " SampleIlluminatedFraction,\n", " WavelengthMonitor,\n", ")\n", - "from ess.reduce.nexus.types import NeXusName\n", + "from ess.reduce.nexus.types import GravityVector, NeXusName\n", "from ess.reflectometry.types import (\n", " BeamSize,\n", " Filename,\n", @@ -52,7 +52,9 @@ "source": [ "## Select the runs\n", "\n", - "Set the sample and no-sample filenames. These paths are relative to the notebook directory. Choose an incident wavelength monitor upstream of the sample, so its intensity is independent of reflectivity. If no monitor is available, use `RunNormalization.none`; the two runs must then already share an exposure and flux scale." + "Set the sample and no-sample filenames. These paths are relative to the notebook directory. Choose an incident wavelength monitor upstream of the sample, so its intensity is independent of reflectivity. If no monitor is available, use `RunNormalization.none`; the two runs must then already share an exposure and flux scale.\n", + "\n", + "These runs were simulated without gravity, so `GravityVector` is set to zero. Use the workflow default for data that includes gravity." ] }, { @@ -63,9 +65,10 @@ "outputs": [], "source": [ "workflow = FreiaMcStasWorkflow(run_norm=RunNormalization.monitor_histogram)\n", - "workflow[Filename[SampleRun]] = '../../../265149.h5'\n", - "workflow[Filename[ReferenceRun]] = '../../../265148.h5'\n", + "workflow[Filename[SampleRun]] = '../../../265224.h5'\n", + "workflow[Filename[ReferenceRun]] = '../../../265222.h5'\n", "workflow[NeXusName[IncidentMonitor]] = 'GuideexitLambda'\n", + "workflow[GravityVector] = sc.vector([0.0, 0.0, 0.0], unit='m/s^2')\n", "\n", "workflow[WavelengthBins] = sc.linspace('wavelength', 2.0, 10.0, 81, unit='angstrom')\n", "workflow[QBins] = sc.geomspace('Q', 0.07, 0.4, 21, unit='1/angstrom')" @@ -78,7 +81,7 @@ "source": [ "## Inspect the reflected and direct peaks\n", "\n", - "The angle θ is measured above the sample surface and includes gravity correction. Reflected beams have positive angles; the direct beams lie below the sample plane. The specular assumption gives $Q = 4\\pi\\sin(\\theta)/\\lambda$.\n", + "The angle θ is measured above the sample surface and is corrected for gravity when enabled. Its sign distinguishes the two sides of the sample plane. The specular assumption gives $Q = 4\\pi\\sin(\\theta)/\\lambda$.\n", "\n", "Inspect both distributions before selecting corresponding peaks." ] @@ -108,7 +111,7 @@ "source": [ "## Select matching regions of interest\n", "\n", - "Select the reflected peak and its corresponding direct peak separately. The example uses the beam near 3.5°. Additional bounds can be set on pixel coordinates such as `height` or `pixel_id`. Wavelength and Q bins are deliberately coarse because these runs have few events." + "Select the reflected peak and its corresponding direct peak separately. The example uses the beam near 3.5°. Additional bounds can be set on pixel coordinates such as `height` or `pixel_id`. Choose wavelength and Q bins to match the available statistics." ] }, { @@ -209,7 +212,7 @@ "source": [ "## Compute reflectivity\n", "\n", - "The workflow maps the direct beam to its corresponding specular Q, integrates the two selected peaks into matching Q bins, and divides their intensities. Both counting uncertainties propagate. Bins with no usable direct-beam intensity are masked." + "The reference uses the negative of the direct beam’s θ to compute its corresponding reflected Q. The workflow integrates the two selected peaks into matching Q bins and divides their intensities. Both counting uncertainties propagate. Bins with no usable direct-beam intensity are masked." ] }, { @@ -235,7 +238,7 @@ "id": "15", "metadata": {}, "source": [ - "The sparse statistics and omitted footprint correction limit the interpretation of this curve. Q resolution, background subtraction, finite-sample corrections, and complete ORSO export are not yet included." + "Footprint correction is disabled in this example. Q resolution, background subtraction, finite-sample corrections, and complete ORSO export are not yet included." ] } ], @@ -254,8 +257,7 @@ "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.12.13" + "pygments_lexer": "ipython3" }, "nbsphinx": { "execute": "never" diff --git a/packages/essreflectometry/docs/user-guide/freia/freia-wavelength-lookup-table.ipynb b/packages/essreflectometry/docs/user-guide/freia/freia-wavelength-lookup-table.ipynb index 133005254..8318f01a1 100644 --- a/packages/essreflectometry/docs/user-guide/freia/freia-wavelength-lookup-table.ipynb +++ b/packages/essreflectometry/docs/user-guide/freia/freia-wavelength-lookup-table.ipynb @@ -7,7 +7,7 @@ "source": [ "# FREIA analytical wavelength lookup table\n", "\n", - "Explore which wavelengths can reach the detector at each arrival time for the WFM chopper configuration. This calculation needs geometry and chopper settings; detector events are not required.\n" + "Explore which wavelengths can reach the detector at each arrival time for the WFM chopper configuration. The table is computed from the instrument geometry and chopper settings." ] }, { @@ -18,12 +18,13 @@ "outputs": [], "source": [ "%matplotlib widget\n", + "from pathlib import Path\n", + "\n", "import plopp as pp\n", "import scipp as sc\n", "from scippnexus import NXdetector\n", "\n", "from ess import freia\n", - "from ess.freia import data\n", "from ess.reduce.nexus.types import DiskChoppers, Filename\n", "from ess.reduce.unwrap import (\n", " ChopperFrameSequence,\n", @@ -32,7 +33,7 @@ " LookupTable,\n", " TimeResolution,\n", ")\n", - "from ess.reflectometry.types import SampleRun\n" + "from ess.reflectometry.types import SampleRun" ] }, { @@ -42,7 +43,9 @@ "source": [ "## Use the WFM configuration\n", "\n", - "Set the input file and lookup-table resolution. This example uses 20 µs time resolution and 10 cm flight-path resolution.\n" + "Set the path to a local file, relative to this notebook, and the lookup-table resolution. This example uses 20 µs time resolution and 10 cm flight-path resolution.\n", + "\n", + "For a non-WFM run, insert `freia.mcstas.non_wfm_choppers` into the workflow after creating it: `freia_mcstas.insert(freia.mcstas.non_wfm_choppers)`." ] }, { @@ -52,12 +55,12 @@ "metadata": {}, "outputs": [], "source": [ - "filename = data.freia_mcstas_sample_run()\n", + "filename = Path('../../../265224.h5')\n", "\n", "freia_mcstas = freia.FreiaMcStasWorkflow(wavelength_from='analytical')\n", "freia_mcstas[Filename[SampleRun]] = filename\n", "freia_mcstas[TimeResolution] = sc.scalar(20.0, unit='us')\n", - "freia_mcstas[DistanceResolution] = sc.scalar(0.1, unit='m')\n" + "freia_mcstas[DistanceResolution] = sc.scalar(0.1, unit='m')" ] }, { @@ -67,21 +70,27 @@ "metadata": {}, "outputs": [], "source": [ - "results = freia_mcstas.compute((\n", - " DiskChoppers[SampleRun],\n", - " ChopperFrameSequence[SampleRun],\n", - " DetectorLtotal[SampleRun],\n", - " LookupTable[SampleRun, NXdetector],\n", - "))\n", + "results = freia_mcstas.compute(\n", + " (\n", + " DiskChoppers[SampleRun],\n", + " ChopperFrameSequence[SampleRun],\n", + " DetectorLtotal[SampleRun],\n", + " LookupTable[SampleRun, NXdetector],\n", + " )\n", + ")\n", "choppers = results[DiskChoppers[SampleRun]]\n", "frames = results[ChopperFrameSequence[SampleRun]]\n", "flight_paths = results[DetectorLtotal[SampleRun]]\n", "lookup = results[LookupTable[SampleRun, NXdetector]]\n", "\n", "print(f'Input: {filename.name}')\n", - "print(f'{len(choppers)} disks; {len(frames.frames[-1].subframes)} transmitted subframes')\n", - "print(f'Detector flight paths: {flight_paths.min().value:.4f} to {flight_paths.max().value:.4f} m')\n", - "print(f'Table shape: {lookup.array.sizes}')\n" + "print(\n", + " f'{len(choppers)} disks; {len(frames.frames[-1].subframes)} transmitted subframes'\n", + ")\n", + "print(\n", + " f'Detector flight paths: {flight_paths.min().value:.4f} to {flight_paths.max().value:.4f} m'\n", + ")\n", + "print(f'Table shape: {lookup.array.sizes}')" ] }, { @@ -91,9 +100,7 @@ "source": [ "## Wavelength as a function of arrival time and flight path\n", "\n", - "The heatmap shows the wavelength assigned to each arrival time and flight path. White regions have no modeled transmission. The lower panel shows a slice near the middle of the detector's flight-path range; shading spans the modeled wavelength bounds.\n", - "\n", - "No additional relative-uncertainty cutoff has been applied.\n" + "The heatmap shows the wavelength assigned to each arrival time and flight path. White regions have no modeled transmission. The lower panel shows a slice near the middle of the detector's flight-path range; shading spans the modeled wavelength bounds." ] }, { @@ -132,11 +139,15 @@ "figure[0, 0] = heatmap\n", "figure[1, 0] = curve\n", "figure[0, 0].cax.set_ylabel('Wavelength [Å]')\n", - "figure[0, 0].ax.axhline(line.coords['distance'].value, color='black', linestyle='--', linewidth=0.8)\n", - "figure[1, 0].ax.fill_between(time, line.values - half_width, line.values + half_width, alpha=0.3)\n", + "figure[0, 0].ax.axhline(\n", + " line.coords['distance'].value, color='black', linestyle='--', linewidth=0.8\n", + ")\n", + "figure[1, 0].ax.fill_between(\n", + " time, line.values - half_width, line.values + half_width, alpha=0.3\n", + ")\n", "figure[1, 0].ax.set_xlim(0, lookup.pulse_period.to(unit='ms').value)\n", "figure[1, 0].ax.grid(alpha=0.2)\n", - "figure\n" + "figure" ] } ], @@ -148,6 +159,9 @@ }, "language_info": { "name": "python" + }, + "nbsphinx": { + "execute": "never" } }, "nbformat": 4, diff --git a/packages/essreflectometry/docs/user-guide/freia/index.md b/packages/essreflectometry/docs/user-guide/freia/index.md index 4844732c8..fdbc1b1e0 100644 --- a/packages/essreflectometry/docs/user-guide/freia/index.md +++ b/packages/essreflectometry/docs/user-guide/freia/index.md @@ -1,7 +1,9 @@ # FREIA -Explore FREIA detector data, reconstruct wavelengths, and reduce reflectivity -using a direct-beam measurement. The examples use WFM chopper settings. +Explore FREIA detector data, reconstruct wavelengths, and compute reflectivity +using a direct-beam measurement. + +The examples use local McStas files. Set the input paths in each notebook. ```{toctree} :maxdepth: 1 From 6adc0905c4e07d5157157787724f4b95434f645d Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Thu, 17 Sep 2026 11:59:42 +0200 Subject: [PATCH 06/16] refactor: rename domain type to common ess name --- .../src/ess/reflectometry/corrections.py | 10 +++++----- .../essreflectometry/src/ess/reflectometry/types.py | 6 +++--- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/packages/essreflectometry/src/ess/reflectometry/corrections.py b/packages/essreflectometry/src/ess/reflectometry/corrections.py index ebd604ff7..ffd1ec363 100644 --- a/packages/essreflectometry/src/ess/reflectometry/corrections.py +++ b/packages/essreflectometry/src/ess/reflectometry/corrections.py @@ -8,11 +8,11 @@ from .tools import fwhm_to_std from .types import ( + CorrectedDetector, ProtonCharge, RawSampleRotation, ReducibleData, RunType, - RunUnnormalizedData, SampleRotation, SampleRotationOffset, ) @@ -85,14 +85,14 @@ def correct_by_proton_charge( def no_run_normalization( - detector: RunUnnormalizedData[RunType], + detector: CorrectedDetector[RunType], ) -> ReducibleData[RunType]: """Use prepared detector data without applying a run normalization.""" return ReducibleData[RunType](detector) def normalize_by_monitor_histogram( - detector: RunUnnormalizedData[RunType], + detector: CorrectedDetector[RunType], *, monitor: sc.DataArray, uncertainty_broadcast_mode: UncertaintyBroadcastMode, @@ -109,7 +109,7 @@ def normalize_by_monitor_histogram( def normalize_by_monitor_integrated( - detector: RunUnnormalizedData[RunType], + detector: CorrectedDetector[RunType], *, monitor: sc.DataArray, uncertainty_broadcast_mode: UncertaintyBroadcastMode, @@ -126,7 +126,7 @@ def normalize_by_monitor_integrated( def normalize_by_proton_charge( - detector: RunUnnormalizedData[RunType], + detector: CorrectedDetector[RunType], proton_charge: ProtonCharge[RunType], ) -> ReducibleData[RunType]: """Normalize detector data by time-dependent proton charge.""" diff --git a/packages/essreflectometry/src/ess/reflectometry/types.py b/packages/essreflectometry/src/ess/reflectometry/types.py index 0e8613a4e..b53e2e5d2 100644 --- a/packages/essreflectometry/src/ess/reflectometry/types.py +++ b/packages/essreflectometry/src/ess/reflectometry/types.py @@ -40,11 +40,11 @@ class RawChopper(sciline.Scope[RunType, sc.DataGroup], sc.DataGroup): class ReducibleData(sciline.Scope[RunType, sc.DataArray], sc.DataArray): - """Event data with common coordinates added""" + """Detector data ready for reduction after the selected run normalization.""" -class RunUnnormalizedData(sciline.Scope[RunType, sc.DataArray], sc.DataArray): - """Detector data prepared for reduction, before run normalization.""" +class CorrectedDetector(sciline.Scope[RunType, sc.DataArray], sc.DataArray): + """Detector data with coordinates, masks, and detector corrections applied.""" ReducedReference = NewType("ReducedReference", sc.DataArray) From 9d86aa36fd5e48cbc2a900b00a045135597ffcf8 Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Thu, 17 Sep 2026 12:00:47 +0200 Subject: [PATCH 07/16] refactor: rename in estia --- packages/essreflectometry/src/ess/estia/corrections.py | 10 +++++----- packages/essreflectometry/src/ess/estia/mcstas.py | 6 +++--- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/packages/essreflectometry/src/ess/estia/corrections.py b/packages/essreflectometry/src/ess/estia/corrections.py index ff2df30ca..1801189cc 100644 --- a/packages/essreflectometry/src/ess/estia/corrections.py +++ b/packages/essreflectometry/src/ess/estia/corrections.py @@ -9,10 +9,10 @@ from ..reflectometry.types import ( BeamDivergenceLimits, CoordTransformationGraph, + CorrectedDetector, CorrectionsToApply, ReducibleData, RunType, - RunUnnormalizedData, WavelengthBins, WavelengthDetector, YIndexLimits, @@ -24,7 +24,7 @@ def normalize_by_monitor_histogram( - detector: RunUnnormalizedData[RunType], + detector: CorrectedDetector[RunType], *, monitor: WavelengthMonitor[RunType], uncertainty_broadcast_mode: UncertaintyBroadcastMode, @@ -64,7 +64,7 @@ def normalize_by_monitor_histogram( def normalize_by_monitor_integrated( - detector: RunUnnormalizedData[RunType], + detector: CorrectedDetector[RunType], *, monitor: WavelengthMonitor[RunType], uncertainty_broadcast_mode: UncertaintyBroadcastMode, @@ -97,7 +97,7 @@ def add_coords_masks_and_apply_corrections( wbins: WavelengthBins, graph: CoordTransformationGraph[RunType], corrections_to_apply: CorrectionsToApply, -) -> RunUnnormalizedData[RunType]: +) -> CorrectedDetector[RunType]: """ Computes coordinates, masks and corrections that are the same for the sample measurement and the reference measurement. @@ -108,7 +108,7 @@ def add_coords_masks_and_apply_corrections( for correction in corrections_to_apply: da = correction(da) - return RunUnnormalizedData[RunType](da) + return CorrectedDetector[RunType](da) def correct_by_footprint(da: sc.DataArray) -> sc.DataArray: diff --git a/packages/essreflectometry/src/ess/estia/mcstas.py b/packages/essreflectometry/src/ess/estia/mcstas.py index 061a74a14..77bfd754f 100644 --- a/packages/essreflectometry/src/ess/estia/mcstas.py +++ b/packages/essreflectometry/src/ess/estia/mcstas.py @@ -12,13 +12,13 @@ from ..reflectometry.types import ( BeamDivergenceLimits, CoordTransformationGraph, + CorrectedDetector, CorrectionsToApply, DetectorLtotal, DetectorRotation, Filename, RawDetector, RunType, - RunUnnormalizedData, SampleRotation, SampleRotationOffset, WavelengthBins, @@ -294,7 +294,7 @@ def use_mcstas_wavelengths_instead_of_estimates_from_time_of_arrival( wbins: WavelengthBins, graph: CoordTransformationGraph[RunType], corrections_to_apply: CorrectionsToApply, -) -> RunUnnormalizedData[RunType]: +) -> CorrectedDetector[RunType]: out = add_coords_masks_and_apply_corrections( da=da, ylim=ylim, @@ -307,7 +307,7 @@ def use_mcstas_wavelengths_instead_of_estimates_from_time_of_arrival( }, corrections_to_apply=corrections_to_apply, ) - return RunUnnormalizedData[RunType](out) + return CorrectedDetector[RunType](out) providers = ( From 7cd3c9125ae00b84ab0a2d4aa74ab0483e2355a1 Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Thu, 17 Sep 2026 12:02:04 +0200 Subject: [PATCH 08/16] refactor: make structure more similar to other workflows --- .../freia/freia-mcstas-visualization.ipynb | 5 +- .../user-guide/freia/freia-reflectivity.ipynb | 210 +++++++++++------- .../src/ess/freia/conversions.py | 87 +++++--- .../src/ess/freia/corrections.py | 30 ++- .../src/ess/freia/maskings.py | 30 +-- .../src/ess/freia/normalization.py | 14 +- .../essreflectometry/src/ess/freia/types.py | 6 +- .../src/ess/freia/workflow.py | 7 +- .../tests/freia/conversions_test.py | 28 ++- .../tests/freia/mcstas_test.py | 12 +- .../tests/freia/workflow_test.py | 27 ++- 11 files changed, 285 insertions(+), 171 deletions(-) diff --git a/packages/essreflectometry/docs/user-guide/freia/freia-mcstas-visualization.ipynb b/packages/essreflectometry/docs/user-guide/freia/freia-mcstas-visualization.ipynb index 9932bbc57..0d958afce 100644 --- a/packages/essreflectometry/docs/user-guide/freia/freia-mcstas-visualization.ipynb +++ b/packages/essreflectometry/docs/user-guide/freia/freia-mcstas-visualization.ipynb @@ -52,7 +52,7 @@ "outputs": [], "source": [ "freia_mcstas = freia.FreiaMcStasWorkflow(wavelength_from='analytical')\n", - "freia_mcstas[Filename[SampleRun]] = '../../../265224.h5'\n", + "freia_mcstas[Filename[SampleRun]] = '../../../265305.h5'\n", "freia_mcstas[KeepEventTimeOffset] = True\n", "freia_mcstas[TimeResolution] = sc.scalar(20.0, unit='us')\n", "freia_mcstas[DistanceResolution] = sc.scalar(0.1, unit='m')" @@ -209,7 +209,8 @@ "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", - "pygments_lexer": "ipython3" + "pygments_lexer": "ipython3", + "version": "3.12.13" }, "nbsphinx": { "execute": "never" diff --git a/packages/essreflectometry/docs/user-guide/freia/freia-reflectivity.ipynb b/packages/essreflectometry/docs/user-guide/freia/freia-reflectivity.ipynb index e4cca52fc..e8bae3676 100644 --- a/packages/essreflectometry/docs/user-guide/freia/freia-reflectivity.ipynb +++ b/packages/essreflectometry/docs/user-guide/freia/freia-reflectivity.ipynb @@ -7,9 +7,9 @@ "source": [ "# FREIA reflectivity reduction\n", "\n", - "Compute reflectivity from a sample measurement and a direct-beam measurement taken without a sample. The workflow converts events to specular Q, normalizes by an incident wavelength monitor, and divides the selected peaks to obtain R(Q).\n", + "Compute one reflectivity curve for each incident beam using a sample measurement and a direct-beam measurement taken without a sample. Normalize by an incident wavelength monitor, select matching peaks, and compare the three curves with the silicon reflectivity reference.\n", "\n", - "This example uses local McStas files with WFM chopper settings. The sample and direct-beam runs must have matching slit and chopper settings." + "The sample and direct-beam runs must have matching slit and chopper settings. This example uses local files with WFM chopper settings." ] }, { @@ -20,6 +20,7 @@ "outputs": [], "source": [ "%matplotlib inline\n", + "import numpy as np\n", "import scipp as sc\n", "\n", "from ess.freia import FreiaMcStasWorkflow\n", @@ -27,21 +28,21 @@ "from ess.freia.types import (\n", " DetectorRegionOfInterest,\n", " IncidentMonitor,\n", - " QDetector,\n", " SampleIlluminatedFraction,\n", " WavelengthMonitor,\n", ")\n", - "from ess.reduce.nexus.types import GravityVector, NeXusName\n", + "from ess.reduce.nexus.types import NeXusName\n", "from ess.reflectometry.types import (\n", " BeamSize,\n", + " CorrectedDetector,\n", " Filename,\n", " QBins,\n", - " ReducibleData,\n", " ReferenceRun,\n", " ReflectivityOverQ,\n", " SampleRun,\n", " SampleSize,\n", " WavelengthBins,\n", + " WavelengthDetector,\n", ")" ] }, @@ -52,9 +53,17 @@ "source": [ "## Select the runs\n", "\n", - "Set the sample and no-sample filenames. These paths are relative to the notebook directory. Choose an incident wavelength monitor upstream of the sample, so its intensity is independent of reflectivity. If no monitor is available, use `RunNormalization.none`; the two runs must then already share an exposure and flux scale.\n", + "Paths are relative to this notebook. `SampleLambda` measures the incident spectrum before the sample in these files. If no incident monitor is available, use `RunNormalization.none`; the two runs must then already share an exposure and flux scale.\n", "\n", - "These runs were simulated without gravity, so `GravityVector` is set to zero. Use the workflow default for data that includes gravity." + "Gravity correction is enabled for these runs. For data simulated without gravity, set it to zero after creating the workflow:\n", + "\n", + "```python\n", + "from ess.reduce.nexus.types import GravityVector\n", + "\n", + "workflow[GravityVector] = sc.vector([0.0, 0.0, 0.0], unit='m/s^2')\n", + "```\n", + "\n", + "The Q grid covers all three beams. Bins without direct-beam coverage are masked separately for each curve." ] }, { @@ -65,13 +74,12 @@ "outputs": [], "source": [ "workflow = FreiaMcStasWorkflow(run_norm=RunNormalization.monitor_histogram)\n", - "workflow[Filename[SampleRun]] = '../../../265224.h5'\n", - "workflow[Filename[ReferenceRun]] = '../../../265222.h5'\n", - "workflow[NeXusName[IncidentMonitor]] = 'GuideexitLambda'\n", - "workflow[GravityVector] = sc.vector([0.0, 0.0, 0.0], unit='m/s^2')\n", + "workflow[Filename[SampleRun]] = '../../../265305.h5'\n", + "workflow[Filename[ReferenceRun]] = '../../../265301.h5'\n", + "workflow[NeXusName[IncidentMonitor]] = 'SampleLambda'\n", "\n", "workflow[WavelengthBins] = sc.linspace('wavelength', 2.0, 10.0, 81, unit='angstrom')\n", - "workflow[QBins] = sc.geomspace('Q', 0.07, 0.4, 21, unit='1/angstrom')" + "workflow[QBins] = sc.geomspace('Q', 0.003, 0.5, 151, unit='1/angstrom')" ] }, { @@ -79,11 +87,7 @@ "id": "4", "metadata": {}, "source": [ - "## Inspect the reflected and direct peaks\n", - "\n", - "The angle θ is measured above the sample surface and is corrected for gravity when enabled. Its sign distinguishes the two sides of the sample plane. The specular assumption gives $Q = 4\\pi\\sin(\\theta)/\\lambda$.\n", - "\n", - "Inspect both distributions before selecting corresponding peaks." + "Load the events, reconstruct wavelengths, and load the monitors once. Reuse these inputs when selecting each beam." ] }, { @@ -93,15 +97,16 @@ "metadata": {}, "outputs": [], "source": [ - "detectors = workflow.compute((QDetector[SampleRun], QDetector[ReferenceRun]))\n", - "angle_bins = sc.linspace('theta', -4.0, 4.0, 161, unit='deg').to(unit='rad')\n", - "profiles = {}\n", - "for label, run in [('Sample', SampleRun), ('Direct beam', ReferenceRun)]:\n", - " detector = detectors[QDetector[run]]\n", - " profile = detector.hist(theta=angle_bins, dim=detector.dims)\n", - " profile.coords['theta'] = profile.coords['theta'].to(unit='deg')\n", - " profiles[label] = profile\n", - "sc.plot(profiles, norm='log', title='Reflected and direct beams', vmin=1e4, vmax=1e8)" + "inputs = workflow.compute(\n", + " (\n", + " WavelengthDetector[SampleRun],\n", + " WavelengthDetector[ReferenceRun],\n", + " WavelengthMonitor[SampleRun],\n", + " WavelengthMonitor[ReferenceRun],\n", + " )\n", + ")\n", + "for key, value in inputs.items():\n", + " workflow[key] = value" ] }, { @@ -109,9 +114,13 @@ "id": "6", "metadata": {}, "source": [ - "## Select matching regions of interest\n", + "## Inspect the reflected and direct peaks\n", + "\n", + "The signed scattering angle γ is the elevation of the outgoing ray above the laboratory x–z plane, corrected for gravity. It is available for both runs.\n", "\n", - "Select the reflected peak and its corresponding direct peak separately. The example uses the beam near 3.5°. Additional bounds can be set on pixel coordinates such as `height` or `pixel_id`. Choose wavelength and Q bins to match the available statistics." + "For the sample, θ is the angle between the outgoing ray and the sample surface. Assuming specular reflection, this also gives the incidence angle and hence $Q = 4\\pi\\sin(\\theta)/\\lambda$. The conversion uses the full three-dimensional direction; θ = γ − μ, with sample tilt μ, applies when the ray and sample tilt lie in the vertical scattering plane.\n", + "\n", + "Inspect both distributions before selecting corresponding peaks. These runs have reflected beams near +0.3°, +1°, and +3.5°, with direct beams at the corresponding negative angles." ] }, { @@ -121,12 +130,21 @@ "metadata": {}, "outputs": [], "source": [ - "workflow[DetectorRegionOfInterest[SampleRun]] = {\n", - " 'theta': (sc.scalar(3.2, unit='deg'), sc.scalar(3.9, unit='deg')),\n", - "}\n", - "workflow[DetectorRegionOfInterest[ReferenceRun]] = {\n", - " 'theta': (sc.scalar(-3.9, unit='deg'), sc.scalar(-3.2, unit='deg')),\n", - "}" + "workflow[DetectorRegionOfInterest[SampleRun]] = {}\n", + "workflow[DetectorRegionOfInterest[ReferenceRun]] = {}\n", + "detectors = workflow.compute(\n", + " (CorrectedDetector[SampleRun], CorrectedDetector[ReferenceRun])\n", + ")\n", + "angle_bins = sc.linspace('scattering_angle', -4.2, 4.2, 421, unit='deg').to(unit='rad')\n", + "profiles = {}\n", + "for label, run in [('Sample', SampleRun), ('Direct beam', ReferenceRun)]:\n", + " detector = detectors[CorrectedDetector[run]]\n", + " profile = detector.hist(scattering_angle=angle_bins, dim=detector.dims)\n", + " profile.coords['scattering_angle'] = profile.coords['scattering_angle'].to(\n", + " unit='deg'\n", + " )\n", + " profiles[label] = profile\n", + "sc.plot(profiles, logy=True, title='Reflected and direct beams', vmin=1e0)" ] }, { @@ -136,7 +154,7 @@ "source": [ "## Inspect the wavelength normalization\n", "\n", - "The incident monitor corrects wavelength-dependent flux differences between runs. Its wavelength range must cover the detector data." + "The incident monitor corrects wavelength-dependent flux differences between runs. Its wavelength range must cover the selected detector data." ] }, { @@ -146,99 +164,134 @@ "metadata": {}, "outputs": [], "source": [ - "monitors = workflow.compute(\n", - " (WavelengthMonitor[SampleRun], WavelengthMonitor[ReferenceRun])\n", - ")\n", "sc.plot(\n", " {\n", - " 'Sample monitor': monitors[WavelengthMonitor[SampleRun]],\n", - " 'Direct-beam monitor': monitors[WavelengthMonitor[ReferenceRun]],\n", + " 'Sample monitor': inputs[WavelengthMonitor[SampleRun]],\n", + " 'Direct-beam monitor': inputs[WavelengthMonitor[ReferenceRun]],\n", " },\n", " title='Incident wavelength spectra',\n", ")" ] }, + { + "cell_type": "markdown", + "id": "10", + "metadata": {}, + "source": [ + "## Select the beams and configure footprint correction\n", + "\n", + "Each entry below gives separate scattering-angle bounds for the reflected and direct peaks, in degrees. Adjust them after inspecting the profiles when using other runs.\n", + "\n", + "The Gaussian footprint correction uses the sample length along the beam and each beam's FWHM at the sample. It applies only to the reflected run. Beam widths have not yet been established for these data, so the example leaves this correction off. To enable it for a beam, replace its `None` width with a measured value, for example `sc.scalar(width_in_mm, unit='mm')`." + ] + }, { "cell_type": "code", "execution_count": null, - "id": "10", + "id": "11", "metadata": {}, "outputs": [], "source": [ - "selected = workflow.compute((ReducibleData[SampleRun], ReducibleData[ReferenceRun]))\n", - "spectra = {}\n", - "for label, run in [('Sample', SampleRun), ('Direct beam', ReferenceRun)]:\n", - " detector = selected[ReducibleData[run]]\n", - " spectra[label] = detector.hist(\n", - " wavelength=workflow.compute(WavelengthBins), dim=detector.dims\n", - " )\n", - "sc.plot(spectra, title='Selected peaks after monitor normalization')" + "# Beam label: (sample ROI, direct-beam ROI), in degrees.\n", + "beam_rois = {\n", + " '0.3° beam': ((0.2, 0.4), (-0.4, -0.2)),\n", + " '1° beam': ((0.8, 1.2), (-1.2, -0.8)),\n", + " '3.5° beam': ((3.1, 3.9), (-3.9, -3.1)),\n", + "}\n", + "sample_size = sc.scalar(80.0, unit='mm')\n", + "beam_sizes = dict.fromkeys(beam_rois)" ] }, { "cell_type": "markdown", - "id": "11", + "id": "12", "metadata": {}, "source": [ - "## Configure footprint correction\n", + "## Compute one reflectivity curve per beam\n", "\n", - "The Gaussian footprint model uses the sample length along the beam and the beam FWHM at the sample. It corrects only the reflected run. Different incident beams may require different widths.\n", + "For each beam, apply the matching sample and direct-beam ROIs to a copy of the workflow. When building the reference, the workflow reflects the direct-beam direction in the sample plane and uses the resulting reflection angle to compute Q. It integrates both selected peaks into matching Q bins and divides their intensities, propagating counting uncertainties.\n", "\n", - "A beam width has not yet been established for these data, so this example leaves the footprint correction off. Replace `beam_size = None` with a measured width, for example `sc.scalar(width_in_mm, unit='mm')`, to enable it." + "The curves retain their direct-beam normalization; no scale is fitted to the silicon reference or to another beam." ] }, { "cell_type": "code", "execution_count": null, - "id": "12", + "id": "13", "metadata": {}, "outputs": [], "source": [ - "sample_size = sc.scalar(80.0, unit='mm')\n", - "beam_size = None\n", - "\n", - "reduction = workflow.copy()\n", - "if beam_size is None:\n", - " reduction[SampleIlluminatedFraction] = sc.scalar(1.0)\n", - "else:\n", - " reduction[SampleSize[SampleRun]] = sample_size\n", - " reduction[BeamSize[SampleRun]] = beam_size" + "reflectivities = {}\n", + "for label, (sample_roi, direct_roi) in beam_rois.items():\n", + " reduction = workflow.copy()\n", + " for run, bounds in [(SampleRun, sample_roi), (ReferenceRun, direct_roi)]:\n", + " reduction[DetectorRegionOfInterest[run]] = {\n", + " 'scattering_angle': tuple(sc.scalar(edge, unit='deg') for edge in bounds),\n", + " }\n", + " beam_size = beam_sizes[label]\n", + " if beam_size is None:\n", + " reduction[SampleIlluminatedFraction] = sc.scalar(1.0)\n", + " else:\n", + " reduction[SampleSize[SampleRun]] = sample_size\n", + " reduction[BeamSize[SampleRun]] = beam_size\n", + " reflectivity = reduction.compute(ReflectivityOverQ)\n", + " reflectivities[label] = reflectivity\n", + " covered = (~reflectivity.masks['direct_beam']).sum().value\n", + " print(f'{label}: {covered} Q bins have direct-beam coverage.')" ] }, { "cell_type": "markdown", - "id": "13", + "id": "14", "metadata": {}, "source": [ - "## Compute reflectivity\n", + "## Compare with the silicon reference\n", "\n", - "The reference uses the negative of the direct beam’s θ to compute its corresponding reflected Q. The workflow integrates the two selected peaks into matching Q bins and divides their intensities. Both counting uncertainties propagate. Bins with no usable direct-beam intensity are masked." + "`Si-15SiO2-air.txt` is the exact reflectivity table used for the silicon sample in these runs. Its columns are Q in Å⁻¹ and dimensionless reflectivity. Plot it together with all three reduced curves on the same axes." ] }, { "cell_type": "code", "execution_count": null, - "id": "14", + "id": "15", "metadata": {}, "outputs": [], "source": [ - "reflectivity = reduction.compute(ReflectivityOverQ)\n", - "covered = (~reflectivity.masks['direct_beam']).sum().value\n", - "print(f'{covered} of {reflectivity.sizes[\"Q\"]} Q bins have direct-beam coverage.')\n", - "title = (\n", - " 'Reflectivity'\n", - " if beam_size is not None\n", - " else 'Reflectivity (footprint correction omitted)'\n", + "reference_file = '../../../../../Si-15SiO2-air.txt'\n", + "q, r = np.loadtxt(reference_file, unpack=True)\n", + "si_reference = sc.DataArray(\n", + " sc.array(dims=['Q'], values=r, unit='dimensionless'),\n", + " coords={'Q': sc.array(dims=['Q'], values=q, unit='1/angstrom')},\n", ")\n", - "reflectivity.plot(norm='log', title=title)" + "# Display each histogram value at its Q-bin center.\n", + "plot_curves = {\n", + " label: curve.assign_coords(Q=sc.midpoints(curve.coords['Q']))\n", + " for label, curve in reflectivities.items()\n", + "}\n", + "qbins = workflow.compute(QBins)\n", + "comparison = sc.plot(\n", + " {'Si reference': si_reference, **plot_curves},\n", + " logx=True,\n", + " logy=True,\n", + " linestyle={'Si reference': '-'},\n", + " marker={'Si reference': 'none'},\n", + " color={'Si reference': 'black'},\n", + " markersize=4,\n", + " figsize=(8, 5),\n", + " xmin=qbins[0],\n", + " xmax=qbins[-1],\n", + " title='FREIA reflectivity by beam',\n", + " ylabel='Reflectivity',\n", + ")\n", + "comparison" ] }, { "cell_type": "markdown", - "id": "15", + "id": "16", "metadata": {}, "source": [ - "Footprint correction is disabled in this example. Q resolution, background subtraction, finite-sample corrections, and complete ORSO export are not yet included." + "Footprint correction is disabled with the widths above. Background subtraction and Q-resolution averaging are not included, and the reference is shown without resolution broadening. These effects and the counting statistics should be considered when comparing the curves." ] } ], @@ -257,7 +310,8 @@ "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", - "pygments_lexer": "ipython3" + "pygments_lexer": "ipython3", + "version": "3.12.13" }, "nbsphinx": { "execute": "never" diff --git a/packages/essreflectometry/src/ess/freia/conversions.py b/packages/essreflectometry/src/ess/freia/conversions.py index 4a6a04239..70e6e5a53 100644 --- a/packages/essreflectometry/src/ess/freia/conversions.py +++ b/packages/essreflectometry/src/ess/freia/conversions.py @@ -6,21 +6,23 @@ from scippnexus import NXsample, NXsource from ..reflectometry.conversions import reflectometry_q -from ..reflectometry.types import CoordTransformationGraph, RunType, WavelengthDetector -from .types import QDetector, SampleSurfaceNormal +from ..reflectometry.types import ( + CoordTransformationGraph, + RunType, + SampleRun, +) +from .types import SampleSurfaceNormal -def theta( +def outgoing_direction( scattered_beam: sc.Variable, wavelength: sc.Variable, gravity: sc.Variable, - sample_surface_normal: sc.Variable, ) -> sc.Variable: - """Signed, gravity-corrected angle above the sample plane. + """Unit direction of the outgoing ray at the sample, corrected for gravity. Approximate the flight time using the straight sample-to-detector distance, - as in ScippNeutron's gravity correction. Positive angles point toward the - sample surface normal. + as in ScippNeutron's gravity correction. """ flight_time = tof.tof_from_wavelength( wavelength=wavelength, Ltotal=sc.norm(scattered_beam) @@ -28,8 +30,25 @@ def theta( outgoing_beam = scattered_beam - (0.5 * gravity * flight_time**2).to( unit=scattered_beam.unit ) + return outgoing_beam / sc.norm(outgoing_beam) + + +def scattering_angle(outgoing_direction: sc.Variable) -> sc.Variable: + """Signed elevation of the outgoing ray above the laboratory x-z plane.""" + return sc.asin(outgoing_direction.fields.y) + + +def theta( + outgoing_direction: sc.Variable, + sample_surface_normal: sc.Variable, +) -> sc.Variable: + """Specular reflection angle between the outgoing ray and the sample surface. + + Use the full three-dimensional direction. Under specular reflection this + also determines the incidence angle, without requiring the incoming ray. + """ normal = sample_surface_normal / sc.norm(sample_surface_normal) - return sc.asin(sc.dot(outgoing_beam, normal) / sc.norm(outgoing_beam)) + return sc.asin(sc.dot(outgoing_direction, normal)) def coordinate_transformation_graph( @@ -38,7 +57,7 @@ def coordinate_transformation_graph( sample_surface_normal: SampleSurfaceNormal[RunType], gravity: GravityVector, ) -> CoordTransformationGraph[RunType]: - """Build a specular conversion graph.""" + """Build the scattering coordinates shared by sample and direct-beam runs.""" length = sc.norm(sample_surface_normal) if ( not sc.isfinite(length).value @@ -46,9 +65,10 @@ def coordinate_transformation_graph( ): raise ValueError('SampleSurfaceNormal must be a finite, nonzero vector.') return { - **graph.beamline.beamline(scatter=True), - 'theta': theta, - 'Q': reflectometry_q, + **graph.beamline.L1(), + **graph.beamline.L2(), + 'outgoing_direction': outgoing_direction, + 'scattering_angle': scattering_angle, 'source_position': lambda: source_position, 'sample_position': lambda: sample_position, 'sample_surface_normal': lambda: sample_surface_normal, @@ -56,28 +76,27 @@ def coordinate_transformation_graph( } +def sample_coordinate_transformation_graph( + source_position: Position[NXsource, SampleRun], + sample_position: Position[NXsample, SampleRun], + sample_surface_normal: SampleSurfaceNormal[SampleRun], + gravity: GravityVector, +) -> CoordTransformationGraph[SampleRun]: + """Extend the scattering graph with the sample's reflection angle and Q.""" + return coordinate_transformation_graph( + source_position, sample_position, sample_surface_normal, gravity + ) | {'theta': theta, 'Q': reflectometry_q} + + def add_coords( - da: WavelengthDetector[RunType], - graph: CoordTransformationGraph[RunType], -) -> QDetector[RunType]: - """Add specular Q and the gravity-corrected angle to detector events.""" - return QDetector[RunType]( - da.transform_coords( - ( - 'theta', - 'Q', - 'L1', - 'L2', - 'incident_beam', - 'sample_position', - 'sample_surface_normal', - ), - graph, - rename_dims=False, - keep_intermediate=False, - keep_aliases=False, - ) - ) + da: sc.DataArray, + graph: dict, +) -> sc.DataArray: + """Add the scattering coordinates provided by the run's transformation graph.""" + return da.transform_coords(rename_dims=False, **graph) -providers = (coordinate_transformation_graph, add_coords) +providers = ( + coordinate_transformation_graph, + sample_coordinate_transformation_graph, +) diff --git a/packages/essreflectometry/src/ess/freia/corrections.py b/packages/essreflectometry/src/ess/freia/corrections.py index 0cf9885a1..4143a5bbd 100644 --- a/packages/essreflectometry/src/ess/freia/corrections.py +++ b/packages/essreflectometry/src/ess/freia/corrections.py @@ -8,18 +8,38 @@ from ..reflectometry.corrections import RunNormalization from ..reflectometry.types import ( BeamSize, + CoordTransformationGraph, + CorrectedDetector, ReducibleData, RunType, - RunUnnormalizedData, Sample, SampleRun, SampleSize, + WavelengthBins, + WavelengthDetector, ) -from .types import SampleIlluminatedFraction, WavelengthMonitor +from .conversions import add_coords +from .maskings import add_masks +from .types import ( + DetectorRegionOfInterest, + SampleIlluminatedFraction, + WavelengthMonitor, +) + + +def add_coords_and_masks( + da: WavelengthDetector[RunType], + graph: CoordTransformationGraph[RunType], + roi: DetectorRegionOfInterest[RunType], + wavelength_bins: WavelengthBins, +) -> CorrectedDetector[RunType]: + """Transform coordinates and mask events before run normalization.""" + da = add_coords(da, graph) + return CorrectedDetector[RunType](add_masks(da, roi, wavelength_bins)) def normalize_by_monitor_histogram( - detector: RunUnnormalizedData[RunType], + detector: CorrectedDetector[RunType], *, monitor: WavelengthMonitor[RunType], uncertainty_broadcast_mode: UncertaintyBroadcastMode, @@ -59,7 +79,7 @@ def normalize_by_monitor_histogram( def normalize_by_monitor_integrated( - detector: RunUnnormalizedData[RunType], + detector: CorrectedDetector[RunType], *, monitor: WavelengthMonitor[RunType], uncertainty_broadcast_mode: UncertaintyBroadcastMode, @@ -122,4 +142,4 @@ def prepare_sample( return Sample(sample / fraction) -providers = (sample_illuminated_fraction, prepare_sample) +providers = (add_coords_and_masks, sample_illuminated_fraction, prepare_sample) diff --git a/packages/essreflectometry/src/ess/freia/maskings.py b/packages/essreflectometry/src/ess/freia/maskings.py index 84eec42a3..a0a60483d 100644 --- a/packages/essreflectometry/src/ess/freia/maskings.py +++ b/packages/essreflectometry/src/ess/freia/maskings.py @@ -2,15 +2,12 @@ # Copyright (c) 2026 Scipp contributors (https://github.com/scipp) import scipp as sc -from ..reflectometry.types import RunType, RunUnnormalizedData, WavelengthBins -from .types import DetectorRegionOfInterest, QDetector - def add_masks( - da: QDetector[RunType], - roi: DetectorRegionOfInterest[RunType], - wavelength_bins: WavelengthBins, -) -> RunUnnormalizedData[RunType]: + da: sc.DataArray, + roi: dict, + wavelength_bins: sc.Variable, +) -> sc.DataArray: """Mask events outside the ROI and wavelength range.""" masks = {} event_masks = {} @@ -24,16 +21,11 @@ def add_masks( (coord >= low) & (coord <= high) ) wavelength = da.bins.coords['wavelength'] - return RunUnnormalizedData[RunType]( - da.assign_masks(masks).bins.assign_masks( - event_masks, - wavelength=~( - sc.isfinite(wavelength) - & (wavelength >= wavelength_bins[0].to(unit=wavelength.unit)) - & (wavelength < wavelength_bins[-1].to(unit=wavelength.unit)) - ), - ) + return da.assign_masks(masks).bins.assign_masks( + event_masks, + wavelength=~( + sc.isfinite(wavelength) + & (wavelength >= wavelength_bins[0].to(unit=wavelength.unit)) + & (wavelength < wavelength_bins[-1].to(unit=wavelength.unit)) + ), ) - - -providers = (add_masks,) diff --git a/packages/essreflectometry/src/ess/freia/normalization.py b/packages/essreflectometry/src/ess/freia/normalization.py index 9aeb35326..2faf0a962 100644 --- a/packages/essreflectometry/src/ess/freia/normalization.py +++ b/packages/essreflectometry/src/ess/freia/normalization.py @@ -11,16 +11,24 @@ ReflectivityOverQ, Sample, ) +from .conversions import theta def evaluate_direct_beam( direct_beam: ReducibleData[ReferenceRun], ) -> Reference: - """Compute reference Q using the direct beam's incidence angle.""" - theta = -direct_beam.bins.coords['theta'] + """Compute reference Q for the direct beam reflected in the sample plane.""" + normal = direct_beam.coords['sample_surface_normal'] + normal = normal / sc.norm(normal) + outgoing = direct_beam.bins.coords['outgoing_direction'] + reflected = outgoing - 2 * sc.dot(outgoing, normal) * normal + reflection_angle = theta( + outgoing_direction=reflected, + sample_surface_normal=normal, + ) wavelength = direct_beam.bins.coords['wavelength'] return Reference( - direct_beam.bins.assign_coords(Q=reflectometry_q(wavelength, theta)) + direct_beam.bins.assign_coords(Q=reflectometry_q(wavelength, reflection_angle)) ) diff --git a/packages/essreflectometry/src/ess/freia/types.py b/packages/essreflectometry/src/ess/freia/types.py index c5e7a401a..42e96d6d2 100644 --- a/packages/essreflectometry/src/ess/freia/types.py +++ b/packages/essreflectometry/src/ess/freia/types.py @@ -23,15 +23,11 @@ class SampleSurfaceNormal(sciline.Scope[RunType, sc.Variable], sc.Variable): """Normal pointing out of the reflecting surface, in global coordinates.""" -class QDetector(sciline.Scope[RunType, sc.DataArray], sc.DataArray): - """Detector events with specular Q and signed angle to the sample surface.""" - - class DetectorRegionOfInterest(sciline.Scope[RunType, dict], dict): """Pixel or event coordinates mapped to inclusive (lower, upper) bounds. Select corresponding reflected and direct peaks separately for SampleRun and - ReferenceRun, for example using signed ``theta`` and ``height``. + ReferenceRun, for example using ``scattering_angle`` and ``height``. An empty dictionary explicitly selects the entire detector. """ diff --git a/packages/essreflectometry/src/ess/freia/workflow.py b/packages/essreflectometry/src/ess/freia/workflow.py index c2562252b..371d2eef9 100644 --- a/packages/essreflectometry/src/ess/freia/workflow.py +++ b/packages/essreflectometry/src/ess/freia/workflow.py @@ -20,7 +20,6 @@ conversions, corrections, load, - maskings, mcstas, normalization, orso, @@ -31,7 +30,6 @@ *reflectometry_providers, *conversions.providers, *corrections.providers, - *maskings.providers, *normalization.providers, *orso.providers, *load.providers, @@ -104,7 +102,10 @@ def FreiaWorkflow( ) -> sciline.Pipeline: """Workflow for reduction of data for the Freia instrument. - ``QDetector`` provides specular Q with gravity correction. Reflectivity + The coordinate transformation graph computes the signed, gravity-corrected + scattering angle above the laboratory x-z plane for both runs, with reflection + angle and Q for the sample. The direct beam is mapped to Q when building + ``Reference``. Reflectivity requires separate sample/direct-beam ROIs, wavelength and Q bins, and beam and sample sizes for the footprint correction. The reference run must be a measurement without a sample, taken with matching slit and chopper settings. diff --git a/packages/essreflectometry/tests/freia/conversions_test.py b/packages/essreflectometry/tests/freia/conversions_test.py index e61d0909c..094af42bf 100644 --- a/packages/essreflectometry/tests/freia/conversions_test.py +++ b/packages/essreflectometry/tests/freia/conversions_test.py @@ -4,33 +4,40 @@ import scipp as sc from scipp.testing import assert_allclose -from ess.freia.conversions import theta +from ess.freia.conversions import outgoing_direction, scattering_angle, theta -def test_theta_relative_to_sample_surface(): - result = theta( +def test_scattering_and_reflection_angles(): + normal = sc.vector([0.0, np.sqrt(3) / 2, -0.5]) + # The first ray is 15 degrees above the tilted surface. The second lies in + # the surface but outside the vertical scattering plane. + direction = outgoing_direction( scattered_beam=sc.vectors( - dims=['event'], values=[[0, 1, 1], [0, -1, 1]], unit='m' + dims=['event'], values=[[0, 1, 1], [1, 0.5, np.sqrt(3) / 2]], unit='m' ), wavelength=sc.scalar(4.0, unit='angstrom'), gravity=sc.vector([0.0, 0.0, 0.0], unit='m/s^2'), - sample_surface_normal=sc.vector([0.0, np.sqrt(3) / 2, -0.5]), ) assert_allclose( - result, - sc.array(dims=['event'], values=[15.0, -75.0], unit='deg').to(unit='rad'), + scattering_angle(direction)['event', 0], + sc.scalar(45.0, unit='deg').to(unit='rad'), + ) + assert_allclose( + theta(direction, normal), + sc.array(dims=['event'], values=[15.0, 0.0], unit='deg').to(unit='rad'), + atol=sc.scalar(1e-15, unit='rad'), ) -def test_theta_with_gravity(): +def test_scattering_angle_with_gravity(): # Horizontal rays at 1000 and 500 m/s fall 0.49 and 1.96 mm over 10 m. wavelength = ( sc.constants.h / sc.constants.m_n / sc.array(dims=['event'], values=[1000.0, 500.0], unit='m/s') ) - result = theta( + direction = outgoing_direction( scattered_beam=sc.vectors( dims=['event'], values=[[0, -0.0004903325, 10], [0, -0.00196133, 10]], @@ -38,11 +45,10 @@ def test_theta_with_gravity(): ), wavelength=wavelength.to(unit='angstrom'), gravity=sc.vector([0.0, -9.80665, 0.0], unit='m/s^2'), - sample_surface_normal=sc.vector([0.0, 1.0, 0.0]), ) assert_allclose( - result, + scattering_angle(direction), sc.zeros(dims=['event'], shape=[2], unit='rad'), atol=sc.scalar(1e-10, unit='rad'), ) diff --git a/packages/essreflectometry/tests/freia/mcstas_test.py b/packages/essreflectometry/tests/freia/mcstas_test.py index d83073b74..c651aba5b 100644 --- a/packages/essreflectometry/tests/freia/mcstas_test.py +++ b/packages/essreflectometry/tests/freia/mcstas_test.py @@ -11,8 +11,8 @@ from ess.freia import FreiaMcStasWorkflow from ess.freia.mcstas import mcstas_detector_geometry -from ess.freia.types import IncidentMonitor, QDetector, WavelengthMonitor -from ess.reflectometry.types import SampleRun +from ess.freia.types import DetectorRegionOfInterest, IncidentMonitor, WavelengthMonitor +from ess.reflectometry.types import CorrectedDetector, SampleRun, WavelengthBins def _component(components, name, position): @@ -65,8 +65,12 @@ def mcstas_file(tmp_path: Path) -> Path: def test_load_detector_and_compute_q(mcstas_file): workflow = FreiaMcStasWorkflow() workflow[Filename[SampleRun]] = mcstas_file + workflow[DetectorRegionOfInterest[SampleRun]] = {} + workflow[WavelengthBins] = sc.array( + dims=['wavelength'], values=[1.0, 12.0], unit='angstrom' + ) - result = workflow.compute((RawDetector[SampleRun], QDetector[SampleRun])) + result = workflow.compute((RawDetector[SampleRun], CorrectedDetector[SampleRun])) raw = result[RawDetector[SampleRun]] image = sc.sort(raw.bins.sum(), 'pixel_id') @@ -88,7 +92,7 @@ def test_load_detector_and_compute_q(mcstas_file): events.coords['event_time_offset'], sc.full(sizes=events.sizes, value=0.025, unit='s'), ) - q = result[QDetector[SampleRun]].bins.constituents['data'].coords['Q'] + q = result[CorrectedDetector[SampleRun]].bins.constituents['data'].coords['Q'] assert sc.isfinite(q).all().value diff --git a/packages/essreflectometry/tests/freia/workflow_test.py b/packages/essreflectometry/tests/freia/workflow_test.py index 319ab8ce5..71543bd06 100644 --- a/packages/essreflectometry/tests/freia/workflow_test.py +++ b/packages/essreflectometry/tests/freia/workflow_test.py @@ -18,6 +18,7 @@ from ess.reflectometry.types import ( BeamSize, QBins, + ReducibleData, ReferenceRun, ReflectivityOverQ, SampleRun, @@ -43,19 +44,27 @@ def _make_workflow(run_norm): 'pixel_id': sc.array(dims=['event'], values=[0, 0, 1], unit=None), }, ).group('pixel_id') - # The first pixel is at 30 degrees; the second is outside the ROI. + # Sample rotation is 10 degrees. The matching beams make angles of + # 30 degrees with the surface; the second pixel is outside the ROI. + angles = np.deg2rad([10.0 + sign * 30.0, 10.0 + sign * 60.0]) events.coords['position'] = sc.vectors( dims=['pixel_id'], - values=[[0, sign * 0.5, np.sqrt(3) / 2], [0, sign * np.sqrt(3) / 2, 0.5]], + values=[[0, np.sin(angle), np.cos(angle)] for angle in angles], unit='m', ) wf[WavelengthDetector[run]] = events wf[Position[NXsample, run]] = sc.vector([0.0, 0.0, 0.0], unit='m') - wf[Position[NXsource, run]] = sc.vector([0.0, 0.0, -20.0], unit='m') - wf[SampleSurfaceNormal[run]] = sc.vector([0.0, 1.0, 0.0]) - low, high = sorted([sign * 20.0, sign * 40.0]) + # Source position must not set the angular reference frame. + wf[Position[NXsource, run]] = sc.vector([3.0, 2.0, -20.0], unit='m') + wf[SampleSurfaceNormal[run]] = sc.vector( + [0.0, np.cos(np.deg2rad(10.0)), -np.sin(np.deg2rad(10.0))] + ) + low, high = sorted([10.0 + sign * 20.0, 10.0 + sign * 40.0]) wf[DetectorRegionOfInterest[run]] = { - 'theta': (sc.scalar(low, unit='deg'), sc.scalar(high, unit='deg')) + 'scattering_angle': ( + sc.scalar(low, unit='deg'), + sc.scalar(high, unit='deg'), + ) } wf[WavelengthBins] = sc.array( dims=['wavelength'], values=[1.0, 3.0, 5.0], unit='angstrom' @@ -76,7 +85,11 @@ def test_reduce_reflectivity_with_monitor_and_footprint(): coords={'wavelength': wf.compute(WavelengthBins)}, ) - result = wf.compute(ReflectivityOverQ) + results = wf.compute((ReflectivityOverQ, ReducibleData[ReferenceRun])) + result = results[ReflectivityOverQ] + direct_beam = results[ReducibleData[ReferenceRun]] + assert 'theta' not in direct_beam.bins.coords + assert 'Q' not in direct_beam.bins.coords fraction = footprint_on_sample(sc.scalar(30.0, unit='deg'), beam_size, sample_size) # Q bins contain wavelengths 4 and 2 angstrom, respectively. From 8302aadbe6ebdb452a11c5994b31660bbf36e82e Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Thu, 17 Sep 2026 14:36:43 +0200 Subject: [PATCH 09/16] refactor: remove unnecessary loaders and domain types --- .../src/ess/freia/__init__.py | 2 - .../src/ess/freia/beamline.py | 28 -------- .../src/ess/freia/corrections.py | 68 ++++--------------- .../essreflectometry/src/ess/freia/load.py | 27 -------- .../essreflectometry/src/ess/freia/types.py | 4 -- .../src/ess/freia/workflow.py | 32 +++++---- .../tests/freia/workflow_test.py | 4 +- 7 files changed, 34 insertions(+), 131 deletions(-) delete mode 100644 packages/essreflectometry/src/ess/freia/beamline.py delete mode 100644 packages/essreflectometry/src/ess/freia/load.py diff --git a/packages/essreflectometry/src/ess/freia/__init__.py b/packages/essreflectometry/src/ess/freia/__init__.py index 7e64e627d..ef3b9237b 100644 --- a/packages/essreflectometry/src/ess/freia/__init__.py +++ b/packages/essreflectometry/src/ess/freia/__init__.py @@ -5,7 +5,6 @@ from ..reflectometry import supermirror from . import ( conversions, - load, maskings, mcstas, normalization, @@ -51,7 +50,6 @@ "SampleSizeResolution", "WavelengthResolution", "conversions", - "load", "maskings", "mcstas", "normalization", diff --git a/packages/essreflectometry/src/ess/freia/beamline.py b/packages/essreflectometry/src/ess/freia/beamline.py deleted file mode 100644 index fc87b9c4e..000000000 --- a/packages/essreflectometry/src/ess/freia/beamline.py +++ /dev/null @@ -1,28 +0,0 @@ -import sciline -from ess.reduce.nexus.types import DetectorBankSizes -from ess.reduce.unwrap.workflow import GenericUnwrapWorkflow - -from ess.reflectometry.types import ReferenceRun, SampleRun - -from .types import IncidentMonitor - -DETECTOR_BANK_SIZES = { - "multiblade_detector": { - "strip": 64, - "blade": 32, - "wire": 32, - }, -} - - -def LoadNeXusWorkflow(**kwargs) -> sciline.Pipeline: - """ - Workflow for loading NeXus data. - """ - workflow = GenericUnwrapWorkflow( - run_types=[SampleRun, ReferenceRun], - monitor_types=[IncidentMonitor], - **kwargs, - ) - workflow[DetectorBankSizes] = DETECTOR_BANK_SIZES - return workflow diff --git a/packages/essreflectometry/src/ess/freia/corrections.py b/packages/essreflectometry/src/ess/freia/corrections.py index 4143a5bbd..c9ae62231 100644 --- a/packages/essreflectometry/src/ess/freia/corrections.py +++ b/packages/essreflectometry/src/ess/freia/corrections.py @@ -20,11 +20,7 @@ ) from .conversions import add_coords from .maskings import add_masks -from .types import ( - DetectorRegionOfInterest, - SampleIlluminatedFraction, - WavelengthMonitor, -) +from .types import DetectorRegionOfInterest, WavelengthMonitor def add_coords_and_masks( @@ -46,30 +42,8 @@ def normalize_by_monitor_histogram( ) -> ReducibleData[RunType]: """Normalize detector data by a histogrammed monitor. - The detector is normalized according to - - .. math:: - - d_i^\\text{Norm} = \\frac{d_i}{m_i} \\Delta \\lambda_i - - Parameters - ---------- - detector: - Input event data in wavelength. - monitor: - A histogrammed monitor in wavelength. - uncertainty_broadcast_mode: - Choose how uncertainties of the monitor are broadcast to the sample data. - - Returns - ------- - : - `detector` normalized by a monitor. - - See also - -------- - ess.reduce.normalization.normalize_by_monitor_histogram: - For details and the actual implementation. + See :func:`ess.reduce.normalization.normalize_by_monitor_histogram` for the + normalization and uncertainty treatment. """ return common_corrections.normalize_by_monitor_histogram( detector=detector, @@ -104,16 +78,15 @@ def insert_run_normalization( ) -def sample_illuminated_fraction( +def prepare_sample( sample: ReducibleData[SampleRun], beam_size: BeamSize[SampleRun], sample_size: SampleSize[SampleRun], -) -> SampleIlluminatedFraction: - """Use Amor's Gaussian footprint model with the specular incidence angle. +) -> Sample: + """Apply the Gaussian footprint correction using the specular incidence angle. - Beam size is the FWHM at the sample; sample size is its length along the - beam. These must be supplied explicitly, since slit openings alone do not - determine the profile of the beam reaching the sample. + Beam size is the FWHM at the sample; sample size is its length along the beam. + The correction applies only to the reflected run. """ for name, size in [('BeamSize', beam_size), ('SampleSize', sample_size)]: if ( @@ -121,25 +94,12 @@ def sample_illuminated_fraction( or not (size > sc.scalar(0.0, unit=size.unit)).value ): raise ValueError(f'{name} must be finite and positive.') - return SampleIlluminatedFraction( - common_corrections.footprint_on_sample( - sample.bins.coords['theta'], beam_size=beam_size, sample_size=sample_size - ) + fraction = common_corrections.footprint_on_sample( + sample.bins.coords['theta'], beam_size=beam_size, sample_size=sample_size ) + corrected = sample / fraction + invalid = ~sc.isfinite(fraction) | (fraction <= sc.scalar(0.0)) + return Sample(corrected.bins.assign_masks(footprint=invalid)) -def prepare_sample( - sample: ReducibleData[SampleRun], - illuminated_fraction: SampleIlluminatedFraction, -) -> Sample: - """Correct the reflected beam for footprint; the direct beam has no sample.""" - if illuminated_fraction.bins is None: - illuminated_fraction = sc.bins_like(sample, illuminated_fraction) - valid = sc.isfinite(illuminated_fraction) & (illuminated_fraction > sc.scalar(0.0)) - valid &= illuminated_fraction <= sc.scalar(1.0) - sample = sample.bins.assign_masks(footprint=~valid) - fraction = sc.where(valid, illuminated_fraction, sc.scalar(1.0)) - return Sample(sample / fraction) - - -providers = (add_coords_and_masks, sample_illuminated_fraction, prepare_sample) +providers = (add_coords_and_masks, prepare_sample) diff --git a/packages/essreflectometry/src/ess/freia/load.py b/packages/essreflectometry/src/ess/freia/load.py deleted file mode 100644 index 51c851ae9..000000000 --- a/packages/essreflectometry/src/ess/freia/load.py +++ /dev/null @@ -1,27 +0,0 @@ -# SPDX-License-Identifier: BSD-3-Clause -# Copyright (c) 2025 Scipp contributors (https://github.com/scipp) -from ess.reduce.nexus.types import NeXusComponent -from scippnexus import NXdetector, NXsample - -from ..reflectometry.types import ( - DetectorRotation, - RawSampleRotation, - RunType, -) - - -def load_sample_rotation( - sample: NeXusComponent[NXsample, RunType], -) -> RawSampleRotation[RunType]: - """Load sample rotation from the NeXus sample group.""" - return sample['sample_rotation'][0].data - - -def load_detector_rotation( - detector: NeXusComponent[NXdetector, RunType], -) -> DetectorRotation[RunType]: - """Load detector rotation from the NeXus detector group.""" - return detector['transformations']['detector_rotation'].value[0].data - - -providers = (load_sample_rotation, load_detector_rotation) diff --git a/packages/essreflectometry/src/ess/freia/types.py b/packages/essreflectometry/src/ess/freia/types.py index 42e96d6d2..cf49d7534 100644 --- a/packages/essreflectometry/src/ess/freia/types.py +++ b/packages/essreflectometry/src/ess/freia/types.py @@ -30,7 +30,3 @@ class DetectorRegionOfInterest(sciline.Scope[RunType, dict], dict): ReferenceRun, for example using ``scattering_angle`` and ``height``. An empty dictionary explicitly selects the entire detector. """ - - -SampleIlluminatedFraction = NewType('SampleIlluminatedFraction', sc.Variable) -"""Fraction of the incoming beam hitting the sample; set to 1 to skip footprint.""" diff --git a/packages/essreflectometry/src/ess/freia/workflow.py b/packages/essreflectometry/src/ess/freia/workflow.py index 371d2eef9..eec0babd1 100644 --- a/packages/essreflectometry/src/ess/freia/workflow.py +++ b/packages/essreflectometry/src/ess/freia/workflow.py @@ -3,8 +3,10 @@ import sciline import scipp as sc +from ess.reduce.nexus.types import DetectorBankSizes from ess.reduce.uncertainty import UncertaintyBroadcastMode from ess.reduce.unwrap import WavelengthLutMode +from ess.reduce.unwrap.workflow import GenericUnwrapWorkflow from ess.reduce.workflow import register_workflow from ..reflectometry import providers as reflectometry_providers @@ -12,19 +14,12 @@ DetectorSpatialResolution, LookupTableRelativeErrorThreshold, NeXusDetectorName, - RunType, - SampleRotationOffset, -) -from . import ( - beamline, - conversions, - corrections, - load, - mcstas, - normalization, - orso, + ReferenceRun, + SampleRun, ) +from . import conversions, corrections, mcstas, normalization, orso from .corrections import RunNormalization, insert_run_normalization +from .types import IncidentMonitor providers = ( *reflectometry_providers, @@ -32,7 +27,6 @@ *corrections.providers, *normalization.providers, *orso.providers, - *load.providers, ) """List of providers for setting up a Sciline pipeline data. @@ -55,7 +49,9 @@ def default_parameters() -> dict: """Return default parameters for the NeXus Freia workflow.""" return { NeXusDetectorName: "multiblade_detector", - SampleRotationOffset[RunType]: sc.scalar(0.0, unit='deg'), + DetectorBankSizes: { + "multiblade_detector": {"strip": 64, "blade": 32, "wire": 32}, + }, DetectorSpatialResolution: 0.0025 * sc.units.m, LookupTableRelativeErrorThreshold: { "multiblade_detector": float('inf'), @@ -111,6 +107,9 @@ def FreiaWorkflow( measurement without a sample, taken with matching slit and chopper settings. Set its ``SampleSurfaceNormal`` to the sample run's orientation. + To skip footprint correction, set + ``workflow[Sample] = workflow[ReducibleData[SampleRun]]``. + Monitor normalization requires an incident monitor selected through ``NeXusName[IncidentMonitor]``, or supplied as ``WavelengthMonitor[RunType]``. @@ -123,7 +122,12 @@ def FreiaWorkflow( 'analytical', 'simulation', and 'file'. See https://scipp.github.io/ess/reduce/user-guide/unwrap/lut-building-methods.html """ - workflow = beamline.LoadNeXusWorkflow(wavelength_from=wavelength_from, **kwargs) + workflow = GenericUnwrapWorkflow( + run_types=[SampleRun, ReferenceRun], + monitor_types=[IncidentMonitor], + wavelength_from=wavelength_from, + **kwargs, + ) for provider in providers: workflow.insert(provider) insert_run_normalization(workflow, run_norm) diff --git a/packages/essreflectometry/tests/freia/workflow_test.py b/packages/essreflectometry/tests/freia/workflow_test.py index 71543bd06..3c2ac20cf 100644 --- a/packages/essreflectometry/tests/freia/workflow_test.py +++ b/packages/essreflectometry/tests/freia/workflow_test.py @@ -10,7 +10,6 @@ from ess.freia.corrections import RunNormalization from ess.freia.types import ( DetectorRegionOfInterest, - SampleIlluminatedFraction, SampleSurfaceNormal, WavelengthMonitor, ) @@ -21,6 +20,7 @@ ReducibleData, ReferenceRun, ReflectivityOverQ, + Sample, SampleRun, SampleSize, WavelengthBins, @@ -105,7 +105,7 @@ def test_reduce_reflectivity_with_monitor_and_footprint(): def test_rebinning_integrates_before_dividing(): wf = _make_workflow(RunNormalization.none) - wf[SampleIlluminatedFraction] = sc.scalar(1.0) + wf[Sample] = wf[ReducibleData[SampleRun]] wf[QBins] = sc.array(dims=['Q'], values=[1.0, 4.0], unit='1/angstrom') result = wf.compute(ReflectivityOverQ) From 8ef8a81a0f6c1b070c80ff934594facd8ae71313 Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Thu, 17 Sep 2026 15:30:39 +0200 Subject: [PATCH 10/16] ci: upload files to pooch --- .../freia/freia-mcstas-visualization.ipynb | 8 +- .../user-guide/freia/freia-reflectivity.ipynb | 32 ++--- .../essreflectometry/src/ess/freia/data.py | 30 +++- .../tests/freia/mcstas_test.py | 131 ++++++------------ 4 files changed, 82 insertions(+), 119 deletions(-) diff --git a/packages/essreflectometry/docs/user-guide/freia/freia-mcstas-visualization.ipynb b/packages/essreflectometry/docs/user-guide/freia/freia-mcstas-visualization.ipynb index 0d958afce..7130a0535 100644 --- a/packages/essreflectometry/docs/user-guide/freia/freia-mcstas-visualization.ipynb +++ b/packages/essreflectometry/docs/user-guide/freia/freia-mcstas-visualization.ipynb @@ -22,6 +22,7 @@ "import scipp as sc\n", "\n", "from ess import freia\n", + "from ess.freia import data\n", "from ess.reduce.nexus.types import DiskChoppers, Filename, RawDetector\n", "from ess.reduce.unwrap import (\n", " ChopperFrameSequence,\n", @@ -41,7 +42,7 @@ "source": [ "## Select a run\n", "\n", - "Set the path to a local event file, relative to this notebook, and the resolution used to reconstruct wavelengths." + "Download the example sample run and set the resolution used to reconstruct wavelengths. The file is cached locally by `pooch` (`pip install pooch`). To use your own data, replace the download function with a file path." ] }, { @@ -52,7 +53,7 @@ "outputs": [], "source": [ "freia_mcstas = freia.FreiaMcStasWorkflow(wavelength_from='analytical')\n", - "freia_mcstas[Filename[SampleRun]] = '../../../265305.h5'\n", + "freia_mcstas[Filename[SampleRun]] = data.freia_mcstas_reference_run()\n", "freia_mcstas[KeepEventTimeOffset] = True\n", "freia_mcstas[TimeResolution] = sc.scalar(20.0, unit='us')\n", "freia_mcstas[DistanceResolution] = sc.scalar(0.1, unit='m')" @@ -211,9 +212,6 @@ "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.12.13" - }, - "nbsphinx": { - "execute": "never" } }, "nbformat": 4, diff --git a/packages/essreflectometry/docs/user-guide/freia/freia-reflectivity.ipynb b/packages/essreflectometry/docs/user-guide/freia/freia-reflectivity.ipynb index e8bae3676..d9a8bc092 100644 --- a/packages/essreflectometry/docs/user-guide/freia/freia-reflectivity.ipynb +++ b/packages/essreflectometry/docs/user-guide/freia/freia-reflectivity.ipynb @@ -9,7 +9,7 @@ "\n", "Compute one reflectivity curve for each incident beam using a sample measurement and a direct-beam measurement taken without a sample. Normalize by an incident wavelength monitor, select matching peaks, and compare the three curves with the silicon reflectivity reference.\n", "\n", - "The sample and direct-beam runs must have matching slit and chopper settings. This example uses local files with WFM chopper settings." + "The sample and direct-beam runs must have matching slit and chopper settings. This example downloads and caches WFM sample and direct-beam runs using `pooch` (`pip install pooch`)." ] }, { @@ -23,12 +23,11 @@ "import numpy as np\n", "import scipp as sc\n", "\n", - "from ess.freia import FreiaMcStasWorkflow\n", + "from ess.freia import FreiaMcStasWorkflow, data\n", "from ess.freia.corrections import RunNormalization\n", "from ess.freia.types import (\n", " DetectorRegionOfInterest,\n", " IncidentMonitor,\n", - " SampleIlluminatedFraction,\n", " WavelengthMonitor,\n", ")\n", "from ess.reduce.nexus.types import NeXusName\n", @@ -37,8 +36,10 @@ " CorrectedDetector,\n", " Filename,\n", " QBins,\n", + " ReducibleData,\n", " ReferenceRun,\n", " ReflectivityOverQ,\n", + " Sample,\n", " SampleRun,\n", " SampleSize,\n", " WavelengthBins,\n", @@ -53,17 +54,7 @@ "source": [ "## Select the runs\n", "\n", - "Paths are relative to this notebook. `SampleLambda` measures the incident spectrum before the sample in these files. If no incident monitor is available, use `RunNormalization.none`; the two runs must then already share an exposure and flux scale.\n", - "\n", - "Gravity correction is enabled for these runs. For data simulated without gravity, set it to zero after creating the workflow:\n", - "\n", - "```python\n", - "from ess.reduce.nexus.types import GravityVector\n", - "\n", - "workflow[GravityVector] = sc.vector([0.0, 0.0, 0.0], unit='m/s^2')\n", - "```\n", - "\n", - "The Q grid covers all three beams. Bins without direct-beam coverage are masked separately for each curve." + "`SampleLambda` is a monitor in the McStas model that measures the wavelength spectrum before the sample. If no incident monitor is available, use `RunNormalization.none`; the two runs must then already share an exposure and flux scale.\n" ] }, { @@ -74,8 +65,8 @@ "outputs": [], "source": [ "workflow = FreiaMcStasWorkflow(run_norm=RunNormalization.monitor_histogram)\n", - "workflow[Filename[SampleRun]] = '../../../265305.h5'\n", - "workflow[Filename[ReferenceRun]] = '../../../265301.h5'\n", + "workflow[Filename[SampleRun]] = data.freia_mcstas_sample_run()\n", + "workflow[Filename[ReferenceRun]] = data.freia_mcstas_reference_run()\n", "workflow[NeXusName[IncidentMonitor]] = 'SampleLambda'\n", "\n", "workflow[WavelengthBins] = sc.linspace('wavelength', 2.0, 10.0, 81, unit='angstrom')\n", @@ -170,6 +161,8 @@ " 'Direct-beam monitor': inputs[WavelengthMonitor[ReferenceRun]],\n", " },\n", " title='Incident wavelength spectra',\n", + " norm='log',\n", + " vmin=1e0,\n", ")" ] }, @@ -230,7 +223,7 @@ " }\n", " beam_size = beam_sizes[label]\n", " if beam_size is None:\n", - " reduction[SampleIlluminatedFraction] = sc.scalar(1.0)\n", + " reduction[Sample] = reduction[ReducibleData[SampleRun]]\n", " else:\n", " reduction[SampleSize[SampleRun]] = sample_size\n", " reduction[BeamSize[SampleRun]] = beam_size\n", @@ -257,7 +250,7 @@ "metadata": {}, "outputs": [], "source": [ - "reference_file = '../../../../../Si-15SiO2-air.txt'\n", + "reference_file = data.freia_mcstas_silicon_reflectivity()\n", "q, r = np.loadtxt(reference_file, unpack=True)\n", "si_reference = sc.DataArray(\n", " sc.array(dims=['Q'], values=r, unit='dimensionless'),\n", @@ -312,9 +305,6 @@ "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.12.13" - }, - "nbsphinx": { - "execute": "never" } }, "nbformat": 4, diff --git a/packages/essreflectometry/src/ess/freia/data.py b/packages/essreflectometry/src/ess/freia/data.py index 277709ad4..44b91e458 100644 --- a/packages/essreflectometry/src/ess/freia/data.py +++ b/packages/essreflectometry/src/ess/freia/data.py @@ -1,5 +1,9 @@ # SPDX-License-Identifier: BSD-3-Clause # Copyright (c) 2025 Scipp contributors (https://github.com/scipp) +"""Data for tests and documentation with FREIA.""" + +from pathlib import Path + from ess.reduce.data import make_registry from ..reflectometry.types import Filename, ReferenceRun, SampleRun @@ -7,18 +11,32 @@ _registry = make_registry( "ess/freia", version="1", - files={}, + files={ + # WFM runs 265305 (silicon with natural oxide) and 265301 (no sample). + "mcstas-wfm-silicon.h5": "md5:07892a4b705faa769f5be8d3601b529e", + "mcstas-wfm-direct-beam.h5": "md5:b7a2f51615bdaec6957bf3af30751ebf", + "Si-15SiO2-air.txt": "md5:2bac8babb5f7042fd2e3cbd62e22a278", + }, ) def freia_mcstas_sample_run() -> Filename[SampleRun]: - """Return path to the McStas sample events file.""" - return Filename[SampleRun](_registry.get_path("mcstas-sample.h5")) + """Return the WFM McStas run for silicon with natural oxide.""" + return Filename[SampleRun](_registry.get_path("mcstas-wfm-silicon.h5")) def freia_mcstas_reference_run() -> Filename[ReferenceRun]: - """Return path to the McStas reference events file.""" - return Filename[ReferenceRun](_registry.get_path("mcstas-reference.h5")) + """Return the matching WFM McStas direct-beam run without a sample.""" + return Filename[ReferenceRun](_registry.get_path("mcstas-wfm-direct-beam.h5")) + + +def freia_mcstas_silicon_reflectivity() -> Path: + """Return the silicon reflectivity table used in the McStas sample run.""" + return _registry.get_path("Si-15SiO2-air.txt") -__all__ = ["freia_mcstas_reference_run", "freia_mcstas_sample_run"] +__all__ = [ + "freia_mcstas_reference_run", + "freia_mcstas_sample_run", + "freia_mcstas_silicon_reflectivity", +] diff --git a/packages/essreflectometry/tests/freia/mcstas_test.py b/packages/essreflectometry/tests/freia/mcstas_test.py index c651aba5b..bb90bd485 100644 --- a/packages/essreflectometry/tests/freia/mcstas_test.py +++ b/packages/essreflectometry/tests/freia/mcstas_test.py @@ -10,6 +10,7 @@ from scipp.testing import assert_allclose, assert_identical from ess.freia import FreiaMcStasWorkflow +from ess.freia.data import freia_mcstas_reference_run, freia_mcstas_sample_run from ess.freia.mcstas import mcstas_detector_geometry from ess.freia.types import DetectorRegionOfInterest, IncidentMonitor, WavelengthMonitor from ess.reflectometry.types import CorrectedDetector, SampleRun, WavelengthBins @@ -22,49 +23,9 @@ def _component(components, name, position): return group -@pytest.fixture -def mcstas_file(tmp_path: Path) -> Path: - """Minimal file for the McStasToX reader; FREIA has no published event fixture.""" - filename = tmp_path / 'freia.h5' - with h5py.File(filename, 'w') as f: - entry = f.create_group('entry1') - entry.create_group('data') - simulation = entry.create_group('simulation') - simulation.attrs['program'] = np.bytes_('3.7.18, git') - simulation.create_group('Param') - components = entry.create_group('instrument/components') - _component(components, 'Source', [0.0, 0.0, 0.0]) - _component(components, 'Arm_Sample', [0.0, 0.0, 20.0]) - detector = _component(components, 'Multiblade', [0.0, 0.0, 20.0]) - detector['Rotation'][...] = [[0, 1, 0], [-1, 0, 0], [0, 0, 1]] - geometry = detector.create_group('Geometry') - geometry.attrs['Shape identifier'] = np.bytes_('4') - geometry.attrs['radius'] = np.bytes_('3') - bins = detector.create_group('output/BINS') - for key, value in { - 'xvar': 'th', - 'yvar': 'y', - 'xlabel': 'theta', - 'ylabel': 'height', - }.items(): - bins.attrs[key] = np.bytes_(value) - bins['theta'] = [1.0, 2.0] - bins['height'] = [-0.001, 0.001] - bins['pixels'] = [[12, 10], [99, 101]] - output = detector.create_group('output/detector_events') - output.attrs['variables'] = np.bytes_('p t id') - output['events'] = [ - [2.0, 0.025, 10.0], - [3.0, 0.025 + 1 / 14, 10.0], - [0.0, 0.025, 99.0], - [4.0, 0.025, 12.0], - ] - return filename - - -def test_load_detector_and_compute_q(mcstas_file): +def test_load_detector_and_compute_q(): workflow = FreiaMcStasWorkflow() - workflow[Filename[SampleRun]] = mcstas_file + workflow[Filename[SampleRun]] = freia_mcstas_sample_run() workflow[DetectorRegionOfInterest[SampleRun]] = {} workflow[WavelengthBins] = sc.array( dims=['wavelength'], values=[1.0, 12.0], unit='angstrom' @@ -73,40 +34,42 @@ def test_load_detector_and_compute_q(mcstas_file): result = workflow.compute((RawDetector[SampleRun], CorrectedDetector[SampleRun])) raw = result[RawDetector[SampleRun]] - image = sc.sort(raw.bins.sum(), 'pixel_id') - assert_identical( - image.coords['pixel_id'], - sc.array(dims=['pixel_id'], values=[10, 12, 99, 101], unit=None), - ) - assert_identical( - image.data, - sc.array( - dims=['pixel_id'], - values=[5.0, 4.0, 0.0, 0.0], - variances=[13.0, 16.0, 0.0, 0.0], - unit='counts', - ), - ) + assert raw.sizes == {'pixel_id': 2048 * 64} events = raw.bins.constituents['data'] + assert events.sizes == {'event': 555874} assert_allclose( - events.coords['event_time_offset'], - sc.full(sizes=events.sizes, value=0.025, unit='s'), + events.sum().data, + sc.scalar(744336.6691526351, variance=679997917.939667, unit='counts'), ) + assert events.coords['wavelength_from_mcstas'].unit == 'angstrom' + offsets = events.coords['event_time_offset'] + assert offsets.min() >= sc.scalar(0.0, unit='s') + assert offsets.max() < sc.scalar(1 / 14, unit='s') q = result[CorrectedDetector[SampleRun]].bins.constituents['data'].coords['Q'] - assert sc.isfinite(q).all().value + assert sc.isfinite(q).any().value -def test_load_histogram_geometry(mcstas_file): - with h5py.File(mcstas_file, 'r+') as f: - detector = f['entry1/instrument/components/0002_Multiblade'] - del detector['output'] +def test_load_histogram_geometry(tmp_path: Path): + filename = tmp_path / 'histogram.h5' + with h5py.File(filename, 'w') as f: + entry = f.create_group('entry1') + entry.create_group('data') + simulation = entry.create_group('simulation') + simulation.attrs['program'] = np.bytes_('3.7.18, git') + simulation.create_group('Param') + components = entry.create_group('instrument/components') + detector = _component(components, 'Multiblade', [0.0, 0.0, 20.0]) + detector['Rotation'][...] = [[0, 1, 0], [-1, 0, 0], [0, 0, 1]] + geometry = detector.create_group('Geometry') + geometry.attrs['Shape identifier'] = np.bytes_('4') + geometry.attrs['radius'] = np.bytes_('3') histogram = detector.create_group('output/histogram') histogram.attrs['xlabel'] = np.bytes_('theta [deg]') histogram.attrs['ylabel'] = np.bytes_('Height [cm]') histogram['theta__deg_'] = [0.0, 90.0] histogram['Height__cm_'] = [-25, 25] - detector = mcstas_detector_geometry(mcstas_file, 'Multiblade') + detector = mcstas_detector_geometry(filename, 'Multiblade') assert_allclose( detector.coords['position'], @@ -118,35 +81,29 @@ def test_load_histogram_geometry(mcstas_file): ) -def test_load_monitor(mcstas_file): - with h5py.File(mcstas_file, 'r+') as f: - monitor = _component( - f['entry1/instrument/components'], 'IncidentLambda', [0, 0, 19] - ) - histogram = monitor.create_group('output/spectrum') - histogram.attrs['xvar'] = np.bytes_('L') - histogram.attrs['xlimits'] = np.bytes_('1 5') - histogram['data'] = [10.0, 20.0] - histogram['errors'] = [3.0, 4.0] +@pytest.mark.parametrize( + 'filename', [freia_mcstas_sample_run, freia_mcstas_reference_run] +) +def test_load_monitor(filename): + path = filename() workflow = FreiaMcStasWorkflow() - workflow[Filename[SampleRun]] = mcstas_file - workflow[NeXusName[IncidentMonitor]] = 'IncidentLambda' + workflow[Filename[SampleRun]] = path + workflow[NeXusName[IncidentMonitor]] = 'SampleLambda' result = workflow.compute(WavelengthMonitor[SampleRun]) assert_identical( - result, - sc.DataArray( + result.coords['wavelength'], + sc.linspace('wavelength', 0.0, 25.0, 101, unit='angstrom'), + ) + with h5py.File(path, 'r') as f: + monitor = f['entry1/data/SampleLambda_dat'] + assert_identical( + result.data, sc.array( dims=['wavelength'], - values=[10.0, 20.0], - variances=[9.0, 16.0], + values=monitor['data'][:], + variances=monitor['errors'][:] ** 2, unit='counts', ), - coords={ - 'wavelength': sc.array( - dims=['wavelength'], values=[1.0, 3.0, 5.0], unit='angstrom' - ) - }, - ), - ) + ) From 0fa59e71c07a48074fbaf02c76324ff9abc26768 Mon Sep 17 00:00:00 2001 From: jokasimr Date: Thu, 17 Sep 2026 15:41:50 +0200 Subject: [PATCH 11/16] Update packages/essreflectometry/docs/user-guide/freia/freia-wavelength-lookup-table.ipynb Co-authored-by: Neil Vaytet <39047984+nvaytet@users.noreply.github.com> --- .../docs/user-guide/freia/freia-wavelength-lookup-table.ipynb | 1 - 1 file changed, 1 deletion(-) diff --git a/packages/essreflectometry/docs/user-guide/freia/freia-wavelength-lookup-table.ipynb b/packages/essreflectometry/docs/user-guide/freia/freia-wavelength-lookup-table.ipynb index 8318f01a1..35db2d668 100644 --- a/packages/essreflectometry/docs/user-guide/freia/freia-wavelength-lookup-table.ipynb +++ b/packages/essreflectometry/docs/user-guide/freia/freia-wavelength-lookup-table.ipynb @@ -17,7 +17,6 @@ "metadata": {}, "outputs": [], "source": [ - "%matplotlib widget\n", "from pathlib import Path\n", "\n", "import plopp as pp\n", From 4f90a0dc2ccd8725bd2d7ff2de79af174c4c30f6 Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Thu, 17 Sep 2026 15:48:31 +0200 Subject: [PATCH 12/16] fix --- .../docs/user-guide/freia/freia-reflectivity.ipynb | 1 - 1 file changed, 1 deletion(-) diff --git a/packages/essreflectometry/docs/user-guide/freia/freia-reflectivity.ipynb b/packages/essreflectometry/docs/user-guide/freia/freia-reflectivity.ipynb index d9a8bc092..5477748a3 100644 --- a/packages/essreflectometry/docs/user-guide/freia/freia-reflectivity.ipynb +++ b/packages/essreflectometry/docs/user-guide/freia/freia-reflectivity.ipynb @@ -19,7 +19,6 @@ "metadata": {}, "outputs": [], "source": [ - "%matplotlib inline\n", "import numpy as np\n", "import scipp as sc\n", "\n", From c3e329b77e3d3fedc0365730695275b0cf279803 Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Thu, 17 Sep 2026 15:53:03 +0200 Subject: [PATCH 13/16] pixi lock --- pixi.lock | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/pixi.lock b/pixi.lock index 4892b941a..c8a52ee22 100644 --- a/pixi.lock +++ b/pixi.lock @@ -4763,6 +4763,7 @@ environments: - pypi: https://files.pythonhosted.org/packages/14/15/5574111ae50dd6e879456888c0eadd4c5a869959775854e18e18a6b345f3/propcache-0.5.2-cp312-cp312-manylinux2014_x86_64.manylinux_2_17_x86_64.manylinux_2_28_x86_64.whl - pypi: https://files.pythonhosted.org/packages/16/8d/93be7e0f7fa915a576859b3bfac7a7baa3303181c44d7db7eefbd3e8a69f/sphinxcontrib_mermaid-2.0.2-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/18/67/36e9267722cc04a6b9f15c7f3441c2363321a3ea07da7ae0c0707beb2a9c/typing_extensions-4.15.0-py3-none-any.whl + - pypi: https://files.pythonhosted.org/packages/1a/84/8c38176f97ee333c5b0d3f003548718ee87abdac67d4a2d5c00183c872bf/mcstastox-0.0.12-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/1d/33/f1c6a276de27b7d7339a34749cc33fa87f077f921969c47185d34a887ae2/gast-0.7.0-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/1e/5e/d4e9f1a599fb8e573b7b87160658329fbf28d19eac2718f51fc3def3aa5a/idna-3.18-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/1e/77/dc8c558f7593132cf8fefec57c4f60c83b16941c574ac5f619abb3ae7933/dill-0.4.1-py3-none-any.whl @@ -5000,6 +5001,7 @@ environments: - pypi: https://files.pythonhosted.org/packages/12/b3/88c0ef22878c86035f058df0ac6c171319ffd0aa52a406455ed3a3847566/ipympl-0.10.0-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/16/8d/93be7e0f7fa915a576859b3bfac7a7baa3303181c44d7db7eefbd3e8a69f/sphinxcontrib_mermaid-2.0.2-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/18/67/36e9267722cc04a6b9f15c7f3441c2363321a3ea07da7ae0c0707beb2a9c/typing_extensions-4.15.0-py3-none-any.whl + - pypi: https://files.pythonhosted.org/packages/1a/84/8c38176f97ee333c5b0d3f003548718ee87abdac67d4a2d5c00183c872bf/mcstastox-0.0.12-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/1d/33/f1c6a276de27b7d7339a34749cc33fa87f077f921969c47185d34a887ae2/gast-0.7.0-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/1e/5e/d4e9f1a599fb8e573b7b87160658329fbf28d19eac2718f51fc3def3aa5a/idna-3.18-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/1e/77/dc8c558f7593132cf8fefec57c4f60c83b16941c574ac5f619abb3ae7933/dill-0.4.1-py3-none-any.whl @@ -5238,6 +5240,7 @@ environments: - pypi: https://files.pythonhosted.org/packages/16/8d/93be7e0f7fa915a576859b3bfac7a7baa3303181c44d7db7eefbd3e8a69f/sphinxcontrib_mermaid-2.0.2-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/18/67/36e9267722cc04a6b9f15c7f3441c2363321a3ea07da7ae0c0707beb2a9c/typing_extensions-4.15.0-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/19/95/6195171e385007300f0f5574592e467c568becce2d937a0b6804f218bc49/pydantic_core-2.46.4-cp312-cp312-macosx_11_0_arm64.whl + - pypi: https://files.pythonhosted.org/packages/1a/84/8c38176f97ee333c5b0d3f003548718ee87abdac67d4a2d5c00183c872bf/mcstastox-0.0.12-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/1d/33/f1c6a276de27b7d7339a34749cc33fa87f077f921969c47185d34a887ae2/gast-0.7.0-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/1e/5e/d4e9f1a599fb8e573b7b87160658329fbf28d19eac2718f51fc3def3aa5a/idna-3.18-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/1e/77/dc8c558f7593132cf8fefec57c4f60c83b16941c574ac5f619abb3ae7933/dill-0.4.1-py3-none-any.whl @@ -5481,6 +5484,7 @@ environments: - pypi: https://files.pythonhosted.org/packages/18/67/36e9267722cc04a6b9f15c7f3441c2363321a3ea07da7ae0c0707beb2a9c/typing_extensions-4.15.0-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/19/e8/6026ed58a64563186a9ee3f29f41261fd1828f527dd93d33b60feca63352/contourpy-1.3.3-cp312-cp312-win_amd64.whl - pypi: https://files.pythonhosted.org/packages/1a/1f/86b4d15221096cb5500bcd73bf350745749e3ba056cdd7a7f75f126f154e/scipp-26.3.1-cp312-cp312-win_amd64.whl + - pypi: https://files.pythonhosted.org/packages/1a/84/8c38176f97ee333c5b0d3f003548718ee87abdac67d4a2d5c00183c872bf/mcstastox-0.0.12-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/1d/33/f1c6a276de27b7d7339a34749cc33fa87f077f921969c47185d34a887ae2/gast-0.7.0-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/1e/5e/d4e9f1a599fb8e573b7b87160658329fbf28d19eac2718f51fc3def3aa5a/idna-3.18-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/1e/77/dc8c558f7593132cf8fefec57c4f60c83b16941c574ac5f619abb3ae7933/dill-0.4.1-py3-none-any.whl @@ -10047,6 +10051,7 @@ environments: - pypi: https://files.pythonhosted.org/packages/12/b3/88c0ef22878c86035f058df0ac6c171319ffd0aa52a406455ed3a3847566/ipympl-0.10.0-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/14/15/5574111ae50dd6e879456888c0eadd4c5a869959775854e18e18a6b345f3/propcache-0.5.2-cp312-cp312-manylinux2014_x86_64.manylinux_2_17_x86_64.manylinux_2_28_x86_64.whl - pypi: https://files.pythonhosted.org/packages/18/67/36e9267722cc04a6b9f15c7f3441c2363321a3ea07da7ae0c0707beb2a9c/typing_extensions-4.15.0-py3-none-any.whl + - pypi: https://files.pythonhosted.org/packages/1a/84/8c38176f97ee333c5b0d3f003548718ee87abdac67d4a2d5c00183c872bf/mcstastox-0.0.12-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/1d/33/f1c6a276de27b7d7339a34749cc33fa87f077f921969c47185d34a887ae2/gast-0.7.0-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/1e/5e/d4e9f1a599fb8e573b7b87160658329fbf28d19eac2718f51fc3def3aa5a/idna-3.18-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/1e/77/dc8c558f7593132cf8fefec57c4f60c83b16941c574ac5f619abb3ae7933/dill-0.4.1-py3-none-any.whl @@ -10225,6 +10230,7 @@ environments: - pypi: https://files.pythonhosted.org/packages/12/aa/fb2a0649fdeef5ab7072d221e8f4df164098792c813af6c87e2581cfa860/mpltoolbox-26.2.0-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/12/b3/88c0ef22878c86035f058df0ac6c171319ffd0aa52a406455ed3a3847566/ipympl-0.10.0-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/18/67/36e9267722cc04a6b9f15c7f3441c2363321a3ea07da7ae0c0707beb2a9c/typing_extensions-4.15.0-py3-none-any.whl + - pypi: https://files.pythonhosted.org/packages/1a/84/8c38176f97ee333c5b0d3f003548718ee87abdac67d4a2d5c00183c872bf/mcstastox-0.0.12-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/1d/33/f1c6a276de27b7d7339a34749cc33fa87f077f921969c47185d34a887ae2/gast-0.7.0-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/1e/5e/d4e9f1a599fb8e573b7b87160658329fbf28d19eac2718f51fc3def3aa5a/idna-3.18-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/1e/77/dc8c558f7593132cf8fefec57c4f60c83b16941c574ac5f619abb3ae7933/dill-0.4.1-py3-none-any.whl @@ -10402,6 +10408,7 @@ environments: - pypi: https://files.pythonhosted.org/packages/12/b3/88c0ef22878c86035f058df0ac6c171319ffd0aa52a406455ed3a3847566/ipympl-0.10.0-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/18/67/36e9267722cc04a6b9f15c7f3441c2363321a3ea07da7ae0c0707beb2a9c/typing_extensions-4.15.0-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/19/95/6195171e385007300f0f5574592e467c568becce2d937a0b6804f218bc49/pydantic_core-2.46.4-cp312-cp312-macosx_11_0_arm64.whl + - pypi: https://files.pythonhosted.org/packages/1a/84/8c38176f97ee333c5b0d3f003548718ee87abdac67d4a2d5c00183c872bf/mcstastox-0.0.12-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/1d/33/f1c6a276de27b7d7339a34749cc33fa87f077f921969c47185d34a887ae2/gast-0.7.0-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/1e/5e/d4e9f1a599fb8e573b7b87160658329fbf28d19eac2718f51fc3def3aa5a/idna-3.18-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/1e/77/dc8c558f7593132cf8fefec57c4f60c83b16941c574ac5f619abb3ae7933/dill-0.4.1-py3-none-any.whl @@ -10586,6 +10593,7 @@ environments: - pypi: https://files.pythonhosted.org/packages/18/67/36e9267722cc04a6b9f15c7f3441c2363321a3ea07da7ae0c0707beb2a9c/typing_extensions-4.15.0-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/19/e8/6026ed58a64563186a9ee3f29f41261fd1828f527dd93d33b60feca63352/contourpy-1.3.3-cp312-cp312-win_amd64.whl - pypi: https://files.pythonhosted.org/packages/1a/1f/86b4d15221096cb5500bcd73bf350745749e3ba056cdd7a7f75f126f154e/scipp-26.3.1-cp312-cp312-win_amd64.whl + - pypi: https://files.pythonhosted.org/packages/1a/84/8c38176f97ee333c5b0d3f003548718ee87abdac67d4a2d5c00183c872bf/mcstastox-0.0.12-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/1d/33/f1c6a276de27b7d7339a34749cc33fa87f077f921969c47185d34a887ae2/gast-0.7.0-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/1e/5e/d4e9f1a599fb8e573b7b87160658329fbf28d19eac2718f51fc3def3aa5a/idna-3.18-py3-none-any.whl - pypi: https://files.pythonhosted.org/packages/1e/77/dc8c558f7593132cf8fefec57c4f60c83b16941c574ac5f619abb3ae7933/dill-0.4.1-py3-none-any.whl @@ -17534,6 +17542,7 @@ packages: requires_dist: - dask>=2022.1.0 - graphviz>=0.20 + - mcstastox>=0.0.11 - python-dateutil>=2.9.0 - plopp>=26.5.0 - orsopy>=1.2 From 5b6ccc52cbfe828c6856d3f46b129f748e264591 Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Thu, 17 Sep 2026 15:56:13 +0200 Subject: [PATCH 14/16] docs --- .../user-guide/freia/freia-reflectivity.ipynb | 32 ++++++++++--------- .../freia/freia-wavelength-lookup-table.ipynb | 10 ++---- 2 files changed, 20 insertions(+), 22 deletions(-) diff --git a/packages/essreflectometry/docs/user-guide/freia/freia-reflectivity.ipynb b/packages/essreflectometry/docs/user-guide/freia/freia-reflectivity.ipynb index 5477748a3..85b964e90 100644 --- a/packages/essreflectometry/docs/user-guide/freia/freia-reflectivity.ipynb +++ b/packages/essreflectometry/docs/user-guide/freia/freia-reflectivity.ipynb @@ -20,6 +20,7 @@ "outputs": [], "source": [ "import numpy as np\n", + "import plopp as pp\n", "import scipp as sc\n", "\n", "from ess.freia import FreiaMcStasWorkflow, data\n", @@ -255,23 +256,24 @@ " sc.array(dims=['Q'], values=r, unit='dimensionless'),\n", " coords={'Q': sc.array(dims=['Q'], values=q, unit='1/angstrom')},\n", ")\n", - "# Display each histogram value at its Q-bin center.\n", - "plot_curves = {\n", - " label: curve.assign_coords(Q=sc.midpoints(curve.coords['Q']))\n", - " for label, curve in reflectivities.items()\n", - "}\n", - "qbins = workflow.compute(QBins)\n", - "comparison = sc.plot(\n", - " {'Si reference': si_reference, **plot_curves},\n", - " logx=True,\n", + "# Display finite, unmasked values at their Q-bin centers.\n", + "plot_curves = {}\n", + "for label, curve in reflectivities.items():\n", + " points = curve.assign_coords(Q=sc.midpoints(curve.coords['Q']))\n", + " valid = sc.isfinite(points.data) & sc.isfinite(sc.variances(points.data))\n", + " valid &= ~points.masks['direct_beam']\n", + " plot_curves[label] = points[valid]\n", + "q_range = sc.concat([curve.coords['Q'] for curve in plot_curves.values()], dim='Q')\n", + "comparison = pp.plot(\n", + " {**plot_curves, 'Si reference': si_reference},\n", " logy=True,\n", - " linestyle={'Si reference': '-'},\n", - " marker={'Si reference': 'none'},\n", + " ymin=1e-8,\n", + " ymax=1.0,\n", + " xmin=q_range.min(),\n", + " xmax=q_range.max(),\n", + " ls=dict.fromkeys(plot_curves, 'solid') | {'Si reference': 'dashed'},\n", + " marker='none',\n", " color={'Si reference': 'black'},\n", - " markersize=4,\n", - " figsize=(8, 5),\n", - " xmin=qbins[0],\n", - " xmax=qbins[-1],\n", " title='FREIA reflectivity by beam',\n", " ylabel='Reflectivity',\n", ")\n", diff --git a/packages/essreflectometry/docs/user-guide/freia/freia-wavelength-lookup-table.ipynb b/packages/essreflectometry/docs/user-guide/freia/freia-wavelength-lookup-table.ipynb index 35db2d668..1fd5454e7 100644 --- a/packages/essreflectometry/docs/user-guide/freia/freia-wavelength-lookup-table.ipynb +++ b/packages/essreflectometry/docs/user-guide/freia/freia-wavelength-lookup-table.ipynb @@ -17,13 +17,12 @@ "metadata": {}, "outputs": [], "source": [ - "from pathlib import Path\n", - "\n", "import plopp as pp\n", "import scipp as sc\n", "from scippnexus import NXdetector\n", "\n", "from ess import freia\n", + "from ess.freia import data\n", "from ess.reduce.nexus.types import DiskChoppers, Filename\n", "from ess.reduce.unwrap import (\n", " ChopperFrameSequence,\n", @@ -42,7 +41,7 @@ "source": [ "## Use the WFM configuration\n", "\n", - "Set the path to a local file, relative to this notebook, and the lookup-table resolution. This example uses 20 µs time resolution and 10 cm flight-path resolution.\n", + "Download the example WFM sample run and set the lookup-table resolution. The file is cached locally by `pooch` (`pip install pooch`). This example uses 20 µs time resolution and 10 cm flight-path resolution.\n", "\n", "For a non-WFM run, insert `freia.mcstas.non_wfm_choppers` into the workflow after creating it: `freia_mcstas.insert(freia.mcstas.non_wfm_choppers)`." ] @@ -54,7 +53,7 @@ "metadata": {}, "outputs": [], "source": [ - "filename = Path('../../../265224.h5')\n", + "filename = data.freia_mcstas_sample_run()\n", "\n", "freia_mcstas = freia.FreiaMcStasWorkflow(wavelength_from='analytical')\n", "freia_mcstas[Filename[SampleRun]] = filename\n", @@ -158,9 +157,6 @@ }, "language_info": { "name": "python" - }, - "nbsphinx": { - "execute": "never" } }, "nbformat": 4, From 6126c9c6406a6e46c5eef1e33a981f9138c5fdce Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Thu, 17 Sep 2026 16:12:45 +0200 Subject: [PATCH 15/16] docs --- .../user-guide/freia/freia-reflectivity.ipynb | 37 +++++++++++++++---- 1 file changed, 30 insertions(+), 7 deletions(-) diff --git a/packages/essreflectometry/docs/user-guide/freia/freia-reflectivity.ipynb b/packages/essreflectometry/docs/user-guide/freia/freia-reflectivity.ipynb index 85b964e90..2e7af070b 100644 --- a/packages/essreflectometry/docs/user-guide/freia/freia-reflectivity.ipynb +++ b/packages/essreflectometry/docs/user-guide/freia/freia-reflectivity.ipynb @@ -256,14 +256,15 @@ " sc.array(dims=['Q'], values=r, unit='dimensionless'),\n", " coords={'Q': sc.array(dims=['Q'], values=q, unit='1/angstrom')},\n", ")\n", - "# Display finite, unmasked values at their Q-bin centers.\n", + "# Keep the Q bin edges so the reduced curves are drawn as histograms.\n", "plot_curves = {}\n", "for label, curve in reflectivities.items():\n", - " points = curve.assign_coords(Q=sc.midpoints(curve.coords['Q']))\n", - " valid = sc.isfinite(points.data) & sc.isfinite(sc.variances(points.data))\n", - " valid &= ~points.masks['direct_beam']\n", - " plot_curves[label] = points[valid]\n", - "q_range = sc.concat([curve.coords['Q'] for curve in plot_curves.values()], dim='Q')\n", + " valid = sc.isfinite(curve.data) & sc.isfinite(sc.variances(curve.data))\n", + " valid &= ~curve.masks['direct_beam']\n", + " plot_curves[label] = curve.assign(\n", + " sc.where(valid, curve.data, sc.scalar(float('nan'), unit=curve.unit))\n", + " )\n", + "q_range = workflow.compute(QBins)\n", "comparison = pp.plot(\n", " {**plot_curves, 'Si reference': si_reference},\n", " logy=True,\n", @@ -281,9 +282,31 @@ ] }, { - "cell_type": "markdown", + "cell_type": "code", + "execution_count": null, "id": "16", "metadata": {}, + "outputs": [], + "source": [ + "critical_edge = pp.plot(\n", + " {'0.3° beam': plot_curves['0.3° beam'], 'Si reference': si_reference},\n", + " xmin=sc.scalar(0.006, unit='1/angstrom'),\n", + " xmax=sc.scalar(0.016, unit='1/angstrom'),\n", + " ymin=0.0,\n", + " ymax=1.1,\n", + " ls={'Si reference': 'dashed'},\n", + " marker='none',\n", + " color={'Si reference': 'black'},\n", + " title='FREIA critical edge',\n", + " ylabel='Reflectivity',\n", + ")\n", + "critical_edge" + ] + }, + { + "cell_type": "markdown", + "id": "17", + "metadata": {}, "source": [ "Footprint correction is disabled with the widths above. Background subtraction and Q-resolution averaging are not included, and the reference is shown without resolution broadening. These effects and the counting statistics should be considered when comparing the curves." ] From 028aed359a71aa74e427ac1007e7a425db39fac4 Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Thu, 17 Sep 2026 16:17:57 +0200 Subject: [PATCH 16/16] docs: remove unnecessary comment --- .../docs/user-guide/freia/freia-reflectivity.ipynb | 9 --------- 1 file changed, 9 deletions(-) diff --git a/packages/essreflectometry/docs/user-guide/freia/freia-reflectivity.ipynb b/packages/essreflectometry/docs/user-guide/freia/freia-reflectivity.ipynb index 2e7af070b..98995fd6b 100644 --- a/packages/essreflectometry/docs/user-guide/freia/freia-reflectivity.ipynb +++ b/packages/essreflectometry/docs/user-guide/freia/freia-reflectivity.ipynb @@ -256,7 +256,6 @@ " sc.array(dims=['Q'], values=r, unit='dimensionless'),\n", " coords={'Q': sc.array(dims=['Q'], values=q, unit='1/angstrom')},\n", ")\n", - "# Keep the Q bin edges so the reduced curves are drawn as histograms.\n", "plot_curves = {}\n", "for label, curve in reflectivities.items():\n", " valid = sc.isfinite(curve.data) & sc.isfinite(sc.variances(curve.data))\n", @@ -302,14 +301,6 @@ ")\n", "critical_edge" ] - }, - { - "cell_type": "markdown", - "id": "17", - "metadata": {}, - "source": [ - "Footprint correction is disabled with the widths above. Background subtraction and Q-resolution averaging are not included, and the reference is shown without resolution broadening. These effects and the counting statistics should be considered when comparing the curves." - ] } ], "metadata": {