Skip to content

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')
(1500, 2000) time by inverse velocity
9814 of 3000000 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')
grid    4.905000e-03 from 9814 cells
region  3.929119e-03 exactly
grid is 24.84% 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. Region has 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.