Methods for creating wavelength lookup tables#
In this notebook, we go through the different methods that can be used to construct a wavelength lookup table, which is what predicts the most likely wavelength of a neutron arriving at the detector, given its arrival time event_time_offset.
The three methods are the following:
'simulation': create the lookup table by simulating individual neutrons traveling through the chopper system using thetofMonte-Carlo package.'analytical': create the lookup table using analytical formulas to propagate and chop a pulse of neutrons through the chopper cascade.'file': read a pre-computed lookup table from a file.
[1]:
import scipp as sc
import plopp as pp
import scippnexus as snx
from scippneutron.chopper import DiskChopper
from ess.reduce import unwrap
from ess.reduce.nexus.types import AnyRun, Position
Beamline setup#
[2]:
psc1 = DiskChopper(
frequency=sc.scalar(14.0, unit="Hz"),
beam_position=sc.scalar(0.0, unit="deg"),
phase=sc.scalar(286 - 180, unit="deg"),
axle_position=sc.vector(value=[0, 0, -70.405], unit="m"),
slit_begin=sc.array(
dims=["cutout"],
values=[-1.23, 70.49, 84.765, 113.565, 170.29, 271.635, 286.035, 301.17],
unit="deg",
),
slit_end=sc.array(
dims=["cutout"],
values=[1.23, 73.51, 88.035, 116.835, 175.31, 275.565, 289.965, 303.63],
unit="deg",
),
slit_height=sc.scalar(10.0, unit="cm"),
radius=sc.scalar(30.0, unit="cm"),
)
psc2 = DiskChopper(
frequency=sc.scalar(-14.0, unit="Hz"),
beam_position=sc.scalar(0.0, unit="deg"),
phase=sc.scalar(-236, unit="deg"),
axle_position=sc.vector(value=[0, 0, -70.395], unit="m"),
slit_begin=sc.array(
dims=["cutout"],
values=[-1.23, 27.0, 55.8, 142.385, 156.765, 214.115, 257.23, 315.49],
unit="deg",
),
slit_end=sc.array(
dims=["cutout"],
values=[1.23, 30.6, 59.4, 145.615, 160.035, 217.885, 261.17, 318.11],
unit="deg",
),
slit_height=sc.scalar(10.0, unit="cm"),
radius=sc.scalar(30.0, unit="cm"),
)
oc = DiskChopper(
frequency=sc.scalar(14.0, unit="Hz"),
beam_position=sc.scalar(0.0, unit="deg"),
phase=sc.scalar(297 - 180 - 90, unit="deg"),
axle_position=sc.vector(value=[0, 0, -70.376], unit="m"),
slit_begin=sc.array(dims=["cutout"], values=[-27.6 * 0.5], unit="deg"),
slit_end=sc.array(dims=["cutout"], values=[27.6 * 0.5], unit="deg"),
slit_height=sc.scalar(10.0, unit="cm"),
radius=sc.scalar(30.0, unit="cm"),
)
bcc = DiskChopper(
frequency=sc.scalar(112.0, unit="Hz"),
beam_position=sc.scalar(0.0, unit="deg"),
phase=sc.scalar(240 - 180, unit="deg"),
axle_position=sc.vector(value=[0, 0, -66.77], unit="m"),
slit_begin=sc.array(dims=["cutout"], values=[-36.875, 143.125], unit="deg"),
slit_end=sc.array(dims=["cutout"], values=[36.875, 216.875], unit="deg"),
slit_height=sc.scalar(10.0, unit="cm"),
radius=sc.scalar(30.0, unit="cm"),
)
t0 = DiskChopper(
frequency=sc.scalar(28.0, unit="Hz"),
beam_position=sc.scalar(0.0, unit="deg"),
phase=sc.scalar(280 - 180, unit="deg"),
axle_position=sc.vector(value=[0, 0, -63.5], unit="m"),
slit_begin=sc.array(dims=["cutout"], values=[-314.9 * 0.5], unit="deg"),
slit_end=sc.array(dims=["cutout"], values=[314.9 * 0.5], unit="deg"),
slit_height=sc.scalar(10.0, unit="cm"),
radius=sc.scalar(30.0, unit="cm"),
)
disk_choppers = {"psc1": psc1, "psc2": psc2, "oc": oc, "bcc": bcc, "t0": t0}
Method: ‘simulation’#
This first method creates the lookup table by simulating individual neutrons traveling through the chopper system using the tof Monte-Carlo package.
This is slower but can be more accurate if the spread in neutron wavelengths is large at the detector.
Other simulation software such as McStas can also be used to replace tof for even better results.
[3]:
wf = unwrap.GenericUnwrapWorkflow(
run_types=[AnyRun], monitor_types=[], wavelength_from="simulation"
)
wf[unwrap.LtotalRange[AnyRun, snx.NXdetector]] = (
sc.scalar(50.0, unit="m"),
sc.scalar(60.0, unit="m"),
)
wf[unwrap.NumberOfSimulatedNeutrons] = (
200_000 # Increase this number for more reliable results
)
wf[Position[snx.NXsource, AnyRun]] = sc.vector([0, 0, -76.55], unit="m")
wf[unwrap.DiskChoppers[AnyRun]] = disk_choppers
wf[unwrap.DistanceResolution] = sc.scalar(0.1, unit="m")
wf[unwrap.TimeResolution] = sc.scalar(250.0, unit='us')
wf[unwrap.PulsePeriod] = 1.0 / sc.scalar(14.0, unit="Hz")
wf[unwrap.PulseStrideOffset] = None
table = wf.compute(unwrap.LookupTable[AnyRun, snx.NXdetector])
table.array
[3]:
scipp.DataArray (475.04 KB)
- distance: 105
- event_time_offset: 287
- distance(distance)float64m49.8, 49.9, ..., 60.100, 60.200
Values:
array([49.8, 49.9, 50. , 50.1, 50.2, 50.3, 50.4, 50.5, 50.6, 50.7, 50.8, 50.9, 51. , 51.1, 51.2, 51.3, 51.4, 51.5, 51.6, 51.7, 51.8, 51.9, 52. , 52.1, 52.2, 52.3, 52.4, 52.5, 52.6, 52.7, 52.8, 52.9, 53. , 53.1, 53.2, 53.3, 53.4, 53.5, 53.6, 53.7, 53.8, 53.9, 54. , 54.1, 54.2, 54.3, 54.4, 54.5, 54.6, 54.7, 54.8, 54.9, 55. , 55.1, 55.2, 55.3, 55.4, 55.5, 55.6, 55.7, 55.8, 55.9, 56. , 56.1, 56.2, 56.3, 56.4, 56.5, 56.6, 56.7, 56.8, 56.9, 57. , 57.1, 57.2, 57.3, 57.4, 57.5, 57.6, 57.7, 57.8, 57.9, 58. , 58.1, 58.2, 58.3, 58.4, 58.5, 58.6, 58.7, 58.8, 58.9, 59. , 59.1, 59.2, 59.3, 59.4, 59.5, 59.6, 59.7, 59.8, 59.9, 60. , 60.1, 60.2]) - event_time_offset(event_time_offset)float64µs0.0, 249.750, ..., 7.118e+04, 7.143e+04
Values:
array([ 0. , 249.75024975, 499.5004995 , 749.25074925, 999.000999 , 1248.75124875, 1498.5014985 , 1748.25174825, 1998.001998 , 2247.75224775, 2497.5024975 , 2747.25274725, 2997.002997 , 3246.75324675, 3496.5034965 , 3746.25374625, 3996.003996 , 4245.75424575, 4495.5044955 , 4745.25474525, 4995.004995 , 5244.75524476, 5494.50549451, 5744.25574426, 5994.00599401, 6243.75624376, 6493.50649351, 6743.25674326, 6993.00699301, 7242.75724276, 7492.50749251, 7742.25774226, 7992.00799201, 8241.75824176, 8491.50849151, 8741.25874126, 8991.00899101, 9240.75924076, 9490.50949051, 9740.25974026, 9990.00999001, 10239.76023976, 10489.51048951, 10739.26073926, 10989.01098901, 11238.76123876, 11488.51148851, 11738.26173826, 11988.01198801, 12237.76223776, 12487.51248751, 12737.26273726, 12987.01298701, 13236.76323676, 13486.51348651, 13736.26373626, 13986.01398601, 14235.76423576, 14485.51448551, 14735.26473526, 14985.01498501, 15234.76523477, 15484.51548452, 15734.26573427, 15984.01598402, 16233.76623377, 16483.51648352, 16733.26673327, 16983.01698302, 17232.76723277, 17482.51748252, 17732.26773227, 17982.01798202, 18231.76823177, 18481.51848152, 18731.26873127, 18981.01898102, 19230.76923077, 19480.51948052, 19730.26973027, 19980.01998002, 20229.77022977, 20479.52047952, 20729.27072927, 20979.02097902, 21228.77122877, 21478.52147852, 21728.27172827, 21978.02197802, 22227.77222777, 22477.52247752, 22727.27272727, 22977.02297702, 23226.77322677, 23476.52347652, 23726.27372627, 23976.02397602, 24225.77422577, 24475.52447552, 24725.27472527, 24975.02497502, 25224.77522478, 25474.52547453, 25724.27572428, 25974.02597403, 26223.77622378, 26473.52647353, 26723.27672328, 26973.02697303, 27222.77722278, 27472.52747253, 27722.27772228, 27972.02797203, 28221.77822178, 28471.52847153, 28721.27872128, 28971.02897103, 29220.77922078, 29470.52947053, 29720.27972028, 29970.02997003, 30219.78021978, 30469.53046953, 30719.28071928, 30969.03096903, 31218.78121878, 31468.53146853, 31718.28171828, 31968.03196803, 32217.78221778, 32467.53246753, 32717.28271728, 32967.03296703, 33216.78321678, 33466.53346653, 33716.28371628, 33966.03396603, 34215.78421578, 34465.53446553, 34715.28471528, 34965.03496503, 35214.78521479, 35464.53546454, 35714.28571429, 35964.03596404, 36213.78621379, 36463.53646354, 36713.28671329, 36963.03696304, 37212.78721279, 37462.53746254, 37712.28771229, 37962.03796204, 38211.78821179, 38461.53846154, 38711.28871129, 38961.03896104, 39210.78921079, 39460.53946054, 39710.28971029, 39960.03996004, 40209.79020979, 40459.54045954, 40709.29070929, 40959.04095904, 41208.79120879, 41458.54145854, 41708.29170829, 41958.04195804, 42207.79220779, 42457.54245754, 42707.29270729, 42957.04295704, 43206.79320679, 43456.54345654, 43706.29370629, 43956.04395604, 44205.79420579, 44455.54445554, 44705.29470529, 44955.04495504, 45204.7952048 , 45454.54545455, 45704.2957043 , 45954.04595405, 46203.7962038 , 46453.54645355, 46703.2967033 , 46953.04695305, 47202.7972028 , 47452.54745255, 47702.2977023 , 47952.04795205, 48201.7982018 , 48451.54845155, 48701.2987013 , 48951.04895105, 49200.7992008 , 49450.54945055, 49700.2997003 , 49950.04995005, 50199.8001998 , 50449.55044955, 50699.3006993 , 50949.05094905, 51198.8011988 , 51448.55144855, 51698.3016983 , 51948.05194805, 52197.8021978 , 52447.55244755, 52697.3026973 , 52947.05294705, 53196.8031968 , 53446.55344655, 53696.3036963 , 53946.05394605, 54195.8041958 , 54445.55444555, 54695.3046953 , 54945.05494505, 55194.80519481, 55444.55544456, 55694.30569431, 55944.05594406, 56193.80619381, 56443.55644356, 56693.30669331, 56943.05694306, 57192.80719281, 57442.55744256, 57692.30769231, 57942.05794206, 58191.80819181, 58441.55844156, 58691.30869131, 58941.05894106, 59190.80919081, 59440.55944056, 59690.30969031, 59940.05994006, 60189.81018981, 60439.56043956, 60689.31068931, 60939.06093906, 61188.81118881, 61438.56143856, 61688.31168831, 61938.06193806, 62187.81218781, 62437.56243756, 62687.31268731, 62937.06293706, 63186.81318681, 63436.56343656, 63686.31368631, 63936.06393606, 64185.81418581, 64435.56443556, 64685.31468531, 64935.06493506, 65184.81518482, 65434.56543457, 65684.31568432, 65934.06593407, 66183.81618382, 66433.56643357, 66683.31668332, 66933.06693307, 67182.81718282, 67432.56743257, 67682.31768232, 67932.06793207, 68181.81818182, 68431.56843157, 68681.31868132, 68931.06893107, 69180.81918082, 69430.56943057, 69680.31968032, 69930.06993007, 70179.82017982, 70429.57042957, 70679.32067932, 70929.07092907, 71178.82117882, 71428.57142857])
- (distance, event_time_offset)float64Ånan, nan, ..., nan, nanσ = nan, nan, ..., nan, nan
Values:
array([[nan, nan, nan, ..., nan, nan, nan], [nan, nan, nan, ..., nan, nan, nan], [nan, nan, nan, ..., nan, nan, nan], ..., [nan, nan, nan, ..., nan, nan, nan], [nan, nan, nan, ..., nan, nan, nan], [nan, nan, nan, ..., nan, nan, nan]], shape=(105, 287))
Variances (σ²):
array([[nan, nan, nan, ..., nan, nan, nan], [nan, nan, nan, ..., nan, nan, nan], [nan, nan, nan, ..., nan, nan, nan], ..., [nan, nan, nan, ..., nan, nan, nan], [nan, nan, nan, ..., nan, nan, nan], [nan, nan, nan, ..., nan, nan, nan]], shape=(105, 287))
The computed table spans a range of distances and can be plotted using:
[4]:
table.plot()
[4]:
We inspect a horizontal slice of the table around 60m.
The way it was built was to propagate the neutrons in the simulation results read at the position of the last component in the beamline (the t0 chopper) to the specified distance of 60m.
We then compute the event_time_offset of these neutrons, which is simply their arrival time modulo the pulse period of ~71ms.
Finally, we define event_time_offset bin edges (from the TimeResolution parameter) and compute the mean observed wavelength inside each bin.
This is illustrated in the figure below where we plot a map of the neutrons in the (event_time_offset, wavelength) space.
[5]:
sim = wf.compute(unwrap.SimulationResults[AnyRun])
def to_event_time_offset(sim):
# Compute event_time_offset at the detector
eto = (
sim.time_of_arrival
+ ((table.array.coords['distance'][-1] - sim.distance) / sim.speed).to(
unit="us"
)
) % sc.scalar(1e6 / 14.0, unit="us")
return sc.DataArray(
data=sim.weight,
coords={"wavelength": sim.wavelength, "event_time_offset": eto},
)
events = to_event_time_offset(sim.readings["t0"])
hist = events.hist(wavelength=300, event_time_offset=300).plot(norm="log")
fig = pp.tiled(nrows=1, ncols=2)
fig[0, 0] = hist
fig[0, 1] = hist
# Overlay mean on the figure above
table.array["distance", -1].plot(ax=fig[0, 0].ax, color="k", ls="-", marker=None)
table.array["distance", -1].plot(ax=fig[0, 1].ax, color="k", ls="-", marker=None)
xr = 43000, 50000
yr = 2.5, 3.2
fig[0, 0].ax.plot(
[xr[0], xr[1], xr[1], xr[0], xr[0]],
[yr[0], yr[0], yr[1], yr[1], yr[0]],
color='grey',
)
fig[0, 1].canvas.xrange = 43000, 50000
fig[0, 1].canvas.yrange = 2.5, 3.2
fig
[5]:
We can see that the values from lookup table (black line) follow the centerline of the colored areas in the 2d histogram.
Method: ‘analytical’#
This second method creates the lookup table using analytical formulas to propagate and chop a pulse of neutrons through the chopper cascade.
This is fast enough to be performed on-the-fly during data reduction, and should be the default for most reduction workflows.
[6]:
wf = unwrap.GenericUnwrapWorkflow(
run_types=[AnyRun], monitor_types=[], wavelength_from="analytical"
)
wf[unwrap.LtotalRange[AnyRun, snx.NXdetector]] = (
sc.scalar(50.0, unit="m"),
sc.scalar(60.0, unit="m"),
)
wf[Position[snx.NXsource, AnyRun]] = sc.vector([0, 0, -76.55], unit="m")
wf[unwrap.DiskChoppers[AnyRun]] = disk_choppers
wf[unwrap.DistanceResolution] = sc.scalar(0.1, unit="m")
wf[unwrap.TimeResolution] = sc.scalar(250.0, unit='us')
wf[unwrap.PulsePeriod] = 1.0 / sc.scalar(14.0, unit="Hz")
wf[unwrap.PulseStrideOffset] = None
table = wf.compute(unwrap.LookupTable[AnyRun, snx.NXdetector])
table.plot()
[6]:
As above, we inspect more closely the same horizontal slice of the table around 60m.
This time, a single pulse of neutrons represented by a polygon spanning a range of times and wavelengths (blue rectangle) is propagated through the chopper system. It stretches as it travels (slower neutrons arrive later at a given distance) and gets chopped into multiple pieces by the choppers which only let through parts of the original pulse.
In the end, at the detector position, we are left with the pink polygons, which are thin enough to be approximated by a straight line.
[7]:
dist = table.array.coords['distance'][-1]
frames = wf.compute(unwrap.ChopperFrameSequence[AnyRun])
at_detector = frames.propagate_to(dist)
fig, ax = at_detector.draw()
table = wf.compute(unwrap.LookupTable[AnyRun, snx.NXdetector])
# Overlay LUT prediction on the polygons figure
da = table.array["distance", -1]
ax.plot(
da.coords['event_time_offset'].values / 1000,
da.values,
color="k",
ls="-",
marker=None,
)
xr_ms = xr[0] / 1000, xr[1] / 1000
ax.plot(
[xr_ms[0], xr_ms[1], xr_ms[1], xr_ms[0], xr_ms[0]],
[yr[0], yr[0], yr[1], yr[1], yr[0]],
color='grey',
)
[7]:
[<matplotlib.lines.Line2D at 0x7f255cfad390>]
[8]:
ax.set(xlim=xr_ms, ylim=yr)
fig
[8]:
The black line following the center of the polygons is what is stored inside the lookup table.
Before the next section, we save the computed table to disk:
[9]:
table.save_hdf5("wavelength-lut-60m-80m.h5")
Writing type '<class 'NoneType'>' to HDF5 not implemented, skipping.
Method: ‘file’#
The final method is to simply read a pre-computed lookup table from disk.
This is very easy and convenient to use, but one needs to be sure that the chopper settings used in the run to be reduced were the same that were used to compute the lookup table we are loading.
[10]:
wf = unwrap.GenericUnwrapWorkflow(
run_types=[AnyRun], monitor_types=[], wavelength_from="file"
)
wf[unwrap.LookupTableFilename[AnyRun, snx.NXdetector]] = "wavelength-lut-60m-80m.h5"
table = wf.compute(unwrap.LookupTable[AnyRun, snx.NXdetector])
table.plot()
[10]:
