Merge branch 'require-chopper-name' of github.com:scipp/tof into requ…

…ire-chopper-name

docs/ess/index.rst

@@ -5,3 +5,4 @@ ESS instruments
    :maxdepth: 2
 
    dream
+   odin

docs/ess/odin.ipynb

@@ -0,0 +1,385 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "0",
+   "metadata": {},
+   "source": [
+    "# ODIN in WFM mode\n",
+    "\n",
+    "This is a simulation of the ODIN chopper cascade in WFM mode.\n",
+    "We also show how one can convert the neutron arrival times at the detector to wavelength."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "1",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "import scipp as sc\n",
+    "import plopp as pp\n",
+    "import tof\n",
+    "\n",
+    "Hz = sc.Unit(\"Hz\")\n",
+    "deg = sc.Unit(\"deg\")\n",
+    "meter = sc.Unit(\"m\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "2",
+   "metadata": {},
+   "source": [
+    "## Create a source pulse\n",
+    "\n",
+    "We first create a source with 4 pulses containing 800,000 neutrons each,\n",
+    "and whose distribution follows the ESS time and wavelength profiles (both thermal and cold neutrons are included)."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "3",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "source = tof.Source(facility=\"ess\", neutrons=500_000, pulses=4)\n",
+    "source.plot()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "4",
+   "metadata": {},
+   "source": [
+    "## Component set-up\n",
+    "\n",
+    "### Choppers\n",
+    "\n",
+    "The ODIN chopper cascade consists of:\n",
+    "\n",
+    "- 2 WFM choppers\n",
+    "- 5 frame-overlap choppers\n",
+    "- 2 band-control choppers\n",
+    "- 1 T0 chopper"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "5",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "parameters = {\n",
+    "    \"WFMC_1\": {\n",
+    "        \"frequency\": 56.0,\n",
+    "        \"phase\": 93.244,\n",
+    "        \"distance\": 6.85,\n",
+    "        \"open\": [-1.9419, 49.5756, 98.9315, 146.2165, 191.5176, 234.9179],\n",
+    "        \"close\": [1.9419, 55.7157, 107.2332, 156.5891, 203.8741, 249.1752]\n",
+    "    },\n",
+    "    \"WFMC_2\": {\n",
+    "        \"frequency\": 56.0,\n",
+    "        \"phase\": 97.128,\n",
+    "        \"distance\": 7.15,\n",
+    "        \"open\": [-1.9419, 51.8318, 103.3493, 152.7052, 199.9903, 245.2914],\n",
+    "        \"close\": [1.9419, 57.9719, 111.6510, 163.0778, 212.3468, 259.5486]\n",
+    "    },\n",
+    "    \"FOC_1\": {\n",
+    "        \"frequency\": 42.0,\n",
+    "        \"phase\": 81.303297,\n",
+    "        \"distance\": 8.4,\n",
+    "        \"open\": [-5.1362, 42.5536, 88.2425, 132.0144, 173.9497, 216.7867],\n",
+    "        \"close\": [5.1362, 54.2095, 101.2237, 146.2653, 189.417, 230.7582]\n",
+    "    },\n",
+    "    \"BP_1\": {\n",
+    "        \"frequency\": 7.0,\n",
+    "        \"phase\": 31.080,\n",
+    "        \"distance\": 8.45,\n",
+    "        \"open\": [-23.6029],\n",
+    "        \"close\": [23.6029]\n",
+    "    },\n",
+    "    \"FOC_2\": {\n",
+    "        \"frequency\": 42.0,\n",
+    "        \"phase\": 107.013442,\n",
+    "        \"distance\": 12.2,\n",
+    "        \"open\": [-16.3227, 53.7401, 120.8633, 185.1701, 246.7787, 307.0165],\n",
+    "        \"close\": [16.3227, 86.8303, 154.3794, 218.7551, 280.7508, 340.3188]\n",
+    "    },\n",
+    "    \"BP_2\": {\n",
+    "        \"frequency\": 7.0,\n",
+    "        \"phase\": 44.224,\n",
+    "        \"distance\": 12.25,\n",
+    "        \"open\": [-34.4663],\n",
+    "        \"close\": [34.4663]\n",
+    "    },\n",
+    "    \"T0_alpha\": {\n",
+    "        \"frequency\": 14.0,\n",
+    "        \"phase\": 179.672,\n",
+    "        \"distance\": 13.5,\n",
+    "        \"open\": [-167.8986],\n",
+    "        \"close\": [167.8986]\n",
+    "    },\n",
+    "    \"T0_beta\": {\n",
+    "        \"frequency\": 14.0,\n",
+    "        \"phase\": 179.672,\n",
+    "        \"distance\": 13.7,\n",
+    "        \"open\": [-167.8986],\n",
+    "        \"close\": [167.8986]\n",
+    "    },\n",
+    "    \"FOC_3\": {\n",
+    "        \"frequency\": 28.0,\n",
+    "        \"phase\": 92.993,\n",
+    "        \"distance\": 17.0,\n",
+    "        \"open\": [-20.302, 45.247, 108.0457, 168.2095, 225.8489, 282.2199],\n",
+    "        \"close\": [20.302, 85.357, 147.6824, 207.3927, 264.5977, 319.4024]\n",
+    "    },\n",
+    "    \"FOC_4\": {\n",
+    "        \"frequency\": 14.0,\n",
+    "        \"phase\": 61.584,\n",
+    "        \"distance\": 23.69,\n",
+    "        \"open\": [-16.7157, 29.1882, 73.1661, 115.2988, 155.6636, 195.5254],\n",
+    "        \"close\": [16.7157, 61.8217, 105.0352, 146.4355, 186.0987, 224.0978]\n",
+    "    },\n",
+    "    \"FOC_5\": {\n",
+    "        \"frequency\": 14.0,\n",
+    "        \"phase\": 82.581,\n",
+    "        \"distance\": 33.0,\n",
+    "        \"open\": [-25.8514, 38.3239, 99.8064, 160.1254, 217.4321, 272.5426],\n",
+    "        \"close\": [25.8514, 88.4621, 147.4729, 204.0245, 257.7603, 313.7139]\n",
+    "    },\n",
+    "\n",
+    "}\n",
+    "\n",
+    "choppers = [\n",
+    "    tof.Chopper(\n",
+    "        frequency=ch[\"frequency\"] * Hz,\n",
+    "        direction=tof.Clockwise,\n",
+    "        open=sc.array(dims=[\"cutout\"], values=ch[\"open\"], unit=\"deg\"),\n",
+    "        close=sc.array(dims=[\"cutout\"], values=ch[\"close\"], unit=\"deg\"),\n",
+    "        phase=ch[\"phase\"] * deg,\n",
+    "        distance=ch[\"distance\"] * meter,\n",
+    "        name=key,\n",
+    "    )\n",
+    "    for key, ch in parameters.items()\n",
+    "]"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "6",
+   "metadata": {},
+   "source": [
+    "### Detector\n",
+    "\n",
+    "ODIN has a single detector panel 60.5 meters from the source."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "7",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "detectors = [\n",
+    "    tof.Detector(distance=60.5 * meter, name=\"detector\"),\n",
+    "]"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "8",
+   "metadata": {},
+   "source": [
+    "## Run the simulation\n",
+    "\n",
+    "We propagate our pulse of neutrons through the chopper cascade and inspect the results."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "9",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "model = tof.Model(source=source, choppers=choppers, detectors=detectors)\n",
+    "results = model.run()\n",
+    "results.plot(blocked_rays=5000)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "10",
+   "metadata": {},
+   "source": [
+    "We can see that the chopper cascade is implementing WFM and pulse-skipping at the same time!\n",
+    "\n",
+    "## Wavelength as a function of time-of-arrival\n",
+    "\n",
+    "### Plotting wavelength vs time-of-arrival\n",
+    "\n",
+    "Since we know the true wavelength of our neutrons,\n",
+    "as well as the time at which the neutrons arrive at the detector\n",
+    "(coordinate named `toa` in the detector reading),\n",
+    "we can plot an image of the wavelengths as a function of time-of-arrival:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "11",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Squeeze the pulse dimension since we only have one pulse\n",
+    "events = results['detector'].data.flatten(to='event')\n",
+    "# Remove the events that don't make it to the detector\n",
+    "events = events[~events.masks['blocked_by_others']]\n",
+    "# Histogram and plot\n",
+    "events.hist(wavelength=500, toa=500).plot(norm='log', grid=True)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "12",
+   "metadata": {},
+   "source": [
+    "### Defining a conversion from `toa` to `wavelength`\n",
+    "\n",
+    "The image above shows that there is a pretty tight correlation between time-of-arrival and wavelength.\n",
+    "\n",
+    "We compute the mean wavelength inside a given `toa` bin to define a relation between `toa` and `wavelength`."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "13",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "binned = events.bin(toa=500)\n",
+    "\n",
+    "# Weighted mean of wavelength inside each bin\n",
+    "mu = (\n",
+    "    binned.bins.data * binned.bins.coords['wavelength']\n",
+    ").bins.sum() / binned.bins.sum()\n",
+    "\n",
+    "# Variance of wavelengths inside each bin\n",
+    "var = (\n",
+    "    binned.bins.data * (binned.bins.coords['wavelength'] - mu) ** 2\n",
+    ") / binned.bins.sum()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "14",
+   "metadata": {},
+   "source": [
+    "We can now overlay our mean wavelength function on the image above:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "15",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import matplotlib.pyplot as plt\n",
+    "\n",
+    "fig, ax = plt.subplots(2, 1)\n",
+    "\n",
+    "f = events.hist(wavelength=500, tof=500).plot(norm='log', cbar=False, ax=ax[0])\n",
+    "mu.name = 'Wavelength'\n",
+    "mu.plot(ax=ax[0], color='C1', grid=True)\n",
+    "stddev = sc.sqrt(var.hist())\n",
+    "stddev.name = 'Standard deviation'\n",
+    "stddev.plot(ax=ax[1], grid=True)\n",
+    "fig.set_size_inches(6, 8)\n",
+    "fig.tight_layout()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "16",
+   "metadata": {},
+   "source": [
+    "## Computing wavelengths\n",
+    "\n",
+    "We set up an interpolator that will compute wavelengths given an array of `toas`."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "17",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from scipp.scipy.interpolate import interp1d\n",
+    "\n",
+    "# Set up interpolator\n",
+    "y = mu.copy()\n",
+    "y.coords['toa'] = sc.midpoints(y.coords['toa'])\n",
+    "f = interp1d(y, 'toa', bounds_error=False)\n",
+    "\n",
+    "# Compute wavelengths\n",
+    "wavs = f(events.coords['toa'].rename_dims(event='toa'))\n",
+    "wavelengths = sc.DataArray(\n",
+    "    data=sc.ones(sizes=wavs.sizes, unit='counts'), coords={'wavelength': wavs.data}\n",
+    ").rename_dims(toa='event')\n",
+    "wavelengths"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "18",
+   "metadata": {},
+   "source": [
+    "We can now compare our computed wavelengths to the true wavelengths of the neutrons:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "19",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "pp.plot(\n",
+    "    {\n",
+    "        'wfm': wavelengths.hist(wavelength=300),\n",
+    "        'original': events.hist(wavelength=300),\n",
+    "    }\n",
+    ")"
+   ]
+  }
+ ],
+ "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"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}