diff --git a/docs/ess/index.rst b/docs/ess/index.rst index 6cd0de3..402bf12 100644 --- a/docs/ess/index.rst +++ b/docs/ess/index.rst @@ -5,3 +5,4 @@ ESS instruments :maxdepth: 2 dream + odin diff --git a/docs/ess/odin.ipynb b/docs/ess/odin.ipynb new file mode 100644 index 0000000..0871599 --- /dev/null +++ b/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 +}