Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Odin notebook #66

Merged
merged 3 commits into from
Nov 27, 2024
Merged

Odin notebook #66

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
24ad558
add first draft of odin notebook
nvaytet Nov 25, 2024
7600cc5
Merge branch 'main' into odin-notebook
nvaytet Nov 26, 2024
8f91616
fix text in notebook
nvaytet Nov 26, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions docs/ess/index.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,3 +5,4 @@ ESS instruments
:maxdepth: 2

dream
odin
385 changes: 385 additions & 0 deletions docs/ess/odin.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,385 @@
{
"cells": [
{
"cell_type": "markdown",
"id": "6b645603-12c0-4b3d-ab3c-529ad4f373e3",
"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": "899a245a-0ae7-45dd-b684-3b9d1c034cc6",
"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": "3890d24c-93ab-4faa-b653-c9208b5eba23",
"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": "872800a5-4cf7-424e-8d2d-0db81292456c",
"metadata": {},
"outputs": [],
"source": [
"source = tof.Source(facility=\"ess\", neutrons=500_000, pulses=4)\n",
"source.plot()"
]
},
{
"cell_type": "markdown",
"id": "b5ffdab6-10cd-4719-8ff2-86e5e1d10571",
"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": "34314bcc-978c-468c-b5a1-d84ce33e53d5",
"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": "bd6e1399-7f05-4d5c-83e3-7d249d9e0a61",
"metadata": {},
"source": [
"### Detector\n",
"\n",
"ODIN has a single detector panel 60.5 meters from the source."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "a377e75c-51a5-4994-af9c-92139d27bee3",
"metadata": {},
"outputs": [],
"source": [
"detectors = [\n",
" tof.Detector(distance=60.5 * meter, name=\"detector\"),\n",
"]"
]
},
{
"cell_type": "markdown",
"id": "4deb2f50-92b7-4838-b180-17e3ce43743d",
"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": "e5812eda-2e56-4f72-89bc-4e778046243f",
"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": "51e45cee-a630-4886-bc9f-a744a20cba28",
"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": "c3f48d70-7836-4e6a-aa9c-e278303af30c",
"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": "aadc3250-f49d-462f-8aac-2d3500596b6b",
"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": "9e8f056e-8c8e-4c1b-8dbf-5b350bba5c94",
"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": "fb5cddc2-9ebf-427f-b38f-cf25b2c5cdfc",
"metadata": {},
"source": [
"We can now overlay our mean wavelength function on the image above:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "b0a83ab4-af22-40a4-851a-4a80ce358f45",
"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": "d3675286-e36e-4d65-9d6c-012f5b4c39e5",
"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": "0f231769-8e79-4692-b4f0-a1e9e11782e1",
"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": "a6512644-d9b7-47e1-861b-9dde24fa012a",
"metadata": {},
"source": [
"We can now compare our computed wavelengths to the true wavelengths of the neutrons:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "1f4aad58-a1b0-4c90-9215-ab2e52191bbb",
"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
}