Skip to content

Commit f574d3b

Browse files
nvaytetnvaytet
authored
Merge pull request #66 from scipp/odin-notebook
Odin notebook
2 parents ec9a76f + 8f91616 commit f574d3b

File tree

2 files changed

+386
-0
lines changed

2 files changed

+386
-0
lines changed

‎docs/ess/index.rst

+1
Original file line numberDiff line numberDiff line change
@@ -5,3 +5,4 @@ ESS instruments
55
:maxdepth: 2
66

77
dream
8+
odin

‎docs/ess/odin.ipynb

+385
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,385 @@
1+
{
2+
"cells": [
3+
{
4+
"cell_type": "markdown",
5+
"id": "6b645603-12c0-4b3d-ab3c-529ad4f373e3",
6+
"metadata": {},
7+
"source": [
8+
"# ODIN in WFM mode\n",
9+
"\n",
10+
"This is a simulation of the ODIN chopper cascade in WFM mode.\n",
11+
"We also show how one can convert the neutron arrival times at the detector to wavelength."
12+
]
13+
},
14+
{
15+
"cell_type": "code",
16+
"execution_count": null,
17+
"id": "899a245a-0ae7-45dd-b684-3b9d1c034cc6",
18+
"metadata": {},
19+
"outputs": [],
20+
"source": [
21+
"import numpy as np\n",
22+
"import scipp as sc\n",
23+
"import plopp as pp\n",
24+
"import tof\n",
25+
"\n",
26+
"Hz = sc.Unit(\"Hz\")\n",
27+
"deg = sc.Unit(\"deg\")\n",
28+
"meter = sc.Unit(\"m\")"
29+
]
30+
},
31+
{
32+
"cell_type": "markdown",
33+
"id": "3890d24c-93ab-4faa-b653-c9208b5eba23",
34+
"metadata": {},
35+
"source": [
36+
"## Create a source pulse\n",
37+
"\n",
38+
"We first create a source with 4 pulses containing 800,000 neutrons each,\n",
39+
"and whose distribution follows the ESS time and wavelength profiles (both thermal and cold neutrons are included)."
40+
]
41+
},
42+
{
43+
"cell_type": "code",
44+
"execution_count": null,
45+
"id": "872800a5-4cf7-424e-8d2d-0db81292456c",
46+
"metadata": {},
47+
"outputs": [],
48+
"source": [
49+
"source = tof.Source(facility=\"ess\", neutrons=500_000, pulses=4)\n",
50+
"source.plot()"
51+
]
52+
},
53+
{
54+
"cell_type": "markdown",
55+
"id": "b5ffdab6-10cd-4719-8ff2-86e5e1d10571",
56+
"metadata": {},
57+
"source": [
58+
"## Component set-up\n",
59+
"\n",
60+
"### Choppers\n",
61+
"\n",
62+
"The ODIN chopper cascade consists of:\n",
63+
"\n",
64+
"- 2 WFM choppers\n",
65+
"- 5 frame-overlap choppers\n",
66+
"- 2 band-control choppers\n",
67+
"- 1 T0 chopper"
68+
]
69+
},
70+
{
71+
"cell_type": "code",
72+
"execution_count": null,
73+
"id": "34314bcc-978c-468c-b5a1-d84ce33e53d5",
74+
"metadata": {},
75+
"outputs": [],
76+
"source": [
77+
"parameters = {\n",
78+
" \"WFMC_1\": {\n",
79+
" \"frequency\": 56.0,\n",
80+
" \"phase\": 93.244,\n",
81+
" \"distance\": 6.85,\n",
82+
" \"open\": [-1.9419, 49.5756, 98.9315, 146.2165, 191.5176, 234.9179],\n",
83+
" \"close\": [1.9419, 55.7157, 107.2332, 156.5891, 203.8741, 249.1752]\n",
84+
" },\n",
85+
" \"WFMC_2\": {\n",
86+
" \"frequency\": 56.0,\n",
87+
" \"phase\": 97.128,\n",
88+
" \"distance\": 7.15,\n",
89+
" \"open\": [-1.9419, 51.8318, 103.3493, 152.7052, 199.9903, 245.2914],\n",
90+
" \"close\": [1.9419, 57.9719, 111.6510, 163.0778, 212.3468, 259.5486]\n",
91+
" },\n",
92+
" \"FOC_1\": {\n",
93+
" \"frequency\": 42.0,\n",
94+
" \"phase\": 81.303297,\n",
95+
" \"distance\": 8.4,\n",
96+
" \"open\": [-5.1362, 42.5536, 88.2425, 132.0144, 173.9497, 216.7867],\n",
97+
" \"close\": [5.1362, 54.2095, 101.2237, 146.2653, 189.417, 230.7582]\n",
98+
" },\n",
99+
" \"BP_1\": {\n",
100+
" \"frequency\": 7.0,\n",
101+
" \"phase\": 31.080,\n",
102+
" \"distance\": 8.45,\n",
103+
" \"open\": [-23.6029],\n",
104+
" \"close\": [23.6029]\n",
105+
" },\n",
106+
" \"FOC_2\": {\n",
107+
" \"frequency\": 42.0,\n",
108+
" \"phase\": 107.013442,\n",
109+
" \"distance\": 12.2,\n",
110+
" \"open\": [-16.3227, 53.7401, 120.8633, 185.1701, 246.7787, 307.0165],\n",
111+
" \"close\": [16.3227, 86.8303, 154.3794, 218.7551, 280.7508, 340.3188]\n",
112+
" },\n",
113+
" \"BP_2\": {\n",
114+
" \"frequency\": 7.0,\n",
115+
" \"phase\": 44.224,\n",
116+
" \"distance\": 12.25,\n",
117+
" \"open\": [-34.4663],\n",
118+
" \"close\": [34.4663]\n",
119+
" },\n",
120+
" \"T0_alpha\": {\n",
121+
" \"frequency\": 14.0,\n",
122+
" \"phase\": 179.672,\n",
123+
" \"distance\": 13.5,\n",
124+
" \"open\": [-167.8986],\n",
125+
" \"close\": [167.8986]\n",
126+
" },\n",
127+
" \"T0_beta\": {\n",
128+
" \"frequency\": 14.0,\n",
129+
" \"phase\": 179.672,\n",
130+
" \"distance\": 13.7,\n",
131+
" \"open\": [-167.8986],\n",
132+
" \"close\": [167.8986]\n",
133+
" },\n",
134+
" \"FOC_3\": {\n",
135+
" \"frequency\": 28.0,\n",
136+
" \"phase\": 92.993,\n",
137+
" \"distance\": 17.0,\n",
138+
" \"open\": [-20.302, 45.247, 108.0457, 168.2095, 225.8489, 282.2199],\n",
139+
" \"close\": [20.302, 85.357, 147.6824, 207.3927, 264.5977, 319.4024]\n",
140+
" },\n",
141+
" \"FOC_4\": {\n",
142+
" \"frequency\": 14.0,\n",
143+
" \"phase\": 61.584,\n",
144+
" \"distance\": 23.69,\n",
145+
" \"open\": [-16.7157, 29.1882, 73.1661, 115.2988, 155.6636, 195.5254],\n",
146+
" \"close\": [16.7157, 61.8217, 105.0352, 146.4355, 186.0987, 224.0978]\n",
147+
" },\n",
148+
" \"FOC_5\": {\n",
149+
" \"frequency\": 14.0,\n",
150+
" \"phase\": 82.581,\n",
151+
" \"distance\": 33.0,\n",
152+
" \"open\": [-25.8514, 38.3239, 99.8064, 160.1254, 217.4321, 272.5426],\n",
153+
" \"close\": [25.8514, 88.4621, 147.4729, 204.0245, 257.7603, 313.7139]\n",
154+
" },\n",
155+
"\n",
156+
"}\n",
157+
"\n",
158+
"choppers = [\n",
159+
" tof.Chopper(\n",
160+
" frequency=ch[\"frequency\"] * Hz,\n",
161+
" direction=tof.Clockwise,\n",
162+
" open=sc.array(dims=[\"cutout\"], values=ch[\"open\"], unit=\"deg\"),\n",
163+
" close=sc.array(dims=[\"cutout\"], values=ch[\"close\"], unit=\"deg\"),\n",
164+
" phase=ch[\"phase\"] * deg,\n",
165+
" distance=ch[\"distance\"] * meter,\n",
166+
" name=key,\n",
167+
" )\n",
168+
" for key, ch in parameters.items()\n",
169+
"]"
170+
]
171+
},
172+
{
173+
"cell_type": "markdown",
174+
"id": "bd6e1399-7f05-4d5c-83e3-7d249d9e0a61",
175+
"metadata": {},
176+
"source": [
177+
"### Detector\n",
178+
"\n",
179+
"ODIN has a single detector panel 60.5 meters from the source."
180+
]
181+
},
182+
{
183+
"cell_type": "code",
184+
"execution_count": null,
185+
"id": "a377e75c-51a5-4994-af9c-92139d27bee3",
186+
"metadata": {},
187+
"outputs": [],
188+
"source": [
189+
"detectors = [\n",
190+
" tof.Detector(distance=60.5 * meter, name=\"detector\"),\n",
191+
"]"
192+
]
193+
},
194+
{
195+
"cell_type": "markdown",
196+
"id": "4deb2f50-92b7-4838-b180-17e3ce43743d",
197+
"metadata": {},
198+
"source": [
199+
"## Run the simulation\n",
200+
"\n",
201+
"We propagate our pulse of neutrons through the chopper cascade and inspect the results."
202+
]
203+
},
204+
{
205+
"cell_type": "code",
206+
"execution_count": null,
207+
"id": "e5812eda-2e56-4f72-89bc-4e778046243f",
208+
"metadata": {},
209+
"outputs": [],
210+
"source": [
211+
"model = tof.Model(source=source, choppers=choppers, detectors=detectors)\n",
212+
"results = model.run()\n",
213+
"results.plot(blocked_rays=5000)"
214+
]
215+
},
216+
{
217+
"cell_type": "markdown",
218+
"id": "51e45cee-a630-4886-bc9f-a744a20cba28",
219+
"metadata": {},
220+
"source": [
221+
"We can see that the chopper cascade is implementing WFM and pulse-skipping at the same time!\n",
222+
"\n",
223+
"## Wavelength as a function of time-of-arrival\n",
224+
"\n",
225+
"### Plotting wavelength vs time-of-arrival\n",
226+
"\n",
227+
"Since we know the true wavelength of our neutrons,\n",
228+
"as well as the time at which the neutrons arrive at the detector\n",
229+
"(coordinate named `toa` in the detector reading),\n",
230+
"we can plot an image of the wavelengths as a function of time-of-arrival:"
231+
]
232+
},
233+
{
234+
"cell_type": "code",
235+
"execution_count": null,
236+
"id": "c3f48d70-7836-4e6a-aa9c-e278303af30c",
237+
"metadata": {},
238+
"outputs": [],
239+
"source": [
240+
"# Squeeze the pulse dimension since we only have one pulse\n",
241+
"events = results['detector'].data.flatten(to='event')\n",
242+
"# Remove the events that don't make it to the detector\n",
243+
"events = events[~events.masks['blocked_by_others']]\n",
244+
"# Histogram and plot\n",
245+
"events.hist(wavelength=500, toa=500).plot(norm='log', grid=True)"
246+
]
247+
},
248+
{
249+
"cell_type": "markdown",
250+
"id": "aadc3250-f49d-462f-8aac-2d3500596b6b",
251+
"metadata": {},
252+
"source": [
253+
"### Defining a conversion from `toa` to `wavelength`\n",
254+
"\n",
255+
"The image above shows that there is a pretty tight correlation between time-of-arrival and wavelength.\n",
256+
"\n",
257+
"We compute the mean wavelength inside a given `toa` bin to define a relation between `toa` and `wavelength`."
258+
]
259+
},
260+
{
261+
"cell_type": "code",
262+
"execution_count": null,
263+
"id": "9e8f056e-8c8e-4c1b-8dbf-5b350bba5c94",
264+
"metadata": {},
265+
"outputs": [],
266+
"source": [
267+
"binned = events.bin(toa=500)\n",
268+
"\n",
269+
"# Weighted mean of wavelength inside each bin\n",
270+
"mu = (\n",
271+
" binned.bins.data * binned.bins.coords['wavelength']\n",
272+
").bins.sum() / binned.bins.sum()\n",
273+
"\n",
274+
"# Variance of wavelengths inside each bin\n",
275+
"var = (\n",
276+
" binned.bins.data * (binned.bins.coords['wavelength'] - mu) ** 2\n",
277+
") / binned.bins.sum()"
278+
]
279+
},
280+
{
281+
"cell_type": "markdown",
282+
"id": "fb5cddc2-9ebf-427f-b38f-cf25b2c5cdfc",
283+
"metadata": {},
284+
"source": [
285+
"We can now overlay our mean wavelength function on the image above:"
286+
]
287+
},
288+
{
289+
"cell_type": "code",
290+
"execution_count": null,
291+
"id": "b0a83ab4-af22-40a4-851a-4a80ce358f45",
292+
"metadata": {},
293+
"outputs": [],
294+
"source": [
295+
"import matplotlib.pyplot as plt\n",
296+
"\n",
297+
"fig, ax = plt.subplots(2, 1)\n",
298+
"\n",
299+
"f = events.hist(wavelength=500, tof=500).plot(norm='log', cbar=False, ax=ax[0])\n",
300+
"mu.name = 'Wavelength'\n",
301+
"mu.plot(ax=ax[0], color='C1', grid=True)\n",
302+
"stddev = sc.sqrt(var.hist())\n",
303+
"stddev.name = 'Standard deviation'\n",
304+
"stddev.plot(ax=ax[1], grid=True)\n",
305+
"fig.set_size_inches(6, 8)\n",
306+
"fig.tight_layout()"
307+
]
308+
},
309+
{
310+
"cell_type": "markdown",
311+
"id": "d3675286-e36e-4d65-9d6c-012f5b4c39e5",
312+
"metadata": {},
313+
"source": [
314+
"## Computing wavelengths\n",
315+
"\n",
316+
"We set up an interpolator that will compute wavelengths given an array of `toas`."
317+
]
318+
},
319+
{
320+
"cell_type": "code",
321+
"execution_count": null,
322+
"id": "0f231769-8e79-4692-b4f0-a1e9e11782e1",
323+
"metadata": {},
324+
"outputs": [],
325+
"source": [
326+
"from scipp.scipy.interpolate import interp1d\n",
327+
"\n",
328+
"# Set up interpolator\n",
329+
"y = mu.copy()\n",
330+
"y.coords['toa'] = sc.midpoints(y.coords['toa'])\n",
331+
"f = interp1d(y, 'toa', bounds_error=False)\n",
332+
"\n",
333+
"# Compute wavelengths\n",
334+
"wavs = f(events.coords['toa'].rename_dims(event='toa'))\n",
335+
"wavelengths = sc.DataArray(\n",
336+
" data=sc.ones(sizes=wavs.sizes, unit='counts'), coords={'wavelength': wavs.data}\n",
337+
").rename_dims(toa='event')\n",
338+
"wavelengths"
339+
]
340+
},
341+
{
342+
"cell_type": "markdown",
343+
"id": "a6512644-d9b7-47e1-861b-9dde24fa012a",
344+
"metadata": {},
345+
"source": [
346+
"We can now compare our computed wavelengths to the true wavelengths of the neutrons:"
347+
]
348+
},
349+
{
350+
"cell_type": "code",
351+
"execution_count": null,
352+
"id": "1f4aad58-a1b0-4c90-9215-ab2e52191bbb",
353+
"metadata": {},
354+
"outputs": [],
355+
"source": [
356+
"pp.plot(\n",
357+
" {\n",
358+
" 'wfm': wavelengths.hist(wavelength=300),\n",
359+
" 'original': events.hist(wavelength=300),\n",
360+
" }\n",
361+
")"
362+
]
363+
}
364+
],
365+
"metadata": {
366+
"kernelspec": {
367+
"display_name": "Python 3 (ipykernel)",
368+
"language": "python",
369+
"name": "python3"
370+
},
371+
"language_info": {
372+
"codemirror_mode": {
373+
"name": "ipython",
374+
"version": 3
375+
},
376+
"file_extension": ".py",
377+
"mimetype": "text/x-python",
378+
"name": "python",
379+
"nbconvert_exporter": "python",
380+
"pygments_lexer": "ipython3"
381+
}
382+
},
383+
"nbformat": 4,
384+
"nbformat_minor": 5
385+
}

0 commit comments

Comments
 (0)