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