From db01fec29de1106c43133c02685c7755449fbf76 Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Tue, 2 Dec 2025 15:02:47 +0100 Subject: [PATCH 01/20] feat: add notebook for loading mcstas data into CODA nexus file --- src/ess/estia/load.py | 29 +++- src/ess/estia/mcstas.py | 42 ++++-- src/ess/estia/workflow.py | 2 +- tests/estia/mcstas_data_test.py | 2 - tools/mcstas_to_nexus.ipynb | 250 ++++++++++++++++++++++++++++++++ 5 files changed, 302 insertions(+), 23 deletions(-) create mode 100644 tools/mcstas_to_nexus.ipynb diff --git a/src/ess/estia/load.py b/src/ess/estia/load.py index 40cd5227..feb3b7ad 100644 --- a/src/ess/estia/load.py +++ b/src/ess/estia/load.py @@ -1,5 +1,7 @@ # SPDX-License-Identifier: BSD-3-Clause # Copyright (c) 2025 Scipp contributors (https://github.com/scipp) +import pathlib + import h5py import scipp as sc @@ -13,16 +15,29 @@ from .mcstas import parse_events_ascii, parse_events_h5 -def load_mcstas_events( +def load_mcstas_provider( filename: Filename[RunType], sample_rotation_offset: SampleRotationOffset[RunType], ) -> RawDetector[RunType]: + """See :py:`load_mcstas`.""" + da = load_mcstas(filename) + da.coords['sample_rotation'] += sample_rotation_offset.to( + unit=da.coords['sample_rotation'].unit + ) + return RawDetector[RunType](da) + + +def load_mcstas( + filename: str | pathlib.Path | sc.DataArray, +) -> sc.DataArray: """ Load event data from a McStas run and reshape it to look like what we would expect if we loaded data from a real experiment file. """ - if h5py.is_hdf5(filename): + if isinstance(filename, sc.DataArray): + da = filename + elif h5py.is_hdf5(filename): with h5py.File(filename) as f: da = parse_events_h5(f) else: @@ -35,9 +50,6 @@ def load_mcstas_events( da.coords['detector_rotation'] = 2 * da.coords['sample_rotation'] + sc.scalar( 1.65, unit='deg' ) - da.coords['sample_rotation'] += sample_rotation_offset.to( - unit=da.coords['sample_rotation'].unit - ) nblades = DETECTOR_BANK_SIZES['multiblade_detector']['blade'] nwires = DETECTOR_BANK_SIZES['multiblade_detector']['wire'] @@ -79,7 +91,9 @@ def load_mcstas_events( da.bins.coords['event_time_offset'] = ( sc.scalar(1, unit='s') * da.bins.coords['t'] ).to(unit='ns') % sc.scalar(1 / 14, unit='s').to(unit='ns') - da.bins.coords['wavelength'] = sc.scalar(1, unit='angstrom') * da.bins.coords['L'] + da.bins.coords['wavelength_from_mcstas'] = ( + sc.scalar(1.0, unit='angstrom') * da.bins.coords['L'] + ) da.coords["sample_size"] = sc.scalar(1.0, unit='m') * float( da.coords['sample_length'].value @@ -96,8 +110,7 @@ def load_mcstas_events( ) da.bins.coords.pop('L') da.bins.coords.pop('t') - da.bins.coords['wavelength_from_mcstas'] = da.bins.coords.pop('wavelength') - return RawDetector[RunType](da) + return da providers = () diff --git a/src/ess/estia/mcstas.py b/src/ess/estia/mcstas.py index f1a07f99..01bb9885 100644 --- a/src/ess/estia/mcstas.py +++ b/src/ess/estia/mcstas.py @@ -1,5 +1,6 @@ # Copyright (c) 2025 Scipp contributors (https://github.com/scipp) import h5py +import numpy as np import pandas as pd import scipp as sc @@ -63,7 +64,7 @@ def parse_events_ascii(lines): raise ValueError('Could not parse the file as a list of events.') -def parse_events_h5(f): +def parse_events_h5(f, nevents_to_sample=None): if isinstance(f, str): with h5py.File(f) as ff: return parse_events_h5(ff) @@ -75,20 +76,37 @@ def parse_events_h5(f): 'NXentry/simulation/Param', ) events = events[()] - da = sc.DataArray( - # The squares on the variances is the correct way to load weighted event data. - # Consult the McStas documentation - # (section 2.2.1) https://www.mcstas.org/documentation/manual/ - # for more information. - sc.array(dims=['events'], values=events[:, 0], variances=events[:, 0] ** 2), - ) + if nevents_to_sample is None: + da = sc.DataArray( + # The squares on the variances is the correct way to load + # weighted event data. + # Consult the McStas documentation + # (section 2.2.1) https://www.mcstas.org/documentation/manual/ + # for more information. + sc.array(dims=['events'], values=events[:, 0], variances=events[:, 0] ** 2), + ) + for i, label in enumerate(data.attrs["ylabel"].decode().strip().split(' ')): + if label == 'p': + continue + da.coords[label] = sc.array(dims=['events'], values=events[:, i]) + else: + probabilities = events[:, 0] + inds = np.random.choice( + np.arange(len(probabilities)), + nevents_to_sample, + p=probabilities / probabilities.sum(), + ) + da = sc.DataArray( + sc.ones(dims=['events'], shape=(nevents_to_sample,), with_variances=True), + ) + for i, label in enumerate(data.attrs["ylabel"].decode().strip().split(' ')): + if label == 'p': + continue + da.coords[label] = sc.array(dims=['events'], values=events[inds, i]) + for name, value in data.attrs.items(): da.coords[name] = sc.scalar(value.decode()) - for i, label in enumerate(data.attrs["ylabel"].decode().strip().split(' ')): - if label == 'p': - continue - da.coords[label] = sc.array(dims=['events'], values=events[:, i]) for k, v in params.items(): v = v[0] if isinstance(v, bytes): diff --git a/src/ess/estia/workflow.py b/src/ess/estia/workflow.py index c1bcf5a2..7b817768 100644 --- a/src/ess/estia/workflow.py +++ b/src/ess/estia/workflow.py @@ -27,7 +27,7 @@ mcstas_providers = ( *_general_providers, *load.providers, - load.load_mcstas_events, + load.load_mcstas_provider, ) """List of providers for setting up a Sciline pipeline for McStas data. diff --git a/tests/estia/mcstas_data_test.py b/tests/estia/mcstas_data_test.py index acd586af..46064e05 100644 --- a/tests/estia/mcstas_data_test.py +++ b/tests/estia/mcstas_data_test.py @@ -18,7 +18,6 @@ estia_mcstas_sample_run, estia_tof_lookup_table, ) -from ess.estia.load import load_mcstas_events from ess.reflectometry import orso from ess.reflectometry.types import ( BeamDivergenceLimits, @@ -39,7 +38,6 @@ @pytest.fixture def estia_mcstas_pipeline() -> sciline.Pipeline: wf = EstiaMcStasWorkflow() - wf.insert(load_mcstas_events) wf[Filename[ReferenceRun]] = estia_mcstas_reference_run() wf[YIndexLimits] = sc.scalar(35), sc.scalar(64) diff --git a/tools/mcstas_to_nexus.ipynb b/tools/mcstas_to_nexus.ipynb new file mode 100644 index 00000000..0ec6337a --- /dev/null +++ b/tools/mcstas_to_nexus.ipynb @@ -0,0 +1,250 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "0", + "metadata": {}, + "source": [ + "# McStas to NeXus for ESTIA\n", + "\n", + "This notebook converts data from a McStas simulation output (`.h5` filetype) to a NeXus file that uses a file for the Estia instrument (written by CODA) as a template for the geometry information.\n", + "\n", + "It adds\n", + "\n", + "* events to the `multiblade_detector` detector,\n", + "* a single value to the `proton_current` log,\n", + "* a single detector rotation value to the `detector_rotation` log,\n", + "* a single sample rotation value to the `sample_rotation` log." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import scipp as sc\n", + "import h5py\n", + "import shutil\n", + "\n", + "from ess.estia.mcstas import parse_events_h5\n", + "from ess.estia.load import load_mcstas\n", + "from ess.estia.data import estia_mcstas_example\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2", + "metadata": {}, + "outputs": [], + "source": [ + "def replace_dataset(entry, name, values):\n", + " attrs = dict(entry[name].attrs)\n", + " del entry[name]\n", + " dset = entry.create_dataset(name, data=values)\n", + " dset.attrs.update(attrs)\n", + "\n", + "\n", + "def mcstas_to_nexus(\n", + " mcstas_data_file: str,\n", + " template_nexus_file: str,\n", + " outfile: str,\n", + " nevents: int = 1_000_000,\n", + "):\n", + " \"\"\"\n", + " Store the events from a McStas Estia simulation in a NeXus CODA file.\n", + "\n", + " Parameters\n", + " ----------\n", + " mcstas_data_file:\n", + " Data file containing simulated McStas events.\n", + " template_nexus_file:\n", + " NeXus file containing geometry and instrument info, used as a template.\n", + " outfile:\n", + " Output file to be written.\n", + " nevents:\n", + " Number of events to have in the output file\n", + " (events are sampled from the probabilities of the mcstas events).\n", + " \"\"\"\n", + " detector_entry_path = '/entry/instrument/multiblade_detector'\n", + "\n", + " with h5py.File(template_nexus_file, \"r\") as f:\n", + " det_numbers = f[f\"{detector_entry_path}/detector_number\"][()]\n", + "\n", + " with h5py.File(mcstas_data_file) as f:\n", + " events = parse_events_h5(f, nevents_to_sample=nevents)\n", + " binned = load_mcstas(events).transpose(['blade', 'wire', 'strip'])\n", + "\n", + " toa = (binned.bins.coords[\"event_time_zero\"] + binned.bins.coords[\"event_time_offset\"]).bins.concat().value\n", + "\n", + " # IMPORTANT! we need to sort the arrays below according to toa,\n", + " # so that the event_index does not get messed up!\n", + " event_id = sc.sort(\n", + " (\n", + " sc.bins_like(binned, sc.array(dims=binned.dims, values=det_numbers.reshape(binned.shape)))\n", + " .bins.concat()\n", + " .value\n", + " ),\n", + " key=toa,\n", + " )\n", + " event_time_zero = sc.sort(binned.bins.coords[\"event_time_zero\"].bins.concat().value, key=toa)\n", + " event_time_offset = sc.sort(binned.bins.coords[\"event_time_offset\"].bins.concat().value, key=toa)\n", + "\n", + " event_index = sc.DataArray(\n", + " data=sc.ones_like(event_time_offset),\n", + " coords={\"event_time_zero\": event_time_zero},\n", + " ).group(\"event_time_zero\")\n", + "\n", + " event_index = sc.cumsum(event_index.bins.size())\n", + " event_index.values = np.concatenate([[0], event_index.values[:-1]])\n", + "\n", + " shutil.copyfile(template_nexus_file, outfile)\n", + "\n", + " with h5py.File(outfile, \"r+\") as f:\n", + "\n", + " detector_rotation = f[f\"{detector_entry_path}/transformations/detector_rotation\"]\n", + " replace_dataset(\n", + " detector_rotation,\n", + " name=\"value\",\n", + " values=[\n", + " binned\n", + " .coords['detector_rotation']\n", + " .to(unit=detector_rotation['value'].attrs['units'])\n", + " .value\n", + " ]\n", + " )\n", + " replace_dataset(\n", + " detector_rotation,\n", + " name=\"time\",\n", + " values=[\n", + " event_time_zero\n", + " .min()\n", + " .to(unit=detector_rotation['time'].attrs['units'])\n", + " .value\n", + " .astype('uint64')\n", + " ]\n", + " )\n", + "\n", + " # TODO: Put sample_rotation below \"transformations\" group\n", + " # TODO: Set correct nexus classes on the sample_rotation group\n", + " sample = f[\"entry/sample\"]\n", + " sample_rotation = sample.create_group('sample_rotation')\n", + " sample_rotation.create_dataset(\n", + " 'value',\n", + " data=[\n", + " binned\n", + " .coords['sample_rotation']\n", + " .to(unit='deg')\n", + " .value\n", + " .astype('float64')\n", + " ],\n", + " )\n", + " sample_rotation['value'].attrs['units'] = 'degrees'\n", + " sample_rotation.create_dataset(\n", + " 'time',\n", + " data=[\n", + " event_time_zero\n", + " .min()\n", + " .to(unit='ns')\n", + " .value\n", + " .astype('uint64')\n", + " ],\n", + " )\n", + " sample_rotation['time'].attrs['units'] = 'ns'\n", + "\n", + " f[\"entry/neutron_prod_info/pulse_charge/value\"][:] = 1.\n", + "\n", + " event_data = f[f\"{detector_entry_path}/event_data\"]\n", + " replace_dataset(event_data, name=\"event_id\", values=event_id.values.astype('uint32'))\n", + " replace_dataset(\n", + " event_data,\n", + " name=\"event_time_offset\",\n", + " values=event_time_offset.to(\n", + " unit=event_data[\"event_time_offset\"].attrs[\"units\"], copy=False\n", + " ).values.astype('uint32'),\n", + " )\n", + " replace_dataset(event_data, name=\"event_index\", values=event_index.values.astype('uint64'))\n", + " replace_dataset(\n", + " event_data,\n", + " name=\"event_time_zero\",\n", + " values=event_index.coords[\"event_time_zero\"]\n", + " .to(unit=event_data[\"event_time_zero\"].attrs[\"units\"], copy=False)\n", + " .values.astype('uint64'),\n", + " )" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3", + "metadata": {}, + "outputs": [], + "source": [ + "\n", + "mcstas_to_nexus(\n", + " mcstas_data_file=estia_mcstas_example('reference'),\n", + " template_nexus_file=\"~/Downloads/estia_999999_00008786.hdf\",\n", + " outfile='out.nx',\n", + " nevents=1000_000,\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4", + "metadata": {}, + "outputs": [], + "source": [ + "from ess.reflectometry.types import *\n", + "from ess.estia import EstiaWorkflow, EstiaMcStasWorkflow\n", + "\n", + "wfm = EstiaMcStasWorkflow()\n", + "wfm[Filename[SampleRun]] = estia_mcstas_example('reference')\n", + "h = sc.round(sc.values(wfm.compute(RawDetector[SampleRun]).flatten(['blade', 'wire'], 'x').bins.sum()))\n", + "p1 = (1e6 * h / h.sum().value).plot(norm='log', vmin=1, vmax=3e2)\n", + "\n", + "wf = EstiaWorkflow()\n", + "wf[Filename[SampleRun]] = 'out.nx'\n", + "p2 = wf.compute(RawDetector[SampleRun]).flatten(['blade', 'wire'], 'x').transpose().bins.sum().plot(norm='log', vmin=1, vmax=3e2)\n", + "\n", + "p1 + p2" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5", + "metadata": {}, + "outputs": [], + "source": [ + "wfm.compute(RawDetector[SampleRun]).bins.concat().hist(event_time_offset=200).plot() + wf.compute(RawDetector[SampleRun]).bins.concat().hist(event_time_offset=200).plot()" + ] + } + ], + "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.11.13" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} From be64ad4217c563298bc9b3d83f044574cfac4895 Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Wed, 3 Dec 2025 09:24:49 +0100 Subject: [PATCH 02/20] feat: Nexus file writer for McStas data --- tools/mcstas_to_nexus.ipynb | 33 ++++++++++++++++++++++++++++----- 1 file changed, 28 insertions(+), 5 deletions(-) diff --git a/tools/mcstas_to_nexus.ipynb b/tools/mcstas_to_nexus.ipynb index 0ec6337a..97068902 100644 --- a/tools/mcstas_to_nexus.ipynb +++ b/tools/mcstas_to_nexus.ipynb @@ -187,9 +187,9 @@ "\n", "mcstas_to_nexus(\n", " mcstas_data_file=estia_mcstas_example('reference'),\n", - " template_nexus_file=\"~/Downloads/estia_999999_00008786.hdf\",\n", + " template_nexus_file=\"/home/johannes/Downloads/estia_999999_00008786.hdf\",\n", " outfile='out.nx',\n", - " nevents=1000_000,\n", + " nevents=20_000_000,\n", ")" ] }, @@ -203,14 +203,16 @@ "from ess.reflectometry.types import *\n", "from ess.estia import EstiaWorkflow, EstiaMcStasWorkflow\n", "\n", + "n = 20_000_000\n", + "\n", "wfm = EstiaMcStasWorkflow()\n", "wfm[Filename[SampleRun]] = estia_mcstas_example('reference')\n", - "h = sc.round(sc.values(wfm.compute(RawDetector[SampleRun]).flatten(['blade', 'wire'], 'x').bins.sum()))\n", - "p1 = (1e6 * h / h.sum().value).plot(norm='log', vmin=1, vmax=3e2)\n", + "h = wfm.compute(RawDetector[SampleRun]).flatten(['blade', 'wire'], 'x').bins.sum()\n", + "p1 = sc.round(sc.values(0.9 * n * h / h.sum().value)).plot(norm='log')\n", "\n", "wf = EstiaWorkflow()\n", "wf[Filename[SampleRun]] = 'out.nx'\n", - "p2 = wf.compute(RawDetector[SampleRun]).flatten(['blade', 'wire'], 'x').transpose().bins.sum().plot(norm='log', vmin=1, vmax=3e2)\n", + "p2 = wf.compute(RawDetector[SampleRun]).flatten(['blade', 'wire'], 'x').transpose().bins.sum().plot(norm='log')\n", "\n", "p1 + p2" ] @@ -224,6 +226,27 @@ "source": [ "wfm.compute(RawDetector[SampleRun]).bins.concat().hist(event_time_offset=200).plot() + wf.compute(RawDetector[SampleRun]).bins.concat().hist(event_time_offset=200).plot()" ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6", + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib ipympl\n", + "sc.round(sc.values(0.9* n * h / h.sum().value)).plot(norm='log')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7", + "metadata": {}, + "outputs": [], + "source": [ + "wf.compute(RawDetector[SampleRun]).flatten(['blade', 'wire'], 'x').transpose().bins.sum().plot(norm='log')" + ] } ], "metadata": { From 7a45ed0287d2cc2e202ad0632c3e89f0f6963bf1 Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Thu, 4 Dec 2025 10:35:21 +0100 Subject: [PATCH 03/20] fix: make coord transformation RunType dependent (same as in esssans and essdiffraction) --- src/ess/amor/conversions.py | 5 +++- src/ess/amor/normalization.py | 3 +- src/ess/amor/workflow.py | 2 +- src/ess/estia/conversions.py | 42 ++++++++++++++++++++++------ src/ess/estia/corrections.py | 2 +- src/ess/estia/normalization.py | 3 +- src/ess/offspec/types.py | 4 --- src/ess/reflectometry/conversions.py | 16 ++++++++++- src/ess/reflectometry/types.py | 3 +- 9 files changed, 60 insertions(+), 20 deletions(-) diff --git a/src/ess/amor/conversions.py b/src/ess/amor/conversions.py index 29935802..a609c091 100644 --- a/src/ess/amor/conversions.py +++ b/src/ess/amor/conversions.py @@ -7,6 +7,7 @@ from ..reflectometry.types import ( BeamDivergenceLimits, CoordTransformationGraph, + RunType, WavelengthBins, YIndexLimits, ZIndexLimits, @@ -135,7 +136,9 @@ def wavelength( return out.to(unit='angstrom', copy=False) -def coordinate_transformation_graph(gravity: GravityToggle) -> CoordTransformationGraph: +def coordinate_transformation_graph( + gravity: GravityToggle, run_type: RunType +) -> CoordTransformationGraph[RunType]: return { "wavelength": wavelength, "theta": theta if gravity else theta_no_gravity, diff --git a/src/ess/amor/normalization.py b/src/ess/amor/normalization.py index 9c0b1d8a..405bae0d 100644 --- a/src/ess/amor/normalization.py +++ b/src/ess/amor/normalization.py @@ -14,6 +14,7 @@ ReducedReference, ReducibleData, Reference, + ReferenceRun, Sample, SampleRun, ) @@ -71,7 +72,7 @@ def evaluate_reference_at_sample_coords( reference: ReducedReference, sample: ReducibleData[SampleRun], detector_spatial_resolution: DetectorSpatialResolution[SampleRun], - graph: CoordTransformationGraph, + graph: CoordTransformationGraph[ReferenceRun], ) -> Reference: """ Adds a :math:`Q` and :math:`Q`-resolution coordinate to each bin of the ideal diff --git a/src/ess/amor/workflow.py b/src/ess/amor/workflow.py index 0c1dcdf3..7c5fa606 100644 --- a/src/ess/amor/workflow.py +++ b/src/ess/amor/workflow.py @@ -26,7 +26,7 @@ def add_coords_masks_and_apply_corrections( bdlim: BeamDivergenceLimits, wbins: WavelengthBins, proton_current: ProtonCurrent[RunType], - graph: CoordTransformationGraph, + graph: CoordTransformationGraph[RunType], ) -> ReducibleData[RunType]: """ Computes coordinates, masks and corrections that are diff --git a/src/ess/estia/conversions.py b/src/ess/estia/conversions.py index f68403b6..7746c29d 100644 --- a/src/ess/estia/conversions.py +++ b/src/ess/estia/conversions.py @@ -2,6 +2,9 @@ # Copyright (c) 2025 Scipp contributors (https://github.com/scipp) import scipp as sc from scippneutron.conversion.tof import wavelength_from_tof +from scippnexus import NXsample, NXsource + +from ess.reduce.nexus.types import DetectorBankSizes, Position from ..reflectometry.conversions import reflectometry_q from ..reflectometry.types import ( @@ -63,26 +66,47 @@ def divergence_angle( def detector_ltotal_from_raw( - da: RawDetector[RunType], graph: CoordTransformationGraph + da: RawDetector[RunType], graph: CoordTransformationGraph[RunType] ) -> DetectorLtotal[RunType]: - return da.transform_coords(['Ltotal'], graph=graph).coords['Ltotal'] - - -def coordinate_transformation_graph() -> CoordTransformationGraph: + return da.transform_coords( + ['Ltotal'], + graph=graph, + ).coords['Ltotal'] + + +def coordinate_transformation_graph( + source_position: Position[NXsource, RunType], + sample_position: Position[NXsample, RunType], + detector_bank_sizes: DetectorBankSizes, +) -> CoordTransformationGraph[RunType]: + bank = detector_bank_sizes['multiblade_detector'] return { "wavelength": wavelength_from_tof, "theta": theta, "divergence_angle": divergence_angle, "Q": reflectometry_q, "L1": lambda source_position, sample_position: sc.norm( - sample_position - source_position + sample_position - source_position.to(unit=sample_position.unit) ), # + extra correction for guides? - "L2": lambda position, sample_position: sc.norm(position - sample_position), - "Ltotal": lambda L1, L2: L1 + L2, + "L2": lambda position, sample_position: sc.norm( + position - sample_position.to(unit=position.unit) + ), + "Ltotal": lambda L1, L2: L1.to(unit=L2.unit) + L2, + 'sample_rotation': lambda: sc.scalar(1.0, unit='deg'), + 'detector_rotation': lambda: sc.scalar(3.65, unit='deg'), + 'source_position': lambda: source_position, + 'sample_position': lambda: sample_position, + 'blade': lambda: sc.arange('blade', 0, bank['blade']), + 'wire': lambda: sc.arange('wire', 0, bank['wire']), + 'strip': lambda: sc.arange('strip', 0, bank['strip']), + 'z_index': lambda blade, wire: blade * wire, + 'sample_size': lambda: sc.scalar(20.0, unit='mm'), } -def mcstas_wavelength_coordinate_transformation_graph() -> CoordTransformationGraph: +def mcstas_wavelength_coordinate_transformation_graph( + run_type: RunType, +) -> CoordTransformationGraph[RunType]: return { **coordinate_transformation_graph(), "wavelength": lambda wavelength_from_mcstas: wavelength_from_mcstas, diff --git a/src/ess/estia/corrections.py b/src/ess/estia/corrections.py index b5a93965..be30f19a 100644 --- a/src/ess/estia/corrections.py +++ b/src/ess/estia/corrections.py @@ -29,7 +29,7 @@ def add_coords_masks_and_apply_corrections( bdlim: BeamDivergenceLimits, wbins: WavelengthBins, proton_current: ProtonCurrent[RunType], - graph: CoordTransformationGraph, + graph: CoordTransformationGraph[RunType], ) -> ReducibleData[RunType]: """ Computes coordinates, masks and corrections that are diff --git a/src/ess/estia/normalization.py b/src/ess/estia/normalization.py index 939704e5..82c99def 100644 --- a/src/ess/estia/normalization.py +++ b/src/ess/estia/normalization.py @@ -14,6 +14,7 @@ ReducedReference, ReducibleData, Reference, + ReferenceRun, Sample, SampleRun, ) @@ -66,7 +67,7 @@ def mask_events_where_supermirror_does_not_cover( def evaluate_reference( reference: ReducedReference, sample: ReducibleData[SampleRun], - graph: CoordTransformationGraph, + graph: CoordTransformationGraph[ReferenceRun], detector_spatial_resolution: DetectorSpatialResolution[SampleRun], ) -> Reference: """ diff --git a/src/ess/offspec/types.py b/src/ess/offspec/types.py index 4bb4cf89..8f6e0966 100644 --- a/src/ess/offspec/types.py +++ b/src/ess/offspec/types.py @@ -14,10 +14,6 @@ BackgroundMinWavelength = NewType("BackgroundMinWavelength", sc.Variable) -class CoordTransformationGraph(sciline.Scope[RunType, dict], dict): - """Coordinate transformation for the runtype""" - - class MonitorData(sciline.Scope[RunType, sc.DataArray], sc.DataArray): """ "Monitor data from the run file, with background subtracted""" diff --git a/src/ess/reflectometry/conversions.py b/src/ess/reflectometry/conversions.py index 7c05b278..f72e5bb7 100644 --- a/src/ess/reflectometry/conversions.py +++ b/src/ess/reflectometry/conversions.py @@ -70,7 +70,21 @@ def add_coords( ) -> sc.DataArray: "Adds scattering coordinates to the raw detector data." return da.transform_coords( - ("wavelength", "theta", "divergence_angle", "Q", "L1", "L2"), + ( + "wavelength", + "theta", + "divergence_angle", + "Q", + "L1", + "L2", + "blade", + "wire", + "strip", + "z_index", + "sample_rotation", + "detector_rotation", + "sample_size", + ), graph, rename_dims=False, keep_intermediate=False, diff --git a/src/ess/reflectometry/types.py b/src/ess/reflectometry/types.py index 51c6cad3..905f48e0 100644 --- a/src/ess/reflectometry/types.py +++ b/src/ess/reflectometry/types.py @@ -29,7 +29,8 @@ DetectorLtotal = tof_t.DetectorLtotal -CoordTransformationGraph = NewType("CoordTransformationGraph", dict) +class CoordTransformationGraph(sciline.Scope[RunType, dict], dict): + pass class RawChopper(sciline.Scope[RunType, sc.DataGroup], sc.DataGroup): From 93ce3c97aa7191f88b942594ac2af4aa27dd81a5 Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Thu, 4 Dec 2025 10:35:51 +0100 Subject: [PATCH 04/20] WIP: writing mcstas to nexus --- tools/mcstas_to_nexus.ipynb | 132 ++++++++++++++++++++++++++++++++++-- 1 file changed, 125 insertions(+), 7 deletions(-) diff --git a/tools/mcstas_to_nexus.ipynb b/tools/mcstas_to_nexus.ipynb index 97068902..09881cd4 100644 --- a/tools/mcstas_to_nexus.ipynb +++ b/tools/mcstas_to_nexus.ipynb @@ -77,15 +77,17 @@ "\n", " with h5py.File(mcstas_data_file) as f:\n", " events = parse_events_h5(f, nevents_to_sample=nevents)\n", - " binned = load_mcstas(events).transpose(['blade', 'wire', 'strip'])\n", + " binned = load_mcstas(events).transpose(['strip', 'blade', 'wire', ])\n", "\n", " toa = (binned.bins.coords[\"event_time_zero\"] + binned.bins.coords[\"event_time_offset\"]).bins.concat().value\n", "\n", + " det_numbers_matching_order_of_binned = det_numbers.reshape(binned.shape)\n", + " det_numbers_matching_order_of_binned = det_numbers_matching_order_of_binned[::-1, ::-1, :]\n", " # IMPORTANT! we need to sort the arrays below according to toa,\n", " # so that the event_index does not get messed up!\n", " event_id = sc.sort(\n", " (\n", - " sc.bins_like(binned, sc.array(dims=binned.dims, values=det_numbers.reshape(binned.shape)))\n", + " sc.bins_like(binned, sc.array(dims=binned.dims, values=det_numbers_matching_order_of_binned))\n", " .bins.concat()\n", " .value\n", " ),\n", @@ -188,8 +190,15 @@ "mcstas_to_nexus(\n", " mcstas_data_file=estia_mcstas_example('reference'),\n", " template_nexus_file=\"/home/johannes/Downloads/estia_999999_00008786.hdf\",\n", - " outfile='out.nx',\n", - " nevents=20_000_000,\n", + " outfile='reference.nx',\n", + " nevents=10_000_000,\n", + ")\n", + "\n", + "mcstas_to_nexus(\n", + " mcstas_data_file=estia_mcstas_example('Ni/Ti-multilayer')[1],\n", + " template_nexus_file=\"/home/johannes/Downloads/estia_999999_00008786.hdf\",\n", + " outfile='niti_ml_1.nx',\n", + " nevents=10_000_000,\n", ")" ] }, @@ -199,6 +208,37 @@ "id": "4", "metadata": {}, "outputs": [], + "source": [ + "import scippneutron as scn\n", + "from ess.reflectometry.types import *\n", + "from ess.estia import EstiaWorkflow, EstiaMcStasWorkflow\n", + "\n", + "wfm = EstiaMcStasWorkflow()\n", + "wfm[Filename[SampleRun]] = estia_mcstas_example('Ni/Ti-multilayer')[1]\n", + "\n", + "scn.instrument_view(wfm.compute(RawDetector[SampleRun]).hist(), pixel_size=0.003, norm='log')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5", + "metadata": {}, + "outputs": [], + "source": [ + "import scippneutron as scn\n", + "wf = EstiaWorkflow()\n", + "wf[Filename[SampleRun]] = 'niti_ml_1.nx'\n", + "\n", + "scn.instrument_view(wf.compute(RawDetector[SampleRun]).hist(), pixel_size=3, norm='log')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6", + "metadata": {}, + "outputs": [], "source": [ "from ess.reflectometry.types import *\n", "from ess.estia import EstiaWorkflow, EstiaMcStasWorkflow\n", @@ -220,7 +260,7 @@ { "cell_type": "code", "execution_count": null, - "id": "5", + "id": "7", "metadata": {}, "outputs": [], "source": [ @@ -230,7 +270,7 @@ { "cell_type": "code", "execution_count": null, - "id": "6", + "id": "8", "metadata": {}, "outputs": [], "source": [ @@ -241,12 +281,90 @@ { "cell_type": "code", "execution_count": null, - "id": "7", + "id": "9", "metadata": {}, "outputs": [], "source": [ "wf.compute(RawDetector[SampleRun]).flatten(['blade', 'wire'], 'x').transpose().bins.sum().plot(norm='log')" ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "10", + "metadata": {}, + "outputs": [], + "source": [ + "import scipp as sc\n", + "\n", + "from ess.estia.data import estia_mcstas_example, estia_tof_lookup_table\n", + "from ess.estia import EstiaWorkflow\n", + "from ess.reflectometry.types import *\n", + "from ess.reflectometry import supermirror\n", + "\n", + "wf = EstiaWorkflow()\n", + "wf[Filename[SampleRun]] = 'niti_ml_1.nx'\n", + "wf[Filename[ReferenceRun]] = 'reference.nx'\n", + "# Select a region of interest:\n", + "wf[YIndexLimits] = sc.scalar(35), sc.scalar(64)\n", + "wf[ZIndexLimits] = sc.scalar(0), sc.scalar(48 * 32)\n", + "wf[BeamDivergenceLimits] = sc.scalar(-0.75, unit='deg'), sc.scalar(0.75, unit='deg')\n", + "\n", + "wf[supermirror.MValue] = sc.scalar(5, unit=sc.units.dimensionless)\n", + "wf[supermirror.CriticalEdge] = sc.scalar(float('inf'), unit='1/angstrom')\n", + "wf[supermirror.Alpha] = sc.scalar(0.25 / 0.088, unit=sc.units.angstrom)\n", + "\n", + "wf[DetectorSpatialResolution[RunType]] = 0.0025 * sc.units.m\n", + "\n", + "# Configure the binning of intermediate and final results:\n", + "wf[WavelengthBins] = sc.geomspace('wavelength', 3.5, 12, 2001, unit='angstrom')\n", + "wf[QBins] = 1000\n", + "\n", + "wf[TimeOfFlightLookupTableFilename] = estia_tof_lookup_table()\n", + "\n", + "wf[ProtonCurrent[SampleRun]] = sc.DataArray(\n", + " sc.array(dims=('time',), values=[]),\n", + " coords={'time': sc.array(dims=('time',), values=[], unit='s')})\n", + "wf[ProtonCurrent[ReferenceRun]] = sc.DataArray(\n", + " sc.array(dims=('time',), values=[]),\n", + " coords={'time': sc.array(dims=('time',), values=[], unit='s')})\n", + "\n", + "wf.visualize(ReflectivityOverQ)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "11", + "metadata": {}, + "outputs": [], + "source": [ + "da = wf.compute(ReducibleData[SampleRun])\n", + "da" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "12", + "metadata": {}, + "outputs": [], + "source": [ + "import scippneutron as scn\n", + "\n", + "scn.instrument_view(da.hist(), pixel_size=3, norm='log')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "13", + "metadata": {}, + "outputs": [], + "source": [ + "from ess.reduce.nexus.types import DetectorBankSizes\n", + "wf.compute(DetectorBankSizes)" + ] } ], "metadata": { From ed220dc7dfc7318cce47a99b5e1c286f99a2247c Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Thu, 4 Dec 2025 15:24:11 +0100 Subject: [PATCH 05/20] WIP --- tools/mcstas_to_nexus.ipynb | 56 +++++++++++++++++++++++++++++++++++-- 1 file changed, 54 insertions(+), 2 deletions(-) diff --git a/tools/mcstas_to_nexus.ipynb b/tools/mcstas_to_nexus.ipynb index 09881cd4..640493f8 100644 --- a/tools/mcstas_to_nexus.ipynb +++ b/tools/mcstas_to_nexus.ipynb @@ -108,6 +108,10 @@ "\n", " with h5py.File(outfile, \"r+\") as f:\n", "\n", + " # TODO: remove this when the files is fixed\n", + " # Fix error in Nexus file\n", + " f['/entry/instrument/source/transformations/location'][()] = -35.0512\n", + "\n", " detector_rotation = f[f\"{detector_entry_path}/transformations/detector_rotation\"]\n", " replace_dataset(\n", " detector_rotation,\n", @@ -227,6 +231,7 @@ "outputs": [], "source": [ "import scippneutron as scn\n", + "\n", "wf = EstiaWorkflow()\n", "wf[Filename[SampleRun]] = 'niti_ml_1.nx'\n", "\n", @@ -338,6 +343,53 @@ "id": "11", "metadata": {}, "outputs": [], + "source": [ + "da = wf.compute(ReducibleData[ReferenceRun])\n", + "da.masks.clear()\n", + "da.bins.masks.clear()\n", + "da.bin(divergence_angle=100).bins.size().plot()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "12", + "metadata": {}, + "outputs": [], + "source": [ + "da = wf.compute(ReducibleData[ReferenceRun])\n", + "da.masks.clear()\n", + "da.bins.masks.clear()\n", + "da.bin(theta=100).bins.size().plot()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "13", + "metadata": {}, + "outputs": [], + "source": [ + "da = wf.compute(ReflectivityOverQ)\n", + "da.bins.concat().value" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "14", + "metadata": {}, + "outputs": [], + "source": [ + "1.5 / 360 * 2 * 3.14" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "15", + "metadata": {}, + "outputs": [], "source": [ "da = wf.compute(ReducibleData[SampleRun])\n", "da" @@ -346,7 +398,7 @@ { "cell_type": "code", "execution_count": null, - "id": "12", + "id": "16", "metadata": {}, "outputs": [], "source": [ @@ -358,7 +410,7 @@ { "cell_type": "code", "execution_count": null, - "id": "13", + "id": "17", "metadata": {}, "outputs": [], "source": [ From 5738b836f6b0cf17436788d98d09896d0cdb4fe9 Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Wed, 10 Dec 2025 14:10:35 +0100 Subject: [PATCH 06/20] fix: sample number of events per weight instead of number of events --- src/ess/estia/mcstas.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/ess/estia/mcstas.py b/src/ess/estia/mcstas.py index 01bb9885..16a6cc4b 100644 --- a/src/ess/estia/mcstas.py +++ b/src/ess/estia/mcstas.py @@ -64,7 +64,7 @@ def parse_events_ascii(lines): raise ValueError('Could not parse the file as a list of events.') -def parse_events_h5(f, nevents_to_sample=None): +def parse_events_h5(f, events_to_sample_per_unit_weight=None): if isinstance(f, str): with h5py.File(f) as ff: return parse_events_h5(ff) @@ -76,7 +76,7 @@ def parse_events_h5(f, nevents_to_sample=None): 'NXentry/simulation/Param', ) events = events[()] - if nevents_to_sample is None: + if events_to_sample_per_unit_weight is None: da = sc.DataArray( # The squares on the variances is the correct way to load # weighted event data. @@ -91,10 +91,12 @@ def parse_events_h5(f, nevents_to_sample=None): da.coords[label] = sc.array(dims=['events'], values=events[:, i]) else: probabilities = events[:, 0] + p_total = probabilities.sum() + nevents_to_sample = round(events_to_sample_per_unit_weight * p_total) inds = np.random.choice( np.arange(len(probabilities)), nevents_to_sample, - p=probabilities / probabilities.sum(), + p=probabilities / p_total, ) da = sc.DataArray( sc.ones(dims=['events'], shape=(nevents_to_sample,), with_variances=True), From 57adaba7f92a65e51583be048c7aaa517e81376c Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Wed, 10 Dec 2025 14:11:15 +0100 Subject: [PATCH 07/20] fix: order matching estia Nexus file --- src/ess/estia/beamline.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ess/estia/beamline.py b/src/ess/estia/beamline.py index f4ec4656..e09b536f 100644 --- a/src/ess/estia/beamline.py +++ b/src/ess/estia/beamline.py @@ -6,9 +6,9 @@ DETECTOR_BANK_SIZES = { "multiblade_detector": { + "strip": 64, "blade": 48, "wire": 32, - "strip": 64, }, } From 0f3edfd29a5951b0d921c2b6ac6118dfffdca4b1 Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Wed, 10 Dec 2025 14:11:58 +0100 Subject: [PATCH 08/20] feat: add loaders for sample- and detector rotation --- src/ess/estia/conversions.py | 14 +++++++++----- src/ess/estia/load.py | 19 ++++++++++++++++++- src/ess/estia/workflow.py | 1 + 3 files changed, 28 insertions(+), 6 deletions(-) diff --git a/src/ess/estia/conversions.py b/src/ess/estia/conversions.py index 7746c29d..222b5818 100644 --- a/src/ess/estia/conversions.py +++ b/src/ess/estia/conversions.py @@ -10,8 +10,10 @@ from ..reflectometry.types import ( CoordTransformationGraph, DetectorLtotal, + DetectorRotation, RawDetector, RunType, + SampleRotation, ) @@ -77,6 +79,8 @@ def detector_ltotal_from_raw( 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, ) -> CoordTransformationGraph[RunType]: bank = detector_bank_sizes['multiblade_detector'] @@ -92,13 +96,13 @@ def coordinate_transformation_graph( position - sample_position.to(unit=position.unit) ), "Ltotal": lambda L1, L2: L1.to(unit=L2.unit) + L2, - 'sample_rotation': lambda: sc.scalar(1.0, unit='deg'), - 'detector_rotation': lambda: sc.scalar(3.65, unit='deg'), + 'sample_rotation': lambda: sample_rotation, + 'detector_rotation': lambda: detector_rotation, 'source_position': lambda: source_position, 'sample_position': lambda: sample_position, - 'blade': lambda: sc.arange('blade', 0, bank['blade']), - 'wire': lambda: sc.arange('wire', 0, bank['wire']), - 'strip': lambda: sc.arange('strip', 0, bank['strip']), + '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_size': lambda: sc.scalar(20.0, unit='mm'), } diff --git a/src/ess/estia/load.py b/src/ess/estia/load.py index feb3b7ad..77962947 100644 --- a/src/ess/estia/load.py +++ b/src/ess/estia/load.py @@ -4,10 +4,15 @@ import h5py import scipp as sc +from scippnexus import NXdetector, NXsample + +from ess.reduce.nexus.types import NeXusComponent from ..reflectometry.types import ( + DetectorRotation, Filename, RawDetector, + RawSampleRotation, RunType, SampleRotationOffset, ) @@ -113,4 +118,16 @@ def load_mcstas( return da -providers = () +def load_sample_rotation( + sample: NeXusComponent[NXsample, RunType], +) -> RawSampleRotation[RunType]: + return sample['sample_rotation']['value'] + + +def load_detector_rotation( + detector: NeXusComponent[NXdetector, RunType], +) -> DetectorRotation[RunType]: + return detector['transformations']['detector_rotation'].value[0].data + + +providers = (load_sample_rotation, load_detector_rotation) diff --git a/src/ess/estia/workflow.py b/src/ess/estia/workflow.py index 7b817768..ff8a6fc1 100644 --- a/src/ess/estia/workflow.py +++ b/src/ess/estia/workflow.py @@ -22,6 +22,7 @@ *maskings.providers, *normalization.providers, *orso.providers, + *load.providers, ) mcstas_providers = ( From 1b22f6018d9e7c3fcf4ce168636390dc5e0fb937 Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Wed, 10 Dec 2025 15:09:01 +0100 Subject: [PATCH 09/20] fix: add providers for geometry data to mcstas workflow --- src/ess/estia/conversions.py | 26 ++++++++++++++++++-------- src/ess/estia/load.py | 33 +++++++++++++++++++++++++++++++-- src/ess/estia/workflow.py | 3 +-- 3 files changed, 50 insertions(+), 12 deletions(-) diff --git a/src/ess/estia/conversions.py b/src/ess/estia/conversions.py index 222b5818..59991d8e 100644 --- a/src/ess/estia/conversions.py +++ b/src/ess/estia/conversions.py @@ -85,7 +85,6 @@ def coordinate_transformation_graph( ) -> CoordTransformationGraph[RunType]: bank = detector_bank_sizes['multiblade_detector'] return { - "wavelength": wavelength_from_tof, "theta": theta, "divergence_angle": divergence_angle, "Q": reflectometry_q, @@ -96,23 +95,34 @@ def coordinate_transformation_graph( position - sample_position.to(unit=position.unit) ), "Ltotal": lambda L1, L2: L1.to(unit=L2.unit) + L2, - 'sample_rotation': lambda: sample_rotation, - 'detector_rotation': lambda: detector_rotation, - 'source_position': lambda: source_position, - 'sample_position': lambda: sample_position, + '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_size': lambda: sc.scalar(20.0, unit='mm'), + "wavelength": wavelength_from_tof, + 'sample_rotation': lambda: sample_rotation, + 'detector_rotation': lambda: detector_rotation, + 'source_position': lambda: source_position, + 'sample_position': lambda: sample_position, } def mcstas_wavelength_coordinate_transformation_graph( - run_type: RunType, + source_position: Position[NXsource, RunType], + sample_position: Position[NXsample, RunType], + sample_rotation: SampleRotation[RunType], + detector_rotation: DetectorRotation[RunType], + detector_bank_sizes: DetectorBankSizes, ) -> CoordTransformationGraph[RunType]: return { - **coordinate_transformation_graph(), + **coordinate_transformation_graph( + source_position, + sample_position, + sample_rotation, + detector_rotation, + detector_bank_sizes, + ), "wavelength": lambda wavelength_from_mcstas: wavelength_from_mcstas, } diff --git a/src/ess/estia/load.py b/src/ess/estia/load.py index 77962947..e45ba46c 100644 --- a/src/ess/estia/load.py +++ b/src/ess/estia/load.py @@ -4,9 +4,9 @@ import h5py import scipp as sc -from scippnexus import NXdetector, NXsample +from scippnexus import NXdetector, NXsample, NXsource -from ess.reduce.nexus.types import NeXusComponent +from ess.reduce.nexus.types import NeXusComponent, Position from ..reflectometry.types import ( DetectorRotation, @@ -118,6 +118,28 @@ def load_mcstas( return da +def mcstas_load_sample_rotation(da: RawDetector[RunType]) -> RawSampleRotation[RunType]: + return da.coords['sample_rotation'] + + +def mcstas_load_detector_rotation( + da: RawDetector[RunType], +) -> DetectorRotation[RunType]: + return da.coords['detector_rotation'] + + +def mcstas_load_source_position( + da: RawDetector[RunType], +) -> Position[NXsource, RunType]: + return da.coords['source_position'] + + +def mcstas_load_sample_position( + da: RawDetector[RunType], +) -> Position[NXsample, RunType]: + return da.coords['sample_position'] + + def load_sample_rotation( sample: NeXusComponent[NXsample, RunType], ) -> RawSampleRotation[RunType]: @@ -131,3 +153,10 @@ def load_detector_rotation( providers = (load_sample_rotation, load_detector_rotation) +mcstas_providers = ( + load_mcstas_provider, + mcstas_load_sample_position, + mcstas_load_source_position, + mcstas_load_detector_rotation, + mcstas_load_sample_rotation, +) diff --git a/src/ess/estia/workflow.py b/src/ess/estia/workflow.py index ff8a6fc1..92c44a0b 100644 --- a/src/ess/estia/workflow.py +++ b/src/ess/estia/workflow.py @@ -27,8 +27,7 @@ mcstas_providers = ( *_general_providers, - *load.providers, - load.load_mcstas_provider, + *load.mcstas_providers, ) """List of providers for setting up a Sciline pipeline for McStas data. From 8927c6e0f97b119c852522076ddd33bbbcc6bc23 Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Wed, 10 Dec 2025 16:06:43 +0100 Subject: [PATCH 10/20] fix: add default sample rotation offset --- src/ess/estia/workflow.py | 1 + 1 file changed, 1 insertion(+) diff --git a/src/ess/estia/workflow.py b/src/ess/estia/workflow.py index 92c44a0b..eb317eea 100644 --- a/src/ess/estia/workflow.py +++ b/src/ess/estia/workflow.py @@ -60,6 +60,7 @@ def mcstas_default_parameters() -> dict: def default_parameters() -> dict: return { NeXusDetectorName: "multiblade_detector", + SampleRotationOffset[RunType]: sc.scalar(0.0, unit='deg'), } From 1e9bb6c8833b071c4fa46abc7b76c6087d3e5025 Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Wed, 10 Dec 2025 16:08:28 +0100 Subject: [PATCH 11/20] feat: mcstas to nexus notebook --- tools/mcstas_to_nexus.ipynb | 200 +++++++++--------------------------- 1 file changed, 48 insertions(+), 152 deletions(-) diff --git a/tools/mcstas_to_nexus.ipynb b/tools/mcstas_to_nexus.ipynb index 640493f8..e4a1261b 100644 --- a/tools/mcstas_to_nexus.ipynb +++ b/tools/mcstas_to_nexus.ipynb @@ -30,9 +30,7 @@ "import shutil\n", "\n", "from ess.estia.mcstas import parse_events_h5\n", - "from ess.estia.load import load_mcstas\n", - "from ess.estia.data import estia_mcstas_example\n", - "\n" + "from ess.estia.load import load_mcstas" ] }, { @@ -53,7 +51,7 @@ " mcstas_data_file: str,\n", " template_nexus_file: str,\n", " outfile: str,\n", - " nevents: int = 1_000_000,\n", + " nevents_per_weight: float = 1/3,\n", "):\n", " \"\"\"\n", " Store the events from a McStas Estia simulation in a NeXus CODA file.\n", @@ -76,13 +74,13 @@ " det_numbers = f[f\"{detector_entry_path}/detector_number\"][()]\n", "\n", " with h5py.File(mcstas_data_file) as f:\n", - " events = parse_events_h5(f, nevents_to_sample=nevents)\n", + " events = parse_events_h5(f, events_to_sample_per_unit_weight=nevents_per_weight)\n", " binned = load_mcstas(events).transpose(['strip', 'blade', 'wire', ])\n", "\n", " toa = (binned.bins.coords[\"event_time_zero\"] + binned.bins.coords[\"event_time_offset\"]).bins.concat().value\n", "\n", " det_numbers_matching_order_of_binned = det_numbers.reshape(binned.shape)\n", - " det_numbers_matching_order_of_binned = det_numbers_matching_order_of_binned[::-1, ::-1, :]\n", + " det_numbers_matching_order_of_binned = det_numbers_matching_order_of_binned[::-1, ::-1, ::-1]\n", " # IMPORTANT! we need to sort the arrays below according to toa,\n", " # so that the event_index does not get messed up!\n", " event_id = sc.sort(\n", @@ -184,26 +182,11 @@ ] }, { - "cell_type": "code", - "execution_count": null, + "cell_type": "markdown", "id": "3", "metadata": {}, - "outputs": [], "source": [ - "\n", - "mcstas_to_nexus(\n", - " mcstas_data_file=estia_mcstas_example('reference'),\n", - " template_nexus_file=\"/home/johannes/Downloads/estia_999999_00008786.hdf\",\n", - " outfile='reference.nx',\n", - " nevents=10_000_000,\n", - ")\n", - "\n", - "mcstas_to_nexus(\n", - " mcstas_data_file=estia_mcstas_example('Ni/Ti-multilayer')[1],\n", - " template_nexus_file=\"/home/johannes/Downloads/estia_999999_00008786.hdf\",\n", - " outfile='niti_ml_1.nx',\n", - " nevents=10_000_000,\n", - ")" + "### Generate nexus files from a Nexus template and McStas events" ] }, { @@ -213,29 +196,31 @@ "metadata": {}, "outputs": [], "source": [ - "import scippneutron as scn\n", - "from ess.reflectometry.types import *\n", - "from ess.estia import EstiaWorkflow, EstiaMcStasWorkflow\n", + "from ess.estia.data import estia_mcstas_example\n", "\n", - "wfm = EstiaMcStasWorkflow()\n", - "wfm[Filename[SampleRun]] = estia_mcstas_example('Ni/Ti-multilayer')[1]\n", "\n", - "scn.instrument_view(wfm.compute(RawDetector[SampleRun]).hist(), pixel_size=0.003, norm='log')" + "mcstas_to_nexus(\n", + " mcstas_data_file=estia_mcstas_example('reference'),\n", + " template_nexus_file=\"/home/johannes/Downloads/estia_999999_00008786.hdf\",\n", + " outfile='reference.nx',\n", + " nevents_per_weight=1,\n", + ")\n", + "\n", + "for i, example in enumerate(estia_mcstas_example('Ni/Ti-multilayer')):\n", + " mcstas_to_nexus(\n", + " mcstas_data_file=example,\n", + " template_nexus_file=\"/home/johannes/Downloads/estia_999999_00008786.hdf\",\n", + " outfile=f'niti_ml_{i}.nx',\n", + " nevents_per_weight=1,\n", + " )" ] }, { - "cell_type": "code", - "execution_count": null, + "cell_type": "markdown", "id": "5", "metadata": {}, - "outputs": [], "source": [ - "import scippneutron as scn\n", - "\n", - "wf = EstiaWorkflow()\n", - "wf[Filename[SampleRun]] = 'niti_ml_1.nx'\n", - "\n", - "scn.instrument_view(wf.compute(RawDetector[SampleRun]).hist(), pixel_size=3, norm='log')" + "## Detector view from the created nexus file" ] }, { @@ -245,31 +230,21 @@ "metadata": {}, "outputs": [], "source": [ + "from ess.estia import EstiaWorkflow\n", "from ess.reflectometry.types import *\n", - "from ess.estia import EstiaWorkflow, EstiaMcStasWorkflow\n", - "\n", - "n = 20_000_000\n", - "\n", - "wfm = EstiaMcStasWorkflow()\n", - "wfm[Filename[SampleRun]] = estia_mcstas_example('reference')\n", - "h = wfm.compute(RawDetector[SampleRun]).flatten(['blade', 'wire'], 'x').bins.sum()\n", - "p1 = sc.round(sc.values(0.9 * n * h / h.sum().value)).plot(norm='log')\n", "\n", "wf = EstiaWorkflow()\n", - "wf[Filename[SampleRun]] = 'out.nx'\n", - "p2 = wf.compute(RawDetector[SampleRun]).flatten(['blade', 'wire'], 'x').transpose().bins.sum().plot(norm='log')\n", + "wf[Filename[SampleRun]] = 'niti_ml_0.nx'\n", "\n", - "p1 + p2" + "scn.instrument_view(wf.compute(RawDetector[SampleRun]).hist(), pixel_size=3, norm='log')" ] }, { - "cell_type": "code", - "execution_count": null, + "cell_type": "markdown", "id": "7", "metadata": {}, - "outputs": [], "source": [ - "wfm.compute(RawDetector[SampleRun]).bins.concat().hist(event_time_offset=200).plot() + wf.compute(RawDetector[SampleRun]).bins.concat().hist(event_time_offset=200).plot()" + "## Compute reflectivity curves from the nexus files" ] }, { @@ -279,36 +254,13 @@ "metadata": {}, "outputs": [], "source": [ - "%matplotlib ipympl\n", - "sc.round(sc.values(0.9* n * h / h.sum().value)).plot(norm='log')" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "9", - "metadata": {}, - "outputs": [], - "source": [ - "wf.compute(RawDetector[SampleRun]).flatten(['blade', 'wire'], 'x').transpose().bins.sum().plot(norm='log')" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "10", - "metadata": {}, - "outputs": [], - "source": [ - "import scipp as sc\n", - "\n", "from ess.estia.data import estia_mcstas_example, estia_tof_lookup_table\n", "from ess.estia import EstiaWorkflow\n", "from ess.reflectometry.types import *\n", "from ess.reflectometry import supermirror\n", "\n", + "\n", "wf = EstiaWorkflow()\n", - "wf[Filename[SampleRun]] = 'niti_ml_1.nx'\n", "wf[Filename[ReferenceRun]] = 'reference.nx'\n", "# Select a region of interest:\n", "wf[YIndexLimits] = sc.scalar(35), sc.scalar(64)\n", @@ -323,10 +275,11 @@ "\n", "# Configure the binning of intermediate and final results:\n", "wf[WavelengthBins] = sc.geomspace('wavelength', 3.5, 12, 2001, unit='angstrom')\n", - "wf[QBins] = 1000\n", + "wf[QBins] = sc.geomspace('Q', 0.01, 0.5, 401, unit='1/angstrom')\n", "\n", "wf[TimeOfFlightLookupTableFilename] = estia_tof_lookup_table()\n", "\n", + "\n", "wf[ProtonCurrent[SampleRun]] = sc.DataArray(\n", " sc.array(dims=('time',), values=[]),\n", " coords={'time': sc.array(dims=('time',), values=[], unit='s')})\n", @@ -334,88 +287,31 @@ " sc.array(dims=('time',), values=[]),\n", " coords={'time': sc.array(dims=('time',), values=[], unit='s')})\n", "\n", - "wf.visualize(ReflectivityOverQ)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "11", - "metadata": {}, - "outputs": [], - "source": [ - "da = wf.compute(ReducibleData[ReferenceRun])\n", - "da.masks.clear()\n", - "da.bins.masks.clear()\n", - "da.bin(divergence_angle=100).bins.size().plot()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "12", - "metadata": {}, - "outputs": [], - "source": [ - "da = wf.compute(ReducibleData[ReferenceRun])\n", - "da.masks.clear()\n", - "da.bins.masks.clear()\n", - "da.bin(theta=100).bins.size().plot()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "13", - "metadata": {}, - "outputs": [], - "source": [ - "da = wf.compute(ReflectivityOverQ)\n", - "da.bins.concat().value" + "wf.visualize(ReflectivityOverQ, graph_attr={\"rankdir\": \"LR\"}, compact=True)" ] }, { "cell_type": "code", "execution_count": null, - "id": "14", - "metadata": {}, - "outputs": [], - "source": [ - "1.5 / 360 * 2 * 3.14" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "15", - "metadata": {}, - "outputs": [], - "source": [ - "da = wf.compute(ReducibleData[SampleRun])\n", - "da" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "16", - "metadata": {}, - "outputs": [], - "source": [ - "import scippneutron as scn\n", - "\n", - "scn.instrument_view(da.hist(), pixel_size=3, norm='log')" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "17", + "id": "9", "metadata": {}, "outputs": [], "source": [ - "from ess.reduce.nexus.types import DetectorBankSizes\n", - "wf.compute(DetectorBankSizes)" + "%matplotlib ipympl\n", + "from ess.reflectometry.tools import batch_compute\n", + "params = {\n", + " str(i): {\n", + " Filename[SampleRun]: f'niti_ml_{i}.nx'\n", + " }\n", + " for i in range(4)\n", + "}\n", + "\n", + "batch_compute(\n", + " wf,\n", + " params,\n", + " ReflectivityOverQ,\n", + " #scale_to_overlap=True\n", + ").hist().plot(norm='log', vmin=1e-7)" ] } ], From 396784b704f800cae8f19e343f0770e16c126ea3 Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Wed, 10 Dec 2025 16:37:30 +0100 Subject: [PATCH 12/20] fix: move geometry information from amor load function to coordtransform graph --- src/ess/amor/conversions.py | 63 +++++++++++++++++++++++++++- src/ess/amor/load.py | 19 --------- src/ess/amor/workflow.py | 3 +- src/ess/estia/conversions.py | 28 +++++++++++++ src/ess/estia/corrections.py | 2 +- src/ess/reflectometry/conversions.py | 28 ------------- 6 files changed, 91 insertions(+), 52 deletions(-) diff --git a/src/ess/amor/conversions.py b/src/ess/amor/conversions.py index a609c091..d531a763 100644 --- a/src/ess/amor/conversions.py +++ b/src/ess/amor/conversions.py @@ -6,14 +6,24 @@ from ..reflectometry.conversions import reflectometry_q from ..reflectometry.types import ( BeamDivergenceLimits, + BeamSize, CoordTransformationGraph, + DetectorRotation, RunType, + SampleRotation, + SampleSize, WavelengthBins, YIndexLimits, ZIndexLimits, ) from .geometry import Detector -from .types import GravityToggle +from .types import ( + ChopperDistance, + ChopperFrequency, + ChopperPhase, + ChopperSeparation, + GravityToggle, +) def theta(wavelength, pixel_divergence_angle, L2, sample_rotation, detector_rotation): @@ -137,7 +147,15 @@ def wavelength( def coordinate_transformation_graph( - gravity: GravityToggle, run_type: RunType + gravity: GravityToggle, + detector_rotation: DetectorRotation[RunType], + sample_rotation: SampleRotation[RunType], + chopper_phase: ChopperPhase[RunType], + chopper_frequency: ChopperFrequency[RunType], + chopper_distance: ChopperDistance[RunType], + chopper_separation: ChopperSeparation[RunType], + sample_size: SampleSize[RunType], + beam_size: BeamSize[RunType], ) -> CoordTransformationGraph[RunType]: return { "wavelength": wavelength, @@ -146,9 +164,50 @@ def coordinate_transformation_graph( "Q": reflectometry_q, "L1": lambda chopper_distance: sc.abs(chopper_distance), "L2": lambda distance_in_detector: distance_in_detector + Detector.distance, + "sample_rotation": lambda: sample_rotation.to(unit='rad'), + "detector_rotation": lambda: detector_rotation.to(unit='rad'), + "chopper_phase": lambda: chopper_phase, + "chopper_frequency": lambda: chopper_frequency, + "chopper_separation": lambda: chopper_separation, + "chopper_distance": lambda: chopper_distance, + "sample_size": lambda: sample_size, + "beam_size": lambda: beam_size, } +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", + "L1", + "L2", + "blade", + "wire", + "strip", + "z_index", + "sample_rotation", + "detector_rotation", + "sample_size", + "beam_size", + "chopper_phase", + "chopper_frequency", + "chopper_separation", + "chopper_distance", + ), + graph, + rename_dims=False, + keep_intermediate=False, + keep_aliases=False, + ) + + def _not_between(v, a, b): return (v < a) | (v > b) diff --git a/src/ess/amor/load.py b/src/ess/amor/load.py index 5d6efe59..3fcd8ffd 100644 --- a/src/ess/amor/load.py +++ b/src/ess/amor/load.py @@ -9,7 +9,6 @@ from ..reflectometry.load import load_nx from ..reflectometry.types import ( Beamline, - BeamSize, DetectorRotation, Filename, Measurement, @@ -20,9 +19,7 @@ RawDetector, RawSampleRotation, RunType, - SampleRotation, SampleRun, - SampleSize, ) from .geometry import pixel_coordinates_in_detector_system from .types import ( @@ -41,14 +38,6 @@ def load_detector( def load_events( detector: NeXusComponent[snx.NXdetector, RunType], - detector_rotation: DetectorRotation[RunType], - sample_rotation: SampleRotation[RunType], - chopper_phase: ChopperPhase[RunType], - chopper_frequency: ChopperFrequency[RunType], - chopper_distance: ChopperDistance[RunType], - chopper_separation: ChopperSeparation[RunType], - sample_size: SampleSize[RunType], - beam_size: BeamSize[RunType], ) -> RawDetector[RunType]: event_data = detector["data"] if 'event_time_zero' in event_data.coords: @@ -69,14 +58,6 @@ def load_events( "data" ].data.values - data.coords["sample_rotation"] = sample_rotation.to(unit='rad') - data.coords["detector_rotation"] = detector_rotation.to(unit='rad') - data.coords["chopper_phase"] = chopper_phase - data.coords["chopper_frequency"] = chopper_frequency - data.coords["chopper_separation"] = chopper_separation - data.coords["chopper_distance"] = chopper_distance - data.coords["sample_size"] = sample_size - data.coords["beam_size"] = beam_size return RawDetector[RunType](data) diff --git a/src/ess/amor/workflow.py b/src/ess/amor/workflow.py index 7c5fa606..576d7359 100644 --- a/src/ess/amor/workflow.py +++ b/src/ess/amor/workflow.py @@ -1,6 +1,5 @@ # Copyright (c) 2025 Scipp contributors (https://github.com/scipp) from ..reflectometry.conversions import ( - add_coords, add_proton_current_coord, add_proton_current_mask, ) @@ -16,7 +15,7 @@ YIndexLimits, ZIndexLimits, ) -from .conversions import add_masks +from .conversions import add_coords, add_masks def add_coords_masks_and_apply_corrections( diff --git a/src/ess/estia/conversions.py b/src/ess/estia/conversions.py index 59991d8e..48b99567 100644 --- a/src/ess/estia/conversions.py +++ b/src/ess/estia/conversions.py @@ -127,6 +127,34 @@ def mcstas_wavelength_coordinate_transformation_graph( } +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", + "L1", + "L2", + "blade", + "wire", + "strip", + "z_index", + "sample_rotation", + "detector_rotation", + "sample_size", + ), + graph, + rename_dims=False, + keep_intermediate=False, + keep_aliases=False, + ) + + providers = ( coordinate_transformation_graph, detector_ltotal_from_raw, diff --git a/src/ess/estia/corrections.py b/src/ess/estia/corrections.py index be30f19a..d440ab30 100644 --- a/src/ess/estia/corrections.py +++ b/src/ess/estia/corrections.py @@ -3,7 +3,6 @@ import scipp as sc from ..reflectometry.conversions import ( - add_coords, add_proton_current_coord, add_proton_current_mask, ) @@ -19,6 +18,7 @@ YIndexLimits, ZIndexLimits, ) +from .conversions import add_coords from .maskings import add_masks diff --git a/src/ess/reflectometry/conversions.py b/src/ess/reflectometry/conversions.py index f72e5bb7..d2866394 100644 --- a/src/ess/reflectometry/conversions.py +++ b/src/ess/reflectometry/conversions.py @@ -64,32 +64,4 @@ def add_proton_current_mask(da: sc.DataArray) -> sc.DataArray: return da -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", - "L1", - "L2", - "blade", - "wire", - "strip", - "z_index", - "sample_rotation", - "detector_rotation", - "sample_size", - ), - graph, - rename_dims=False, - keep_intermediate=False, - keep_aliases=False, - ) - - providers = () From 62c71fc91df569892a3b546a39eb6e75a992e75b Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Thu, 11 Dec 2025 09:08:42 +0100 Subject: [PATCH 13/20] docs: set reference detector to get parameters --- docs/user-guide/estia/simulated-spin-flip-sample.ipynb | 2 ++ 1 file changed, 2 insertions(+) diff --git a/docs/user-guide/estia/simulated-spin-flip-sample.ipynb b/docs/user-guide/estia/simulated-spin-flip-sample.ipynb index ba81692b..fb1e3f0b 100644 --- a/docs/user-guide/estia/simulated-spin-flip-sample.ipynb +++ b/docs/user-guide/estia/simulated-spin-flip-sample.ipynb @@ -353,6 +353,8 @@ "mask = sc.isnan(I0).sum(('blade', 'wire')) > sc.scalar(0, unit=None)\n", "wf[ReducedReference] = I0.assign_masks(nopolcal=mask.data).assign_coords(wavelength=wf.compute(WavelengthBins))\n", "\n", + "# Required to read sample rotation / similar parameters associated with the reference measurement\n", + "wf[RawDetector] = sc.io.load_hdf5(estia_mcstas_spin_flip_example('supermirror', 'offoff'))\n", "\n", "sample_name = 'spin_flip_sample'\n", "\n", From 98b768077dd7f428a59f0e008535391de1d57fb8 Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Thu, 11 Dec 2025 09:23:59 +0100 Subject: [PATCH 14/20] refactor: improve variable naming --- src/ess/estia/mcstas.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/ess/estia/mcstas.py b/src/ess/estia/mcstas.py index 16a6cc4b..df69f3a6 100644 --- a/src/ess/estia/mcstas.py +++ b/src/ess/estia/mcstas.py @@ -90,13 +90,13 @@ def parse_events_h5(f, events_to_sample_per_unit_weight=None): continue da.coords[label] = sc.array(dims=['events'], values=events[:, i]) else: - probabilities = events[:, 0] - p_total = probabilities.sum() - nevents_to_sample = round(events_to_sample_per_unit_weight * p_total) + weights = events[:, 0] + total_weight = weights.sum() + nevents_to_sample = round(events_to_sample_per_unit_weight * total_weight) inds = np.random.choice( - np.arange(len(probabilities)), + np.arange(len(weights)), nevents_to_sample, - p=probabilities / p_total, + p=weights / total_weight, ) da = sc.DataArray( sc.ones(dims=['events'], shape=(nevents_to_sample,), with_variances=True), From 0507320c64b211acb443be2de2ba0e99162045e3 Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Thu, 11 Dec 2025 09:27:58 +0100 Subject: [PATCH 15/20] fix: use coordtransformationgraph type from ess.reflectometry in ess.offspec --- src/ess/offspec/conversions.py | 2 +- src/ess/offspec/load.py | 10 ++++++++-- src/ess/offspec/workflow.py | 2 +- 3 files changed, 10 insertions(+), 4 deletions(-) diff --git a/src/ess/offspec/conversions.py b/src/ess/offspec/conversions.py index 5534e8cf..cc373068 100644 --- a/src/ess/offspec/conversions.py +++ b/src/ess/offspec/conversions.py @@ -3,10 +3,10 @@ from scippneutron.conversion.graph import beamline, tof from ..reflectometry.types import ( + CoordTransformationGraph, ReferenceRun, SampleRun, ) -from .types import CoordTransformationGraph def coordinate_transformation_graph_sample() -> CoordTransformationGraph[SampleRun]: diff --git a/src/ess/offspec/load.py b/src/ess/offspec/load.py index b98a0f75..c617a580 100644 --- a/src/ess/offspec/load.py +++ b/src/ess/offspec/load.py @@ -2,8 +2,14 @@ # Copyright (c) 2025 Scipp contributors (https://github.com/scipp) import scipp as sc -from ..reflectometry.types import Filename, RawDetector, ReferenceRun, RunType -from .types import CoordTransformationGraph, MonitorData, NeXusMonitorName +from ..reflectometry.types import ( + CoordTransformationGraph, + Filename, + RawDetector, + ReferenceRun, + RunType, +) +from .types import MonitorData, NeXusMonitorName def load_offspec_events( diff --git a/src/ess/offspec/workflow.py b/src/ess/offspec/workflow.py index 741ec6c9..499b083a 100644 --- a/src/ess/offspec/workflow.py +++ b/src/ess/offspec/workflow.py @@ -4,6 +4,7 @@ from ..reflectometry import providers as reflectometry_providers from ..reflectometry.types import ( + CoordTransformationGraph, RawDetector, ReducibleData, ReferenceRun, @@ -16,7 +17,6 @@ from .maskings import add_masks from .types import ( BackgroundMinWavelength, - CoordTransformationGraph, MonitorData, NeXusMonitorName, SpectrumLimits, From ef6ad61d41a880dc6477aa450c54db1d960d777d Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Thu, 11 Dec 2025 10:03:58 +0100 Subject: [PATCH 16/20] fix: make sample rotation group a NXlog --- src/ess/estia/load.py | 2 +- tools/mcstas_to_nexus.ipynb | 5 ++--- 2 files changed, 3 insertions(+), 4 deletions(-) diff --git a/src/ess/estia/load.py b/src/ess/estia/load.py index e45ba46c..37fc9ccb 100644 --- a/src/ess/estia/load.py +++ b/src/ess/estia/load.py @@ -143,7 +143,7 @@ def mcstas_load_sample_position( def load_sample_rotation( sample: NeXusComponent[NXsample, RunType], ) -> RawSampleRotation[RunType]: - return sample['sample_rotation']['value'] + return sample['sample_rotation'][0].data def load_detector_rotation( diff --git a/tools/mcstas_to_nexus.ipynb b/tools/mcstas_to_nexus.ipynb index e4a1261b..4bf1b044 100644 --- a/tools/mcstas_to_nexus.ipynb +++ b/tools/mcstas_to_nexus.ipynb @@ -133,10 +133,9 @@ " ]\n", " )\n", "\n", - " # TODO: Put sample_rotation below \"transformations\" group\n", - " # TODO: Set correct nexus classes on the sample_rotation group\n", " sample = f[\"entry/sample\"]\n", " sample_rotation = sample.create_group('sample_rotation')\n", + " sample_rotation.attrs['NX_class'] = 'NXlog'\n", " sample_rotation.create_dataset(\n", " 'value',\n", " data=[\n", @@ -232,6 +231,7 @@ "source": [ "from ess.estia import EstiaWorkflow\n", "from ess.reflectometry.types import *\n", + "import scippneutron as scn\n", "\n", "wf = EstiaWorkflow()\n", "wf[Filename[SampleRun]] = 'niti_ml_0.nx'\n", @@ -262,7 +262,6 @@ "\n", "wf = EstiaWorkflow()\n", "wf[Filename[ReferenceRun]] = 'reference.nx'\n", - "# Select a region of interest:\n", "wf[YIndexLimits] = sc.scalar(35), sc.scalar(64)\n", "wf[ZIndexLimits] = sc.scalar(0), sc.scalar(48 * 32)\n", "wf[BeamDivergenceLimits] = sc.scalar(-0.75, unit='deg'), sc.scalar(0.75, unit='deg')\n", From 046a8fb0f731edae1a70a1419fc18be84ea1d450 Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Fri, 12 Dec 2025 11:06:15 +0100 Subject: [PATCH 17/20] refactor: move mcstas related code to mcstas module --- src/ess/estia/conversions.py | 16 +--- src/ess/estia/load.py | 140 +------------------------------ src/ess/estia/mcstas.py | 157 ++++++++++++++++++++++++++++++++++- src/ess/estia/workflow.py | 13 ++- 4 files changed, 170 insertions(+), 156 deletions(-) diff --git a/src/ess/estia/conversions.py b/src/ess/estia/conversions.py index 48b99567..6084dc30 100644 --- a/src/ess/estia/conversions.py +++ b/src/ess/estia/conversions.py @@ -9,9 +9,7 @@ from ..reflectometry.conversions import reflectometry_q from ..reflectometry.types import ( CoordTransformationGraph, - DetectorLtotal, DetectorRotation, - RawDetector, RunType, SampleRotation, ) @@ -67,15 +65,6 @@ def divergence_angle( return sc.atan2(y=p.fields.x, x=p.fields.z) - detector_rotation.to(unit='rad') -def detector_ltotal_from_raw( - da: RawDetector[RunType], graph: CoordTransformationGraph[RunType] -) -> DetectorLtotal[RunType]: - return da.transform_coords( - ['Ltotal'], - graph=graph, - ).coords['Ltotal'] - - def coordinate_transformation_graph( source_position: Position[NXsource, RunType], sample_position: Position[NXsample, RunType], @@ -155,7 +144,4 @@ def add_coords( ) -providers = ( - coordinate_transformation_graph, - detector_ltotal_from_raw, -) +providers = (coordinate_transformation_graph,) diff --git a/src/ess/estia/load.py b/src/ess/estia/load.py index 37fc9ccb..e60f6d78 100644 --- a/src/ess/estia/load.py +++ b/src/ess/estia/load.py @@ -1,143 +1,14 @@ # SPDX-License-Identifier: BSD-3-Clause # Copyright (c) 2025 Scipp contributors (https://github.com/scipp) -import pathlib +from scippnexus import NXdetector, NXsample -import h5py -import scipp as sc -from scippnexus import NXdetector, NXsample, NXsource - -from ess.reduce.nexus.types import NeXusComponent, Position +from ess.reduce.nexus.types import NeXusComponent from ..reflectometry.types import ( DetectorRotation, - Filename, - RawDetector, RawSampleRotation, RunType, - SampleRotationOffset, ) -from .beamline import DETECTOR_BANK_SIZES -from .mcstas import parse_events_ascii, parse_events_h5 - - -def load_mcstas_provider( - filename: Filename[RunType], - sample_rotation_offset: SampleRotationOffset[RunType], -) -> RawDetector[RunType]: - """See :py:`load_mcstas`.""" - da = load_mcstas(filename) - da.coords['sample_rotation'] += sample_rotation_offset.to( - unit=da.coords['sample_rotation'].unit - ) - return RawDetector[RunType](da) - - -def load_mcstas( - filename: str | pathlib.Path | sc.DataArray, -) -> sc.DataArray: - """ - Load event data from a McStas run and reshape it - to look like what we would expect if - we loaded data from a real experiment file. - """ - if isinstance(filename, sc.DataArray): - da = filename - elif h5py.is_hdf5(filename): - with h5py.File(filename) as f: - da = parse_events_h5(f) - else: - with open(filename) as f: - da = parse_events_ascii(f) - - da.coords['sample_rotation'] = sc.scalar( - float(da.coords['omegaa'].value), unit='deg' - ) - da.coords['detector_rotation'] = 2 * da.coords['sample_rotation'] + sc.scalar( - 1.65, unit='deg' - ) - - nblades = DETECTOR_BANK_SIZES['multiblade_detector']['blade'] - nwires = DETECTOR_BANK_SIZES['multiblade_detector']['wire'] - nstrips = DETECTOR_BANK_SIZES['multiblade_detector']['strip'] - xbins = sc.linspace('x', -0.25, 0.25, nblades * nwires + 1) - ybins = sc.linspace('y', -0.13, 0.13, nstrips + 1) - da = da.bin(y=ybins, x=xbins).rename_dims({'y': 'strip'}) - da.coords['strip'] = sc.arange('strip', 0, nstrips) - da.coords['z_index'] = sc.arange('x', nblades * nwires - 1, -1, -1) - - # Information is not available in the mcstas output files, therefore it's hardcoded - da.coords['sample_position'] = sc.vector([0.264298, -0.427595, 35.0512], unit='m') - da.coords['source_position'] = sc.vector([0, 0, 0.0], unit='m') - da.coords['detector_position'] = sc.vector( - tuple(map(float, da.coords['position'].value.split(' '))), unit='m' - ) - - rotation_by_detector_rotation = sc.spatial.rotation( - value=[ - sc.scalar(0.0), - sc.sin(da.coords['detector_rotation'].to(unit='rad')), - sc.scalar(0.0), - sc.cos(da.coords['detector_rotation'].to(unit='rad')), - ] - ) - - position = sc.spatial.as_vectors( - x=sc.midpoints(da.coords['x']) * sc.scalar(1.0, unit='m'), - y=sc.midpoints(da.coords['y']) * sc.scalar(1.0, unit='m'), - z=sc.scalar(0.0, unit='m'), - ).transpose(da.dims) - da.coords['position'] = ( - da.coords['detector_position'] + rotation_by_detector_rotation * position - ) - - da.bins.coords['event_time_zero'] = ( - sc.scalar(0, unit='s') * da.bins.coords['t'] - ).to(unit='ns') - da.bins.coords['event_time_offset'] = ( - sc.scalar(1, unit='s') * da.bins.coords['t'] - ).to(unit='ns') % sc.scalar(1 / 14, unit='s').to(unit='ns') - da.bins.coords['wavelength_from_mcstas'] = ( - sc.scalar(1.0, unit='angstrom') * da.bins.coords['L'] - ) - - da.coords["sample_size"] = sc.scalar(1.0, unit='m') * float( - da.coords['sample_length'].value - ) - da.coords["beam_size"] = sc.scalar(2.0, unit='mm') - - da = da.fold( - 'x', - sizes={ - k: v - for k, v in DETECTOR_BANK_SIZES['multiblade_detector'].items() - if k in ('blade', 'wire') - }, - ) - da.bins.coords.pop('L') - da.bins.coords.pop('t') - return da - - -def mcstas_load_sample_rotation(da: RawDetector[RunType]) -> RawSampleRotation[RunType]: - return da.coords['sample_rotation'] - - -def mcstas_load_detector_rotation( - da: RawDetector[RunType], -) -> DetectorRotation[RunType]: - return da.coords['detector_rotation'] - - -def mcstas_load_source_position( - da: RawDetector[RunType], -) -> Position[NXsource, RunType]: - return da.coords['source_position'] - - -def mcstas_load_sample_position( - da: RawDetector[RunType], -) -> Position[NXsample, RunType]: - return da.coords['sample_position'] def load_sample_rotation( @@ -153,10 +24,3 @@ def load_detector_rotation( providers = (load_sample_rotation, load_detector_rotation) -mcstas_providers = ( - load_mcstas_provider, - mcstas_load_sample_position, - mcstas_load_source_position, - mcstas_load_detector_rotation, - mcstas_load_sample_rotation, -) diff --git a/src/ess/estia/mcstas.py b/src/ess/estia/mcstas.py index df69f3a6..edce4e78 100644 --- a/src/ess/estia/mcstas.py +++ b/src/ess/estia/mcstas.py @@ -1,10 +1,26 @@ # Copyright (c) 2025 Scipp contributors (https://github.com/scipp) +import pathlib + import h5py import numpy as np import pandas as pd import scipp as sc +from scippnexus import NXsample, NXsource + +from ess.reduce.nexus.types import Position -from ess.reflectometry.load import load_h5 +from ..reflectometry.load import load_h5 +from ..reflectometry.types import ( + CoordTransformationGraph, + DetectorLtotal, + DetectorRotation, + Filename, + RawDetector, + RawSampleRotation, + RunType, + SampleRotationOffset, +) +from .beamline import DETECTOR_BANK_SIZES def parse_metadata_ascii(lines): @@ -115,3 +131,142 @@ def parse_events_h5(f, events_to_sample_per_unit_weight=None): v = v.decode() da.coords[k] = sc.scalar(v) return da + + +def load_mcstas_provider( + filename: Filename[RunType], + sample_rotation_offset: SampleRotationOffset[RunType], +) -> RawDetector[RunType]: + """See :py:`load_mcstas`.""" + da = load_mcstas(filename) + da.coords['sample_rotation'] += sample_rotation_offset.to( + unit=da.coords['sample_rotation'].unit + ) + return RawDetector[RunType](da) + + +def load_mcstas( + filename: str | pathlib.Path | sc.DataArray, +) -> sc.DataArray: + """ + Load event data from a McStas run and reshape it + to look like what we would expect if + we loaded data from a real experiment file. + """ + if isinstance(filename, sc.DataArray): + da = filename + elif h5py.is_hdf5(filename): + with h5py.File(filename) as f: + da = parse_events_h5(f) + else: + with open(filename) as f: + da = parse_events_ascii(f) + + da.coords['sample_rotation'] = sc.scalar( + float(da.coords['omegaa'].value), unit='deg' + ) + da.coords['detector_rotation'] = 2 * da.coords['sample_rotation'] + sc.scalar( + 1.65, unit='deg' + ) + + nblades = DETECTOR_BANK_SIZES['multiblade_detector']['blade'] + nwires = DETECTOR_BANK_SIZES['multiblade_detector']['wire'] + nstrips = DETECTOR_BANK_SIZES['multiblade_detector']['strip'] + xbins = sc.linspace('x', -0.25, 0.25, nblades * nwires + 1) + ybins = sc.linspace('y', -0.13, 0.13, nstrips + 1) + da = da.bin(y=ybins, x=xbins).rename_dims({'y': 'strip'}) + da.coords['strip'] = sc.arange('strip', 0, nstrips) + da.coords['z_index'] = sc.arange('x', nblades * nwires - 1, -1, -1) + + # Information is not available in the mcstas output files, therefore it's hardcoded + da.coords['sample_position'] = sc.vector([0.264298, -0.427595, 35.0512], unit='m') + da.coords['source_position'] = sc.vector([0, 0, 0.0], unit='m') + da.coords['detector_position'] = sc.vector( + tuple(map(float, da.coords['position'].value.split(' '))), unit='m' + ) + + rotation_by_detector_rotation = sc.spatial.rotation( + value=[ + sc.scalar(0.0), + sc.sin(da.coords['detector_rotation'].to(unit='rad')), + sc.scalar(0.0), + sc.cos(da.coords['detector_rotation'].to(unit='rad')), + ] + ) + + position = sc.spatial.as_vectors( + x=sc.midpoints(da.coords['x']) * sc.scalar(1.0, unit='m'), + y=sc.midpoints(da.coords['y']) * sc.scalar(1.0, unit='m'), + z=sc.scalar(0.0, unit='m'), + ).transpose(da.dims) + da.coords['position'] = ( + da.coords['detector_position'] + rotation_by_detector_rotation * position + ) + + da.bins.coords['event_time_zero'] = ( + sc.scalar(0, unit='s') * da.bins.coords['t'] + ).to(unit='ns') + da.bins.coords['event_time_offset'] = ( + sc.scalar(1, unit='s') * da.bins.coords['t'] + ).to(unit='ns') % sc.scalar(1 / 14, unit='s').to(unit='ns') + da.bins.coords['wavelength_from_mcstas'] = ( + sc.scalar(1.0, unit='angstrom') * da.bins.coords['L'] + ) + + da.coords["sample_size"] = sc.scalar(1.0, unit='m') * float( + da.coords['sample_length'].value + ) + da.coords["beam_size"] = sc.scalar(2.0, unit='mm') + + da = da.fold( + 'x', + sizes={ + k: v + for k, v in DETECTOR_BANK_SIZES['multiblade_detector'].items() + if k in ('blade', 'wire') + }, + ) + da.bins.coords.pop('L') + da.bins.coords.pop('t') + return da + + +def load_sample_rotation(da: RawDetector[RunType]) -> RawSampleRotation[RunType]: + return da.coords['sample_rotation'] + + +def load_detector_rotation( + da: RawDetector[RunType], +) -> DetectorRotation[RunType]: + return da.coords['detector_rotation'] + + +def load_source_position( + da: RawDetector[RunType], +) -> Position[NXsource, RunType]: + return da.coords['source_position'] + + +def load_sample_position( + da: RawDetector[RunType], +) -> Position[NXsample, RunType]: + return da.coords['sample_position'] + + +def detector_ltotal_from_raw( + da: RawDetector[RunType], graph: CoordTransformationGraph[RunType] +) -> DetectorLtotal[RunType]: + return da.transform_coords( + ['Ltotal'], + graph=graph, + ).coords['Ltotal'] + + +providers = ( + load_mcstas_provider, + load_sample_position, + load_source_position, + load_detector_rotation, + load_sample_rotation, + detector_ltotal_from_raw, +) diff --git a/src/ess/estia/workflow.py b/src/ess/estia/workflow.py index eb317eea..afc79eb9 100644 --- a/src/ess/estia/workflow.py +++ b/src/ess/estia/workflow.py @@ -13,7 +13,16 @@ RunType, SampleRotationOffset, ) -from . import beamline, conversions, corrections, load, maskings, normalization, orso +from . import ( + beamline, + conversions, + corrections, + load, + maskings, + mcstas, + normalization, + orso, +) _general_providers = ( *reflectometry_providers, @@ -27,7 +36,7 @@ mcstas_providers = ( *_general_providers, - *load.mcstas_providers, + *mcstas.providers, ) """List of providers for setting up a Sciline pipeline for McStas data. From 5d6042aec952e578b0a7ede79bc7fa478e33942d Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Fri, 12 Dec 2025 11:10:38 +0100 Subject: [PATCH 18/20] refactor: move mcstas related into mcstas module --- src/ess/estia/conversions.py | 19 ------------------- src/ess/estia/mcstas.py | 23 ++++++++++++++++++++++- tests/estia/mcstas_data_test.py | 2 +- 3 files changed, 23 insertions(+), 21 deletions(-) diff --git a/src/ess/estia/conversions.py b/src/ess/estia/conversions.py index 6084dc30..d033cf7b 100644 --- a/src/ess/estia/conversions.py +++ b/src/ess/estia/conversions.py @@ -97,25 +97,6 @@ def coordinate_transformation_graph( } -def mcstas_wavelength_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, -) -> CoordTransformationGraph[RunType]: - return { - **coordinate_transformation_graph( - source_position, - sample_position, - sample_rotation, - detector_rotation, - detector_bank_sizes, - ), - "wavelength": lambda wavelength_from_mcstas: wavelength_from_mcstas, - } - - def add_coords( da: sc.DataArray, graph: dict, diff --git a/src/ess/estia/mcstas.py b/src/ess/estia/mcstas.py index edce4e78..ebb773c5 100644 --- a/src/ess/estia/mcstas.py +++ b/src/ess/estia/mcstas.py @@ -7,7 +7,7 @@ import scipp as sc from scippnexus import NXsample, NXsource -from ess.reduce.nexus.types import Position +from ess.reduce.nexus.types import DetectorBankSizes, Position from ..reflectometry.load import load_h5 from ..reflectometry.types import ( @@ -18,9 +18,11 @@ RawDetector, RawSampleRotation, RunType, + SampleRotation, SampleRotationOffset, ) from .beamline import DETECTOR_BANK_SIZES +from .conversions import coordinate_transformation_graph def parse_metadata_ascii(lines): @@ -262,6 +264,25 @@ def detector_ltotal_from_raw( ).coords['Ltotal'] +def mcstas_wavelength_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, +) -> CoordTransformationGraph[RunType]: + return { + **coordinate_transformation_graph( + source_position, + sample_position, + sample_rotation, + detector_rotation, + detector_bank_sizes, + ), + "wavelength": lambda wavelength_from_mcstas: wavelength_from_mcstas, + } + + providers = ( load_mcstas_provider, load_sample_position, diff --git a/tests/estia/mcstas_data_test.py b/tests/estia/mcstas_data_test.py index 46064e05..415c1db0 100644 --- a/tests/estia/mcstas_data_test.py +++ b/tests/estia/mcstas_data_test.py @@ -12,12 +12,12 @@ from orsopy import fileio from ess.estia import EstiaMcStasWorkflow -from ess.estia.conversions import mcstas_wavelength_coordinate_transformation_graph from ess.estia.data import ( estia_mcstas_reference_run, estia_mcstas_sample_run, estia_tof_lookup_table, ) +from ess.estia.mcstas import mcstas_wavelength_coordinate_transformation_graph from ess.reflectometry import orso from ess.reflectometry.types import ( BeamDivergenceLimits, From bbeb8c740c5576900f03bca28e993ba9da04a2d9 Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Fri, 12 Dec 2025 11:17:29 +0100 Subject: [PATCH 19/20] docs: remove instrument view and reduction from mcstas->nexus notebook --- tools/mcstas_to_nexus.ipynb | 105 ++---------------------------------- 1 file changed, 3 insertions(+), 102 deletions(-) diff --git a/tools/mcstas_to_nexus.ipynb b/tools/mcstas_to_nexus.ipynb index 4bf1b044..91c0ffc4 100644 --- a/tools/mcstas_to_nexus.ipynb +++ b/tools/mcstas_to_nexus.ipynb @@ -51,7 +51,7 @@ " mcstas_data_file: str,\n", " template_nexus_file: str,\n", " outfile: str,\n", - " nevents_per_weight: float = 1/3,\n", + " nevents_per_weight: float = 1,\n", "):\n", " \"\"\"\n", " Store the events from a McStas Estia simulation in a NeXus CODA file.\n", @@ -64,9 +64,9 @@ " NeXus file containing geometry and instrument info, used as a template.\n", " outfile:\n", " Output file to be written.\n", - " nevents:\n", + " nevents_per_weight:\n", " Number of events to have in the output file\n", - " (events are sampled from the probabilities of the mcstas events).\n", + " per unit weight in the McStas data.\n", " \"\"\"\n", " detector_entry_path = '/entry/instrument/multiblade_detector'\n", "\n", @@ -213,105 +213,6 @@ " nevents_per_weight=1,\n", " )" ] - }, - { - "cell_type": "markdown", - "id": "5", - "metadata": {}, - "source": [ - "## Detector view from the created nexus file" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "6", - "metadata": {}, - "outputs": [], - "source": [ - "from ess.estia import EstiaWorkflow\n", - "from ess.reflectometry.types import *\n", - "import scippneutron as scn\n", - "\n", - "wf = EstiaWorkflow()\n", - "wf[Filename[SampleRun]] = 'niti_ml_0.nx'\n", - "\n", - "scn.instrument_view(wf.compute(RawDetector[SampleRun]).hist(), pixel_size=3, norm='log')" - ] - }, - { - "cell_type": "markdown", - "id": "7", - "metadata": {}, - "source": [ - "## Compute reflectivity curves from the nexus files" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "8", - "metadata": {}, - "outputs": [], - "source": [ - "from ess.estia.data import estia_mcstas_example, estia_tof_lookup_table\n", - "from ess.estia import EstiaWorkflow\n", - "from ess.reflectometry.types import *\n", - "from ess.reflectometry import supermirror\n", - "\n", - "\n", - "wf = EstiaWorkflow()\n", - "wf[Filename[ReferenceRun]] = 'reference.nx'\n", - "wf[YIndexLimits] = sc.scalar(35), sc.scalar(64)\n", - "wf[ZIndexLimits] = sc.scalar(0), sc.scalar(48 * 32)\n", - "wf[BeamDivergenceLimits] = sc.scalar(-0.75, unit='deg'), sc.scalar(0.75, unit='deg')\n", - "\n", - "wf[supermirror.MValue] = sc.scalar(5, unit=sc.units.dimensionless)\n", - "wf[supermirror.CriticalEdge] = sc.scalar(float('inf'), unit='1/angstrom')\n", - "wf[supermirror.Alpha] = sc.scalar(0.25 / 0.088, unit=sc.units.angstrom)\n", - "\n", - "wf[DetectorSpatialResolution[RunType]] = 0.0025 * sc.units.m\n", - "\n", - "# Configure the binning of intermediate and final results:\n", - "wf[WavelengthBins] = sc.geomspace('wavelength', 3.5, 12, 2001, unit='angstrom')\n", - "wf[QBins] = sc.geomspace('Q', 0.01, 0.5, 401, unit='1/angstrom')\n", - "\n", - "wf[TimeOfFlightLookupTableFilename] = estia_tof_lookup_table()\n", - "\n", - "\n", - "wf[ProtonCurrent[SampleRun]] = sc.DataArray(\n", - " sc.array(dims=('time',), values=[]),\n", - " coords={'time': sc.array(dims=('time',), values=[], unit='s')})\n", - "wf[ProtonCurrent[ReferenceRun]] = sc.DataArray(\n", - " sc.array(dims=('time',), values=[]),\n", - " coords={'time': sc.array(dims=('time',), values=[], unit='s')})\n", - "\n", - "wf.visualize(ReflectivityOverQ, graph_attr={\"rankdir\": \"LR\"}, compact=True)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "9", - "metadata": {}, - "outputs": [], - "source": [ - "%matplotlib ipympl\n", - "from ess.reflectometry.tools import batch_compute\n", - "params = {\n", - " str(i): {\n", - " Filename[SampleRun]: f'niti_ml_{i}.nx'\n", - " }\n", - " for i in range(4)\n", - "}\n", - "\n", - "batch_compute(\n", - " wf,\n", - " params,\n", - " ReflectivityOverQ,\n", - " #scale_to_overlap=True\n", - ").hist().plot(norm='log', vmin=1e-7)" - ] } ], "metadata": { From ec449edc3e92dfd1cd039dbc46af7d7d15005a69 Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Fri, 12 Dec 2025 11:19:36 +0100 Subject: [PATCH 20/20] fix: import --- docs/user-guide/estia/estia-advanced-mcstas-reduction.ipynb | 4 ++-- docs/user-guide/estia/simulated-spin-flip-sample.ipynb | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/docs/user-guide/estia/estia-advanced-mcstas-reduction.ipynb b/docs/user-guide/estia/estia-advanced-mcstas-reduction.ipynb index abe221af..325eb15f 100644 --- a/docs/user-guide/estia/estia-advanced-mcstas-reduction.ipynb +++ b/docs/user-guide/estia/estia-advanced-mcstas-reduction.ipynb @@ -76,7 +76,7 @@ "\n", "# If you want to use the ground truth wavelengths from the McStas simulation\n", "# instead of the wavelengths from the frame unwrapping, uncomment the lines below:\n", - "#from ess.estia.conversions import mcstas_wavelength_coordinate_transformation_graph\n", + "#from ess.estia.mcstas import mcstas_wavelength_coordinate_transformation_graph\n", "#wf.insert( mcstas_wavelength_coordinate_transformation_graph )\n", "\n", "wf[TimeOfFlightLookupTableFilename] = estia_tof_lookup_table()\n", @@ -448,7 +448,7 @@ "\n", "# If you want to use the ground truth wavelengths from the McStas simulation\n", "# instead of the wavelengths from the frame unwrapping, uncomment the lines below:\n", - "#from ess.estia.conversions import mcstas_wavelength_coordinate_transformation_graph\n", + "#from ess.estia.mcstas import mcstas_wavelength_coordinate_transformation_graph\n", "#wf.insert( mcstas_wavelength_coordinate_transformation_graph )\n", "\n", "wf[TimeOfFlightLookupTableFilename] = estia_tof_lookup_table()\n", diff --git a/docs/user-guide/estia/simulated-spin-flip-sample.ipynb b/docs/user-guide/estia/simulated-spin-flip-sample.ipynb index fb1e3f0b..b9ab6d82 100644 --- a/docs/user-guide/estia/simulated-spin-flip-sample.ipynb +++ b/docs/user-guide/estia/simulated-spin-flip-sample.ipynb @@ -53,7 +53,7 @@ " estia_tof_lookup_table,\n", ")\n", "from ess.reflectometry.figures import wavelength_z_figure\n", - "from ess.estia.conversions import mcstas_wavelength_coordinate_transformation_graph\n", + "from ess.estia.mcstas import mcstas_wavelength_coordinate_transformation_graph\n", "\n", "from ess.polarization import ( # noqa: F401\n", " Up, Down, ReducedSampleDataBySpinChannel, SupermirrorEfficiencyFunction, Polarizer, Analyzer, PolarizationAnalysisWorkflow,\n",