[ESSDIFFRACTION] docs: add user guide for the powder reduction workflow

packages/essdiffraction/docs/user-guide/beer/beer_powder_mcstas_analytical.ipynb

@@ -0,0 +1,345 @@
+{
+ "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. To run the notebook you will have to install `pooch` to access them: `pip install pooch`."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "2",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%matplotlib ipympl\n",
+    "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. 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": [
+    "To keep the events in the output, to allow post reduction filtering or other processing, set `KeepEvents` to `True`."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "8",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "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"
+   ]
+  },
+  {
+   "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
+}

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
 ```