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..7130a0535 --- /dev/null +++ b/packages/essreflectometry/docs/user-guide/freia/freia-mcstas-visualization.ipynb @@ -0,0 +1,219 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "0", + "metadata": {}, + "source": [ + "# FREIA detector data\n", + "\n", + "Visualize detector images, arrival-time distributions, and wavelength spectra. Wavelengths are reconstructed using the WFM chopper settings." + ] + }, + { + "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" + ] + }, + { + "cell_type": "markdown", + "id": "2", + "metadata": {}, + "source": [ + "## Select a run\n", + "\n", + "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." + ] + }, + { + "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_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')" + ] + }, + { + "cell_type": "markdown", + "id": "4", + "metadata": {}, + "source": [ + "## Inspect the WFM chopper cascade\n", + "\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]]`." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5", + "metadata": {}, + "outputs": [], + "source": [ + "choppers = freia_mcstas.compute(DiskChoppers[SampleRun])\n", + "sc.DataGroup(choppers)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6", + "metadata": {}, + "outputs": [], + "source": [ + "frames = freia_mcstas.compute(ChopperFrameSequence[SampleRun])\n", + "frames.draw()" + ] + }, + { + "cell_type": "markdown", + "id": "7", + "metadata": {}, + "source": [ + "## Load and unwrap the detector events\n", + "\n", + "Compute wavelengths from the event arrival times and chopper transmission bands." + ] + }, + { + "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" + ] + }, + { + "cell_type": "markdown", + "id": "9", + "metadata": {}, + "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." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "10", + "metadata": {}, + "outputs": [], + "source": [ + "detector_image = raw.hist(longitude=120, height=64, dim=raw.dims)\n", + "pp.plot(detector_image, norm='log', title='FREIA detector', vmin=1e-1)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "11", + "metadata": {}, + "outputs": [], + "source": [ + "pulse_period = freia_mcstas.compute(PulsePeriod).to(unit='s').value\n", + "arrival_times = raw.hist(\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[\n", + " 'event_time_offset'\n", + "].to(unit='ms')\n", + "pp.plot(arrival_times, title='Arrival time within the source period')" + ] + }, + { + "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." + ] + }, + { + "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')" + ] + }, + { + "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', vmin=1e0)" + ] + } + ], + "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" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} 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..98995fd6b --- /dev/null +++ b/packages/essreflectometry/docs/user-guide/freia/freia-reflectivity.ipynb @@ -0,0 +1,327 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "0", + "metadata": {}, + "source": [ + "# FREIA reflectivity reduction\n", + "\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 downloads and caches WFM sample and direct-beam runs using `pooch` (`pip install pooch`)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1", + "metadata": {}, + "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", + "from ess.freia.corrections import RunNormalization\n", + "from ess.freia.types import (\n", + " DetectorRegionOfInterest,\n", + " IncidentMonitor,\n", + " WavelengthMonitor,\n", + ")\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", + " Sample,\n", + " SampleRun,\n", + " SampleSize,\n", + " WavelengthBins,\n", + " WavelengthDetector,\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "2", + "metadata": {}, + "source": [ + "## Select the runs\n", + "\n", + "`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" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3", + "metadata": {}, + "outputs": [], + "source": [ + "workflow = FreiaMcStasWorkflow(run_norm=RunNormalization.monitor_histogram)\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", + "workflow[QBins] = sc.geomspace('Q', 0.003, 0.5, 151, unit='1/angstrom')" + ] + }, + { + "cell_type": "markdown", + "id": "4", + "metadata": {}, + "source": [ + "Load the events, reconstruct wavelengths, and load the monitors once. Reuse these inputs when selecting each beam." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5", + "metadata": {}, + "outputs": [], + "source": [ + "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" + ] + }, + { + "cell_type": "markdown", + "id": "6", + "metadata": {}, + "source": [ + "## 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", + "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." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7", + "metadata": {}, + "outputs": [], + "source": [ + "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)" + ] + }, + { + "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 selected detector data." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9", + "metadata": {}, + "outputs": [], + "source": [ + "sc.plot(\n", + " {\n", + " 'Sample monitor': inputs[WavelengthMonitor[SampleRun]],\n", + " 'Direct-beam monitor': inputs[WavelengthMonitor[ReferenceRun]],\n", + " },\n", + " title='Incident wavelength spectra',\n", + " norm='log',\n", + " vmin=1e0,\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": "11", + "metadata": {}, + "outputs": [], + "source": [ + "# 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": "12", + "metadata": {}, + "source": [ + "## Compute one reflectivity curve per beam\n", + "\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", + "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": "13", + "metadata": {}, + "outputs": [], + "source": [ + "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[Sample] = reduction[ReducibleData[SampleRun]]\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": "14", + "metadata": {}, + "source": [ + "## Compare with the silicon reference\n", + "\n", + "`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": "15", + "metadata": {}, + "outputs": [], + "source": [ + "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", + " coords={'Q': sc.array(dims=['Q'], values=q, unit='1/angstrom')},\n", + ")\n", + "plot_curves = {}\n", + "for label, curve in reflectivities.items():\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", + " 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", + " title='FREIA reflectivity by beam',\n", + " ylabel='Reflectivity',\n", + ")\n", + "comparison" + ] + }, + { + "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" + ] + } + ], + "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" + } + }, + "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..1fd5454e7 --- /dev/null +++ b/packages/essreflectometry/docs/user-guide/freia/freia-wavelength-lookup-table.ipynb @@ -0,0 +1,164 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "0", + "metadata": {}, + "source": [ + "# FREIA analytical wavelength lookup table\n", + "\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." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1", + "metadata": {}, + "outputs": [], + "source": [ + "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", + " TimeResolution,\n", + ")\n", + "from ess.reflectometry.types import SampleRun" + ] + }, + { + "cell_type": "markdown", + "id": "2", + "metadata": {}, + "source": [ + "## Use the WFM configuration\n", + "\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)`." + ] + }, + { + "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.1, unit='m')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4", + "metadata": {}, + "outputs": [], + "source": [ + "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(\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}')" + ] + }, + { + "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 wavelength bounds." + ] + }, + { + "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(\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" + ] + } + ], + "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..fdbc1b1e0 100644 --- a/packages/essreflectometry/docs/user-guide/freia/index.md +++ b/packages/essreflectometry/docs/user-guide/freia/index.md @@ -1,3 +1,14 @@ # FREIA -FREIA-specific reduction guides will be added here once example workflows are available. +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 + +freia-mcstas-visualization +freia-wavelength-lookup-table +freia-reflectivity +``` diff --git a/packages/essreflectometry/pyproject.toml b/packages/essreflectometry/pyproject.toml index 2a564e2f9..02c852a5c 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/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 = ( diff --git a/packages/essreflectometry/src/ess/freia/__init__.py b/packages/essreflectometry/src/ess/freia/__init__.py index 3860491e8..ef3b9237b 100644 --- a/packages/essreflectometry/src/ess/freia/__init__.py +++ b/packages/essreflectometry/src/ess/freia/__init__.py @@ -3,7 +3,15 @@ import importlib.metadata from ..reflectometry import supermirror -from . import conversions, load, maskings, normalization, orso, resolution, workflow +from . import ( + conversions, + maskings, + mcstas, + normalization, + orso, + resolution, + workflow, +) from .types import ( AngularResolution, SampleSizeResolution, @@ -42,8 +50,8 @@ "SampleSizeResolution", "WavelengthResolution", "conversions", - "load", "maskings", + "mcstas", "normalization", "orso", "resolution", diff --git a/packages/essreflectometry/src/ess/freia/beamline.py b/packages/essreflectometry/src/ess/freia/beamline.py deleted file mode 100644 index e6ef59fa1..000000000 --- a/packages/essreflectometry/src/ess/freia/beamline.py +++ /dev/null @@ -1,26 +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 - -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=[], - **kwargs, - ) - workflow[DetectorBankSizes] = DETECTOR_BANK_SIZES - return workflow diff --git a/packages/essreflectometry/src/ess/freia/conversions.py b/packages/essreflectometry/src/ess/freia/conversions.py index 35ef06657..70e6e5a53 100644 --- a/packages/essreflectometry/src/ess/freia/conversions.py +++ b/packages/essreflectometry/src/ess/freia/conversions.py @@ -1,171 +1,102 @@ # 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 graph, tof from scippnexus import NXsample, NXsource from ..reflectometry.conversions import reflectometry_q from ..reflectometry.types import ( CoordTransformationGraph, - DetectorRotation, RunType, - SampleRotation, + SampleRun, ) +from .types import SampleSurfaceNormal -def reflectometry_q_x( - wavelength: sc.Variable, theta: sc.Variable, sample_rotation: sc.Variable +def outgoing_direction( + scattered_beam: sc.Variable, + wavelength: sc.Variable, + gravity: 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``. - - Source: - `Frédéric Ott, "Off-specular data representations in neutron reflectivity" `_ + """Unit direction of the outgoing ray at the sample, corrected for gravity. - Parameters - ---------- - wavelength: - Wavelength values for the events. - theta: - Angle of reflection for the events. - sample_rotation: - Angle of incidence. - - Returns - ------- - : - Qx-values. + Approximate the flight time using the straight sample-to-detector distance, + as in ScippNeutron's gravity correction. """ - 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 + 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 ) + return outgoing_beam / sc.norm(outgoing_beam) -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' - ) +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 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. +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. """ - 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' - ) + normal = sample_surface_normal / sc.norm(sample_surface_normal) + return sc.asin(sc.dot(outgoing_direction, normal)) 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 the scattering coordinates shared by sample and direct-beam runs.""" + 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, + **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, + 'gravity': lambda: gravity, } +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: 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, - ) + """Add the scattering coordinates provided by the run's transformation graph.""" + return da.transform_coords(rename_dims=False, **graph) -providers = (coordinate_transformation_graph,) +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 a6f5c1238..c9ae62231 100644 --- a/packages/essreflectometry/src/ess/freia/corrections.py +++ b/packages/essreflectometry/src/ess/freia/corrections.py @@ -7,54 +7,43 @@ from ..reflectometry import corrections as common_corrections from ..reflectometry.corrections import RunNormalization from ..reflectometry.types import ( - BeamDivergenceLimits, + BeamSize, CoordTransformationGraph, - CorrectionsToApply, + CorrectedDetector, ReducibleData, RunType, - RunUnnormalizedData, + Sample, + SampleRun, + SampleSize, WavelengthBins, WavelengthDetector, - YIndexLimits, - ZIndexLimits, ) from .conversions import add_coords from .maskings import add_masks -from .types import WavelengthMonitor +from .types import DetectorRegionOfInterest, 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, ) -> 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, @@ -64,7 +53,7 @@ def normalize_by_monitor_histogram( def normalize_by_monitor_integrated( - detector: RunUnnormalizedData[RunType], + detector: CorrectedDetector[RunType], *, monitor: WavelengthMonitor[RunType], uncertainty_broadcast_mode: UncertaintyBroadcastMode, @@ -89,33 +78,28 @@ 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 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], + beam_size: BeamSize[SampleRun], + sample_size: SampleSize[SampleRun], +) -> 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. + The correction applies only to the reflected run. + """ + 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.') + 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)) -default_corrections = {correct_by_footprint} -providers = (add_coords_masks_and_apply_corrections,) +providers = (add_coords_and_masks, prepare_sample) diff --git a/packages/essreflectometry/src/ess/freia/data.py b/packages/essreflectometry/src/ess/freia/data.py index 7c92c113e..44b91e458 100644 --- a/packages/essreflectometry/src/ess/freia/data.py +++ b/packages/essreflectometry/src/ess/freia/data.py @@ -1,10 +1,42 @@ +# 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 + _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", + }, ) -__all__ = [] + +def freia_mcstas_sample_run() -> Filename[SampleRun]: + """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 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", + "freia_mcstas_silicon_reflectivity", +] 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/maskings.py b/packages/essreflectometry/src/ess/freia/maskings.py index 472e8bc93..a0a60483d 100644 --- a/packages/essreflectometry/src/ess/freia/maskings.py +++ b/packages/essreflectometry/src/ess/freia/maskings.py @@ -1,54 +1,31 @@ # 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, -) - - -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, + roi: dict, + wavelength_bins: sc.Variable, ) -> 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', - ), - bdlim[1].to( - unit=da.coords["divergence_angle"].unit, - dtype='float64', - ), + """Mask events outside the ROI and wavelength range.""" + 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 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)) ), ) - da = da.bins.assign_masks( - wavelength=_not_between( - da.bins.coords['wavelength'], - wbins[0], - wbins[-1], - ), - ) - return da - - -providers = () diff --git a/packages/essreflectometry/src/ess/freia/mcstas.py b/packages/essreflectometry/src/ess/freia/mcstas.py index 4c48502ca..17624dd6e 100644 --- a/packages/essreflectometry/src/ess/freia/mcstas.py +++ b/packages/essreflectometry/src/ess/freia/mcstas.py @@ -1,3 +1,451 @@ -# 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 = () +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, + NeXusName, + Position, + RawDetector, + RunType, +) +from ess.reduce.unwrap import PulsePeriod +from scippneutron.chopper import DiskChopper + +from .types import IncidentMonitor, SampleSurfaceNormal, WavelengthMonitor + + +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], + }, +} + +# 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. + + 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 _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( + 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 settings.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 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. + 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') + 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.' + ) + variables = ['p', 't', 'id'] + if 'L' in data.get_component_variables(detector_name): + variables.append('L') + values = data.get_event_data( + variables=variables, 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'), + }, + ) + 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( + filename: Filename[RunType], + detector_name: NeXusDetectorName, + geometry: EmptyDetector[RunType], + pulse_period: PulsePeriod, +) -> RawDetector[RunType]: + """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) + ) + + +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 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. + + 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) + # 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: + # 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}.') + 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_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..2faf0a962 100644 --- a/packages/essreflectometry/src/ess/freia/normalization.py +++ b/packages/essreflectometry/src/ess/freia/normalization.py @@ -1,3 +1,55 @@ -# 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 -providers = () +from ..reflectometry.conversions import reflectometry_q +from ..reflectometry.types import ( + QBins, + ReducibleData, + Reference, + ReferenceRun, + ReflectivityOverQ, + Sample, +) +from .conversions import theta + + +def evaluate_direct_beam( + direct_beam: ReducibleData[ReferenceRun], +) -> Reference: + """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, reflection_angle)) + ) + + +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) + ) + return ReflectivityOverQ( + (numerator / denominator.data).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..cf49d7534 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,16 @@ # 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 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 ``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 327a5f9d9..eec0babd1 100644 --- a/packages/essreflectometry/src/ess/freia/workflow.py +++ b/packages/essreflectometry/src/ess/freia/workflow.py @@ -3,73 +3,45 @@ 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 from ..reflectometry.types import ( - BeamDivergenceLimits, - CorrectionsToApply, DetectorSpatialResolution, LookupTableRelativeErrorThreshold, NeXusDetectorName, - RunType, - SampleRotationOffset, -) -from . import ( - beamline, - conversions, - corrections, - load, - maskings, - mcstas, - normalization, - orso, + ReferenceRun, + SampleRun, ) +from . import conversions, corrections, mcstas, normalization, orso from .corrections import RunNormalization, insert_run_normalization +from .types import IncidentMonitor -_general_providers = ( +providers = ( *reflectometry_providers, *conversions.providers, *corrections.providers, - *maskings.providers, *normalization.providers, *orso.providers, - *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", - 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, + return default_parameters() | { + NeXusDetectorName: "Multiblade", LookupTableRelativeErrorThreshold: { - "detector": 0.06, + "Multiblade": 0.06, }, - UncertaintyBroadcastMode: UncertaintyBroadcastMode.drop, } @@ -77,8 +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'), - CorrectionsToApply: corrections.default_corrections, + DetectorBankSizes: { + "multiblade_detector": {"strip": 64, "blade": 32, "wire": 32}, + }, DetectorSpatialResolution: 0.0025 * sc.units.m, LookupTableRelativeErrorThreshold: { "multiblade_detector": float('inf'), @@ -90,10 +63,13 @@ 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 reducing FREIA McStas events with a no-sample direct beam. + + Loads geometry and uses the default WFM chopper settings. Reduction inputs + and outputs are described in :func:`FreiaWorkflow`. Parameters ---------- @@ -104,10 +80,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 @@ -121,6 +98,21 @@ def FreiaWorkflow( ) -> sciline.Pipeline: """Workflow for reduction of data for the Freia instrument. + 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. + 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]``. + Parameters ---------- run_norm: @@ -130,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/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) diff --git a/packages/essreflectometry/tests/freia/conversions_test.py b/packages/essreflectometry/tests/freia/conversions_test.py new file mode 100644 index 000000000..094af42bf --- /dev/null +++ b/packages/essreflectometry/tests/freia/conversions_test.py @@ -0,0 +1,54 @@ +# 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 outgoing_direction, scattering_angle, 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], [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'), + ) + + assert_allclose( + 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_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') + ) + direction = outgoing_direction( + 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'), + ) + + assert_allclose( + 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 new file mode 100644 index 000000000..bb90bd485 --- /dev/null +++ b/packages/essreflectometry/tests/freia/mcstas_test.py @@ -0,0 +1,109 @@ +# SPDX-License-Identifier: BSD-3-Clause +# Copyright (c) 2026 Scipp contributors (https://github.com/scipp) +from pathlib import Path + +import h5py +import numpy as np +import pytest +import scipp as sc +from ess.reduce.nexus.types import Filename, NeXusName, RawDetector +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 + + +def _component(components, name, position): + group = components.create_group(f'{len(components):04d}_{name}') + group['Position'] = position + group['Rotation'] = np.eye(3) + return group + + +def test_load_detector_and_compute_q(): + workflow = FreiaMcStasWorkflow() + workflow[Filename[SampleRun]] = freia_mcstas_sample_run() + workflow[DetectorRegionOfInterest[SampleRun]] = {} + workflow[WavelengthBins] = sc.array( + dims=['wavelength'], values=[1.0, 12.0], unit='angstrom' + ) + + result = workflow.compute((RawDetector[SampleRun], CorrectedDetector[SampleRun])) + + raw = result[RawDetector[SampleRun]] + assert raw.sizes == {'pixel_id': 2048 * 64} + events = raw.bins.constituents['data'] + assert events.sizes == {'event': 555874} + assert_allclose( + 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).any().value + + +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(filename, 'Multiblade') + + assert_allclose( + 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', + ), + ) + + +@pytest.mark.parametrize( + 'filename', [freia_mcstas_sample_run, freia_mcstas_reference_run] +) +def test_load_monitor(filename): + path = filename() + workflow = FreiaMcStasWorkflow() + workflow[Filename[SampleRun]] = path + workflow[NeXusName[IncidentMonitor]] = 'SampleLambda' + + result = workflow.compute(WavelengthMonitor[SampleRun]) + + assert_identical( + 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=monitor['data'][:], + variances=monitor['errors'][:] ** 2, + unit='counts', + ), + ) diff --git a/packages/essreflectometry/tests/freia/workflow_test.py b/packages/essreflectometry/tests/freia/workflow_test.py index 41bc7a887..3c2ac20cf 100644 --- a/packages/essreflectometry/tests/freia/workflow_test.py +++ b/packages/essreflectometry/tests/freia/workflow_test.py @@ -1,41 +1,118 @@ # SPDX-License-Identifier: BSD-3-Clause # Copyright (c) 2026 Scipp contributors (https://github.com/scipp) -import inspect +import numpy as np +import scipp as sc +from ess.reduce.nexus.types import GravityVector, Position +from scipp.testing import assert_allclose, assert_identical +from scippnexus import NXsample, NXsource -from ess.reduce import workflow as reduce_workflow - -from ess import freia +from ess.freia import FreiaWorkflow from ess.freia.corrections import RunNormalization -from ess.reflectometry.types import CorrectionsToApply, NeXusDetectorName +from ess.freia.types import ( + DetectorRegionOfInterest, + SampleSurfaceNormal, + WavelengthMonitor, +) +from ess.reflectometry.corrections import footprint_on_sample +from ess.reflectometry.types import ( + BeamSize, + QBins, + ReducibleData, + ReferenceRun, + ReflectivityOverQ, + Sample, + SampleRun, + SampleSize, + WavelengthBins, + WavelengthDetector, +) + +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( + 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') + # 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, 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') + # 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]] = { + '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' + ) + wf[QBins] = sc.array(dims=['Q'], values=[1.0, 2.0, 4.0, 6.0], unit='1/angstrom') + return wf -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_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)}, + ) + 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 -def test_freia_workflows_have_expected_run_normalization_defaults(): - assert ( - inspect.signature(freia.FreiaMcStasWorkflow).parameters["run_norm"].default - is RunNormalization.none + 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 ( - inspect.signature(freia.FreiaWorkflow).parameters["run_norm"].default - is RunNormalization.proton_charge + assert_allclose(result['Q', :2].data, expected) + assert_identical( + result.masks['direct_beam'], + sc.array(dims=['Q'], values=[False, False, True]), ) -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 test_rebinning_integrates_before_dividing(): + wf = _make_workflow(RunNormalization.none) + wf[Sample] = wf[ReducibleData[SampleRun]] + wf[QBins] = sc.array(dims=['Q'], values=[1.0, 4.0], unit='1/angstrom') + + result = wf.compute(ReflectivityOverQ) + + assert_allclose( + result.data, + sc.array( + dims=['Q'], values=[60.0 / 200.0], variances=[0.3**2 * (1 / 60 + 1 / 200)] + ), + ) 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