DREAM in WFM mode#
This is a simulation of the DREAM 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')
AA = sc.Unit('angstrom')
Create a source#
We first create an ESS source with 2 pulses containing 500,000 neutrons each.
[2]:
source = tof.Source(facility='ess', neutrons=500_000, pulses=2)
source.plot()
[2]:
Component set-up#
Choppers#
The DREAM chopper cascade consists of:
Two counter-rotating pulse-shaping choppers (PSC) that are very close to each other, located ~6m from the source
An overlap chopper placed right after the two PSCs
A band control chopper
A T0 chopper
[3]:
choppers = [
tof.Chopper(
frequency=14 * Hz,
direction=tof.AntiClockwise,
centers=sc.array(
dims=['cutout'],
values=[0, 72, 86.4, 115.2, 172.8, 273.6, 288.0, 302.4],
unit='deg',
),
widths=sc.array(
dims=['cutout'],
values=[2.46, 3.02, 3.27, 3.27, 5.02, 3.93, 3.93, 2.46],
unit='deg',
),
phase=(286 - 180) * deg,
distance=6.145 * meter,
name="PSC1",
),
tof.Chopper(
frequency=14 * Hz,
direction=tof.Clockwise,
centers=sc.array(
dims=['cutout'],
values=[0, 28.8, 57.6, 144, 158.4, 216, 259.2, 316.8],
unit='deg',
),
widths=sc.array(
dims=['cutout'],
values=[2.46, 3.60, 3.60, 3.23, 3.27, 3.77, 3.94, 2.62],
unit='deg',
),
phase=236 * deg,
distance=6.155 * meter,
name="PSC2",
),
tof.Chopper(
frequency=14 * Hz,
direction=tof.AntiClockwise,
centers=sc.array(
dims=['cutout'],
values=[0.0],
unit='deg',
),
widths=sc.array(
dims=['cutout'],
values=[27.6],
unit='deg',
),
phase=(297 - 180 - 90) * deg,
distance=6.174 * meter,
name="OC",
),
tof.Chopper(
frequency=112 * Hz,
direction=tof.AntiClockwise,
centers=sc.array(
dims=['cutout'],
values=[0.0, 180.0],
unit='deg',
),
widths=sc.array(
dims=['cutout'],
values=[73.75, 73.75],
unit='deg',
),
phase=(240 - 180) * deg,
distance=9.78 * meter,
name="BC",
),
tof.Chopper(
frequency=28 * Hz,
direction=tof.AntiClockwise,
centers=sc.array(
dims=['cutout'],
values=[0.0],
unit='deg',
),
widths=sc.array(
dims=['cutout'],
values=[314.9],
unit='deg',
),
phase=(280 - 180) * deg,
distance=13.05 * meter,
name="T0",
),
]
Detector banks and monitors#
DREAM has 5 detector banks: the Mantle, two End-caps, a High-resolution detector and a SANS detector.
For each detector bank, we use a single mean distance (in practice, one could have a different distance for each pixel).
[4]:
sample_position = 76.55 * meter
detectors = [
tof.Detector(distance=sample_position + 1.125 * meter, name='mantle'),
tof.Detector(distance=sample_position + 1.125 * meter, name='end-cap'),
tof.Detector(distance=sample_position + 2.5 * meter, name='high-resolution'),
tof.Detector(distance=sample_position + 2.5 * meter, name='sans'),
]
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>)
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]:
events = sc.DataGroup()
for key, da in results.detectors.items():
bank = da.data.flatten(to='event')
events[key] = bank[~bank.masks['blocked_by_others']]
# Histogram and plot
events['mantle'].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(tof=500)
# Weighted mean of wavelength inside each bin
mu = sc.DataGroup(
{
key: (da.bins.data * da.bins.coords['wavelength']).bins.sum() / da.bins.sum()
for key, da in binned.items()
}
)
mu.plot(grid=True)
[7]:
Computing wavelengths#
We set up an interpolator that will compute wavelengths given an array of toas.
[8]:
from scipp.scipy.interpolate import interp1d
wavelengths = sc.DataGroup()
for key in mu:
# Set up interpolator
y = mu[key].copy()
y.coords['tof'] = sc.midpoints(y.coords['tof'])
f = interp1d(y, 'tof', bounds_error=False)
# Compute wavelengths
wavs = f(events[key].coords['tof'].rename_dims(event='tof'))
wavelengths[key] = sc.DataArray(
data=sc.ones(sizes=wavs.sizes, unit='counts'), coords={'wavelength': wavs.data}
).rename_dims(tof='event')
wavelengths
[8]:
- mantlescippDataArray(event: 90263)float64counts1.0, 1.0, ..., 1.0, 1.0
- end-capscippDataArray(event: 90263)float64counts1.0, 1.0, ..., 1.0, 1.0
- high-resolutionscippDataArray(event: 90263)float64counts1.0, 1.0, ..., 1.0, 1.0
- sansscippDataArray(event: 90263)float64counts1.0, 1.0, ..., 1.0, 1.0
We can now compare our computed wavelengths to the true wavelengths of the neutrons.
[9]:
pp.plot(
{
'wfm': wavelengths['mantle'].hist(wavelength=300),
'original': events['mantle'].hist(wavelength=300),
}
)
[9]:
