Work on a grid instead¶
Sometimes you want an array to plot rather than a region to reason about. The mask is that array — and it is the one part of chopcal that needs numpy.
Build a mask¶
inverse_velocities = np.linspace(wavelength_to_inverse_velocity(0.4),
wavelength_to_inverse_velocity(45.0), 2001)
times = np.linspace(0.0, 3e-3, 1501)
mask, allowed = inverse_velocity_time_mask(train, inverse_velocities, times)
print(mask.shape, 'time by inverse velocity') # note the order
print(f'{allowed} of {mask.size} bins allowed')
The shape is (time, inverse velocity)
Time is the slow axis — the reverse of the argument order. That is how the C stores it, and it is easy to get backwards.
inverse_velocities and times are bin edges, so the mask is one smaller in each
direction. Pass lists if that is what you have; they are converted.
grow=n expands each allowed region by n bins in every direction. It is a blunt
instrument — it admits inverse velocities no disc ever passes — so prefer
Chopper.aperture, which opens the windows in time
only, where the beam's width actually acts.
Transmission through it¶
from chopcal.lib import unmasked_probability
signal = np.ones(mask.shape)
print(f'flat signal transmits {unmasked_probability(signal, mask):.6e}')
unmasked_probability is the allowed fraction of a weighted signal — the transmission
an instrument sees. It is not the weight correction a sampler needs; see
Sample without wasting rays.
What the grid costs¶
from chopcal.lib import MaskSampler
grid = MaskSampler(mask, inverse_velocities, times,
inverse_velocities[0],
inverse_velocities[-1] - inverse_velocities[0],
times[0], times[-1] - times[0])
source = Region.from_wavelengths(0.4, 45.0, 0.0, 3e-3)
exact = source.transmit(train).area / source.area
print(f'grid {grid.acceptance:.6e} from {grid.count} cells')
print(f'region {exact:.6e} exactly')
print(f'grid is {grid.acceptance / exact - 1:.2%} high')
The error is one-signed — a partly covered cell is counted whole, so a grid can only over-accept — and it falls only as fast as the cell count, so halving it costs four times the memory. Three million cells is still a quarter high here, on a deliberately wide source; on a tight one it does much better, but it is never exact.
That matters most for MaskSampler.acceptance, which multiplies every drawn ray's weight.
A grid figure that is 25% high makes every weight 25% wrong.
When the grid is still the right tool¶
- You want a picture.
Regionhas no pixels. - You are reading a mask written by a McStas component.
- The over-estimate is harmless and you would rather not carry a polygon set around.
Otherwise reach for Region.