Skip to content

Weight the region by a source

Region says which neutrons a train passes, not how many: whether a neutron gets through depends only on where it is in (inverse velocity, emission time), not on how intense the source is there. So the transmitted intensity is the source intensity integrated over the region, and the region only has to be built once.

A gridded source

Give the source as a histogram: the amount in each (time, inverse velocity) cell, spread evenly over it — what a Monte Carlo source or a tabulated moderator spectrum gives you. Shape it as the masks are, (time, inverse velocity), with bin edges in s and s/m.

# bin edges over the source rectangle: inverse velocity in s/m, emission time in s
corners = [u for u, _ in source.polygons[0].vertices]
inverse_velocities = np.linspace(min(corners), max(corners), 401)
times = np.linspace(0.0, 3e-3, 61)

# the amount in each (time, inverse velocity) cell -- here a stand-in that peaks early
centres = 0.5 * (times[1:] + times[:-1])
signal = np.exp(-centres / 1e-3)[:, None] * np.ones(len(inverse_velocities) - 1)

coverage = region.coverage(inverse_velocities, times)
print(coverage.shape, f'{coverage.max():.3f}')
print(f'transmitted {region.integrate(signal, inverse_velocities, times):.4f} '
      f'of {signal.sum():.1f}')

(60, 400) 1.000
transmitted 28.1285 of 7600.9

coverage is the fraction of each cell the region covers, found by clipping each polygon to each cell it reaches. integrate is the signal times the coverage, summed. Both are exact for a histogram, where inverse_velocity_time_mask and unmasked_probability count a partly covered cell whole.

At a monitor or the sample

A neutron emitted at inverse velocity a and time t reaches distance L at t + L*a. That shear keeps areas, so arrivals carries each cell's share of the region into arrival-time bins exactly: the density times the area between two lines of constant arrival time.

frame = 1 / 14
arrival_times = np.linspace(0.0, frame, 715)        # 0.1 ms bins over one frame
at_sample = region.arrivals(signal, inverse_velocities, times, 162.0, arrival_times,
                            period=frame)
profile = at_sample.sum(axis=1)                     # over inverse velocity
peak = arrival_times[np.argmax(profile)]
print(at_sample.shape, f'peak at {peak * 1e3:.1f} ms')

(714, 400) peak at 51.8 ms

The result is (arrival time, inverse velocity), because the shear leaves inverse velocity alone: sum the last axis for the time profile, or keep it for the wavelength as well. period folds arrival times into one frame, as a detector's event_time_offset does; without it, give arrival bins that reach the slowest neutron. Distance 0 gives the emission times back.

The mean and spread in each arrival bin

A wavelength lookup table stores, for each arrival time at each distance, the mean wavelength of what arrives and its spread. arrival_moments gives exactly that, one value per arrival bin:

unfolded = np.linspace(0.10, 0.18, 801)               # 0.1 ms bins, no folding
moments = region.arrival_moments(signal, inverse_velocities, times, 162.0, unfolded)
busiest = np.nanargmax(moments.amount)
print(f'busiest bin {unfolded[busiest] * 1e3:.1f} ms: '
      f'{moments.wavelength[busiest]:.4f} ± {moments.sigma[busiest]:.4f} Å')

busiest bin 123.2 ms: 3.0061 ± 0.0023 Å

Each piece of the region contributes its own first and second moments, integrated over the clipped polygon, rather than its cell's centre. So the mean and spread are exact for the gridded source on any grid, and a coarse grid loses only how well it describes the source. Without choppers, at 28.7 m from an ESS source, binning the source at 0.2 Å by 100 µs instead of tof's own 0.02 Å by 6 µs moves the mean by at most 0.002 Å.

Every neutron is given one flight path. A guide's spread of paths smears arrival times, and is not included here.

A finite beam

With disc apertures, which chopcal.bifrost sets by default, a Region's boundary is where transmission falls to zero: everything inside counts as passing, including neutrons that only the edge of the beam carries through. Inside, an edge sweeping a uniformly lit beam ramps the transmission linearly, and transmit_weighted carries that ramp:

from chopcal import constants as c
from chopcal.lib import Chopper, beam_aperture

# the ramp is set by the angle the beam subtends at its centre, not at its inner corner
narrow = beam_aperture(c.DISK_RADIUS, c.NARROW_SLIT_HEIGHT, c.NARROW_WINDOW_WIDTH,
                       c.NARROW_WINDOW_HEIGHT, at='centre')
wide = beam_aperture(c.DISK_RADIUS, c.BANDWIDTH_SLIT_HEIGHT, c.BANDWIDTH_WINDOW_WIDTH,
                     c.BANDWIDTH_WINDOW_HEIGHT, at='centre')
ramped = [Chopper(d.speed, d.delay, d.beam, list(d.edges), d.path,
                  wide if name.startswith('bw') else narrow)
          for name, d in chopcal.bifrost(wavelength_max=3.0).items()]
# the pulse-shaping pair and the bandwidth pair each share a housing
weighted = source.transmit_weighted(ramped, groups=[0, 0, 1, 2, 3, 3])
point = source.transmit(list(chopcal.bifrost(wavelength_max=3.0, apertures=False).values()))
print(f'support {region.area:.4e}, point beam {point.area:.4e}, '
      f'weighted {weighted.integral:.4e} s^2/m')
print(f'transmitted {weighted.integrate(signal, inverse_velocities, times):.4f}')
softer = weighted.arrivals(signal, inverse_velocities, times, 162.0, arrival_times,
                           period=frame)

support 1.3289e-07, point beam 9.2192e-08, weighted 8.9711e-08 s^2/m
transmitted 17.1654

The result is a WeightedRegion, with the same coverage, integrate and arrivals; each piece counts by the integral of its transmission rather than its area. Here the source passes 17.2 rather than the 28.1 the support counts.

The choppers handed to transmit_weighted carry beam_aperture(..., at='centre'). The region's edge is where the last of the beam is cleared, at the window's inner corner, which is what chopcal.bifrost and the default at='inner' give. The ramp in between is set by the beam centre: while an edge spans the window's height, the fraction it has passed grows linearly at the rate the centre radius sets. Against McStas with the real window geometry, the centre angle fits a BIFROST train's arrival times to within counting error; the inner one leaves the tails of the profile 6-18% low.

groups says which discs share a housing. A neutron crosses those at the same place across the beam, so a group passes the part of the beam open in all of them. For a pair turning the same way that is the smaller of the two transmissions. For a counter-rotating pair, like BIFROST's bandwidth discs, equal openings sweeping from opposite sides stay shut until the edges meet at the half-way time and then open twice as fast. Different groups multiply, as reflections in a guide scramble where across the beam a neutron is. Without groups every disc is its own group.

The point beam, apertures=False, has its boundary where each ramp is half open. Its area is the first-order answer, within 2.8 % of the weighted one here; the difference is where ramps of different discs overlap.