From bf13b6ef6b314c3b8b4065283a4a6b0276cc3d42 Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Tue, 23 Jun 2026 16:19:53 +0200 Subject: [PATCH 1/2] docs: add user guide for the powder reduction workflow --- .../beer/beer_powder_mcstas_analytical.ipynb | 344 ++++++++++++++++++ .../docs/user-guide/beer/index.md | 15 + packages/essdiffraction/src/ess/beer/data.py | 27 +- 3 files changed, 385 insertions(+), 1 deletion(-) create mode 100644 packages/essdiffraction/docs/user-guide/beer/beer_powder_mcstas_analytical.ipynb diff --git a/packages/essdiffraction/docs/user-guide/beer/beer_powder_mcstas_analytical.ipynb b/packages/essdiffraction/docs/user-guide/beer/beer_powder_mcstas_analytical.ipynb new file mode 100644 index 000000000..8d6973fb6 --- /dev/null +++ b/packages/essdiffraction/docs/user-guide/beer/beer_powder_mcstas_analytical.ipynb @@ -0,0 +1,344 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "0", + "metadata": {}, + "source": [ + "# BEER McStas powder data reduction\n", + "\n", + "- Audience: BEER users reducing McStas pulse-shaping powder data\n", + "- Prerequisites: Basic knowledge of Scipp and Sciline workflows\n", + "\n", + "This guide reduces a McStas silicon-in-vanadium-can sample run using matching vanadium and empty-can runs. The result is `EmptyCanSubtractedIofDspacing`, the normalized intensity as a function of $d$-spacing after subtracting the empty can." + ] + }, + { + "cell_type": "markdown", + "id": "1", + "metadata": {}, + "source": [ + "The example files are downloaded through the ESSdiffraction data registry. Accessing them requires the `pooch` package." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2", + "metadata": {}, + "outputs": [], + "source": [ + "import scipp as sc\n", + "\n", + "from ess import beer, powder\n", + "import ess.beer.data # noqa: F401\n", + "from ess.beer.types import DetectorBank\n", + "from ess.diffraction.peaks import dspacing_peaks_from_cif\n", + "from ess.powder.types import *" + ] + }, + { + "cell_type": "markdown", + "id": "3", + "metadata": {}, + "source": [ + "## Configure the workflow\n", + "\n", + "Start with the BEER powder McStas workflow and set the run files. We use monitor-histogram normalization here because the McStas files include a wavelength monitor." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4", + "metadata": {}, + "outputs": [], + "source": [ + "workflow = beer.BeerPowderMcStasWorkflowAnalytical(\n", + " run_norm=powder.RunNormalization.monitor_histogram\n", + ")\n", + "\n", + "workflow[Filename[SampleRun]] = beer.data.mcstas_powder_silicon_in_vanadium_can()\n", + "workflow[Filename[VanadiumRun]] = beer.data.mcstas_powder_vanadium()\n", + "workflow[Filename[EmptyCanRun]] = beer.data.mcstas_powder_empty_can()" + ] + }, + { + "cell_type": "markdown", + "id": "5", + "metadata": {}, + "source": [ + "Next set the detector bank and the $d$-spacing range. This is the first, coarse region-of-interest choice: the notebook only reduces the north bank and only histograms the range from 0.6 Å to 2.0 Å." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6", + "metadata": {}, + "outputs": [], + "source": [ + "workflow[DetectorBank] = DetectorBank.north\n", + "workflow[DspacingBins] = sc.linspace(\"dspacing\", 0.6, 2.0, 701, unit=\"angstrom\")" + ] + }, + { + "cell_type": "markdown", + "id": "7", + "metadata": {}, + "source": [ + "The subtraction is done before the final histogram, so keep the sample, vanadium, and empty-can data as events through the correction and subtraction steps." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8", + "metadata": {}, + "outputs": [], + "source": [ + "workflow[KeepEvents[SampleRun]] = KeepEvents[SampleRun](True)\n", + "workflow[KeepEvents[VanadiumRun]] = KeepEvents[VanadiumRun](True)\n", + "workflow[KeepEvents[EmptyCanRun]] = KeepEvents[EmptyCanRun](True)\n", + "\n", + "workflow[MaskedDetectorIDs] = MaskedDetectorIDs({})\n", + "workflow[UncertaintyBroadcastMode] = UncertaintyBroadcastMode.drop" + ] + }, + { + "cell_type": "markdown", + "id": "9", + "metadata": {}, + "source": [ + "A finer region of interest can be applied with the standard powder masks. Mask callables return `True` for values to remove. This example keeps a broad two-theta range for the selected detector data and leaves the other masks disabled." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "10", + "metadata": {}, + "outputs": [], + "source": [ + "two_theta_min = sc.scalar(76.0, unit=\"deg\").to(unit=\"rad\")\n", + "two_theta_max = sc.scalar(104.0, unit=\"deg\").to(unit=\"rad\")\n", + "\n", + "workflow[TwoThetaMask] = lambda two_theta: (two_theta < two_theta_min) | (\n", + " two_theta > two_theta_max\n", + ")\n", + "workflow[TofMask] = None\n", + "workflow[WavelengthMask] = None" + ] + }, + { + "cell_type": "markdown", + "id": "11", + "metadata": {}, + "source": [ + "Some events can have ambiguous wavelength assignments. We first keep all detector events so the baseline spectrum shows the full result." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "12", + "metadata": {}, + "outputs": [], + "source": [ + "workflow[LookupTableRelativeErrorThreshold] = {\n", + " \"detector\": float(\"inf\"),\n", + " \"monitor_bunker\": float(\"inf\"),\n", + " \"monitor_cave\": float(\"inf\"),\n", + "}" + ] + }, + { + "cell_type": "markdown", + "id": "13", + "metadata": {}, + "source": [ + "## Compute the empty-can-subtracted spectrum" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "14", + "metadata": {}, + "outputs": [], + "source": [ + "iof_dspacing = workflow.compute(EmptyCanSubtractedIofDspacing)\n", + "histogram = iof_dspacing.hist()\n", + "histogram" + ] + }, + { + "cell_type": "markdown", + "id": "15", + "metadata": {}, + "source": [ + "Compute the expected silicon peak positions from the CIF and draw the peaks that fall in the plotted $d$ range." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "16", + "metadata": {}, + "outputs": [], + "source": [ + "silicon_peaks = dspacing_peaks_from_cif(\"codid::1526655\").coords[\"dspacing\"]\n", + "silicon_peaks = silicon_peaks[\n", + " (silicon_peaks > histogram.coords[\"dspacing\"].min())\n", + " & (silicon_peaks < histogram.coords[\"dspacing\"].max())\n", + "]\n", + "\n", + "fig = histogram.plot(\n", + " title=\"Intensity over dspacing\",\n", + " xlabel=r\"$d$-spacing [$\\AA$]\",\n", + " ylabel=f\"Intensity [{histogram.unit}]\",\n", + " xmin=sc.scalar(0.6, unit=\"angstrom\"),\n", + " xmax=sc.scalar(2.0, unit=\"angstrom\"),\n", + ")\n", + "\n", + "ymin, ymax = fig.ax.get_ylim()\n", + "for peak in silicon_peaks.values:\n", + " fig.ax.axvline(peak, color=\"black\", alpha=0.25, linewidth=1)\n", + "\n", + "fig.ax.set_ylim(ymin, ymax)\n", + "fig" + ] + }, + { + "cell_type": "markdown", + "id": "17", + "metadata": {}, + "source": [ + "## Reduce frame-overlap artifacts\n", + "\n", + "The detector wavelength-assignment relative-error threshold removes events in regions where the wavelength assignment is ambiguous. The thresholds are set separately for the detector and the monitors." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "18", + "metadata": {}, + "outputs": [], + "source": [ + "thresholded_workflow = workflow.copy()\n", + "thresholded_workflow[LookupTableRelativeErrorThreshold] = {\n", + " \"detector\": 0.005,\n", + " \"monitor_bunker\": float(\"inf\"),\n", + " \"monitor_cave\": float(\"inf\"),\n", + "}\n", + "\n", + "thresholded_histogram = thresholded_workflow.compute(\n", + " EmptyCanSubtractedIofDspacing\n", + ").hist()" + ] + }, + { + "cell_type": "markdown", + "id": "19", + "metadata": {}, + "source": [ + "The threshold removes events in specific regions of scattering angle and wavelength. We can see this directly by comparing `CorrectedDetector` before and after applying the threshold." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "20", + "metadata": {}, + "outputs": [], + "source": [ + "unfiltered_detector = workflow.compute(CorrectedDetector[SampleRun])\n", + "filtered_detector = thresholded_workflow.compute(CorrectedDetector[SampleRun])\n", + "\n", + "two_theta_bins = sc.linspace(\"two_theta\", 76.0, 104.0, 201, unit=\"deg\").to(\n", + " unit=\"rad\"\n", + ")\n", + "wavelength_bins = sc.linspace(\"wavelength\", 0.8, 3.2, 241, unit=\"angstrom\")\n", + "\n", + "unfiltered_map = unfiltered_detector.hist(\n", + " two_theta=two_theta_bins,\n", + " wavelength=wavelength_bins,\n", + " dim=unfiltered_detector.dims,\n", + ")\n", + "filtered_map = filtered_detector.hist(\n", + " two_theta=two_theta_bins,\n", + " wavelength=wavelength_bins,\n", + " dim=filtered_detector.dims,\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "21", + "metadata": {}, + "outputs": [], + "source": [ + "filtered_map.plot(title=\"Kept by threshold\", norm=\"log\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "22", + "metadata": {}, + "outputs": [], + "source": [ + "(unfiltered_map - filtered_map).plot(title=\"Removed by threshold\", norm=\"log\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "23", + "metadata": {}, + "outputs": [], + "source": [ + "fig = sc.plot(\n", + " {\n", + " \"baseline\": histogram,\n", + " \"detector relative error < 0.02\": thresholded_histogram,\n", + " },\n", + " xlabel=r\"$d$-spacing [$\\AA$]\",\n", + " ylabel=f\"Intensity [{histogram.unit}]\",\n", + " xmin=sc.scalar(0.6, unit=\"angstrom\"),\n", + " xmax=sc.scalar(2.0, unit=\"angstrom\"),\n", + ")\n", + "\n", + "ymin, ymax = fig.ax.get_ylim()\n", + "for peak in silicon_peaks.values:\n", + " fig.ax.axvline(peak, color=\"black\", alpha=0.25, linewidth=1)\n", + "\n", + "fig.ax.set_ylim(ymin, ymax)\n", + "fig" + ] + } + ], + "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 +} diff --git a/packages/essdiffraction/docs/user-guide/beer/index.md b/packages/essdiffraction/docs/user-guide/beer/index.md index dd6242b3a..8465fcc7d 100644 --- a/packages/essdiffraction/docs/user-guide/beer/index.md +++ b/packages/essdiffraction/docs/user-guide/beer/index.md @@ -18,6 +18,20 @@ ``` ::: +:::{grid-item-card} Powder McStas data reduction +:link: beer_powder_mcstas_analytical.ipynb +:text-align: center + +```{image} ../../_static/thumbnails/beer_mcstas_light.svg +:class: only-light +:width: 75% +``` +```{image} ../../_static/thumbnails/beer_mcstas_dark.svg +:class: only-dark +:width: 75% +``` +::: + :::: ```{toctree} @@ -26,4 +40,5 @@ hidden: --- beer_modulation_mcstas +beer_powder_mcstas_analytical ``` diff --git a/packages/essdiffraction/src/ess/beer/data.py b/packages/essdiffraction/src/ess/beer/data.py index b8938c9e9..64f09f26a 100644 --- a/packages/essdiffraction/src/ess/beer/data.py +++ b/packages/essdiffraction/src/ess/beer/data.py @@ -8,7 +8,13 @@ from ess.reduce.data import make_registry -__all__ = ["mcstas_duplex", "mcstas_silicon_medium_resolution"] +__all__ = [ + "mcstas_duplex", + "mcstas_powder_empty_can", + "mcstas_powder_silicon_in_vanadium_can", + "mcstas_powder_vanadium", + "mcstas_silicon_medium_resolution", +] _registry = make_registry( "ess/beer", @@ -46,6 +52,10 @@ # - only used to verify we can load the 3D geometry. "few_neutrons_3d_detector_example.h5": "md5:88cbe29cb539c8acebf9fd7cee9d3c57", "silicon-mode9-3d-more-neutrons.h5": "md5:a6afdc4ee9827a57f88c2a3c6ef27383", + # Pulse-shaping powder diffraction simulation tutorial data. + "beer-powder-si-in-vcan-mode4.h5": "md5:7c1ba0fb0b2bb812ff59d0e4e1021ef7", + "beer-powder-vanadium-mode4.h5": "md5:17eed50ddf979789259b212d0be15f37", + "beer-powder-empty-can-mode4.h5": "md5:c374d5a8379749fb807b8f835291ce18", }, ) @@ -79,6 +89,21 @@ def mcstas_silicon_new_model(mode: int) -> Path: return _registry.get_path(f'silicon-mode{mode}-new-model.h5') +def mcstas_powder_silicon_in_vanadium_can() -> Path: + """Return the BEER McStas silicon-in-vanadium-can powder sample file.""" + return _registry.get_path('beer-powder-si-in-vcan-mode4.h5') + + +def mcstas_powder_vanadium() -> Path: + """Return the BEER McStas powder vanadium normalization file.""" + return _registry.get_path('beer-powder-vanadium-mode4.h5') + + +def mcstas_powder_empty_can() -> Path: + """Return the BEER McStas powder empty-can file.""" + return _registry.get_path('beer-powder-empty-can-mode4.h5') + + def mcstas_few_neutrons_3d_detector_example(): return _registry.get_path('few_neutrons_3d_detector_example.h5') From 2464ec78db6ea5c75a82da12a24b727d180a3395 Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Wed, 24 Jun 2026 14:47:52 +0200 Subject: [PATCH 2/2] docs: fix wordings --- .../beer/beer_powder_mcstas_analytical.ipynb | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/packages/essdiffraction/docs/user-guide/beer/beer_powder_mcstas_analytical.ipynb b/packages/essdiffraction/docs/user-guide/beer/beer_powder_mcstas_analytical.ipynb index 8d6973fb6..2b80b7eb6 100644 --- a/packages/essdiffraction/docs/user-guide/beer/beer_powder_mcstas_analytical.ipynb +++ b/packages/essdiffraction/docs/user-guide/beer/beer_powder_mcstas_analytical.ipynb @@ -18,7 +18,7 @@ "id": "1", "metadata": {}, "source": [ - "The example files are downloaded through the ESSdiffraction data registry. Accessing them requires the `pooch` package." + "The example files are downloaded through the ESSdiffraction data registry. To run the notebook you will have to install `pooch` to access them: `pip install pooch`." ] }, { @@ -28,6 +28,7 @@ "metadata": {}, "outputs": [], "source": [ + "%matplotlib ipympl\n", "import scipp as sc\n", "\n", "from ess import beer, powder\n", @@ -68,7 +69,7 @@ "id": "5", "metadata": {}, "source": [ - "Next set the detector bank and the $d$-spacing range. This is the first, coarse region-of-interest choice: the notebook only reduces the north bank and only histograms the range from 0.6 Å to 2.0 Å." + "Next set the detector bank and the $d$-spacing range. The notebook only reduces the north bank and only histograms the range from 0.6 Å to 2.0 Å." ] }, { @@ -87,7 +88,7 @@ "id": "7", "metadata": {}, "source": [ - "The subtraction is done before the final histogram, so keep the sample, vanadium, and empty-can data as events through the correction and subtraction steps." + "To keep the events in the output, to allow post reduction filtering or other processing, set `KeepEvents` to `True`." ] }, { @@ -97,9 +98,9 @@ "metadata": {}, "outputs": [], "source": [ - "workflow[KeepEvents[SampleRun]] = KeepEvents[SampleRun](True)\n", - "workflow[KeepEvents[VanadiumRun]] = KeepEvents[VanadiumRun](True)\n", - "workflow[KeepEvents[EmptyCanRun]] = KeepEvents[EmptyCanRun](True)\n", + "workflow[KeepEvents[SampleRun]] = KeepEvents(True)\n", + "workflow[KeepEvents[VanadiumRun]] = KeepEvents(True)\n", + "workflow[KeepEvents[EmptyCanRun]] = KeepEvents(True)\n", "\n", "workflow[MaskedDetectorIDs] = MaskedDetectorIDs({})\n", "workflow[UncertaintyBroadcastMode] = UncertaintyBroadcastMode.drop"