Chopcal

Description

Exposes chopper calculations to Python which were otherwise hidden in McStas instruments.

Supported components

Placing the band

BIFROST passes a fixed bandwidth of about 1.91 Å, so one number decides where that band sits. Give it as either the longest wavelength or the lowest energy you want to reach the sample — the band runs from there to about 1.77 Å shorter.

>>> import chopcal
>>> chopcal.bifrost(wavelength_max=3.0)
name  speed [Hz]  delay [ms]  beam [deg]  opening [deg]  open [ms]  aperture [deg]  path [m]
 ps1         196      3.8781           0            170      2.409           5.651     6.342
 ps2         196      6.0874           0            170      2.409           5.651     6.362
 fo1          14     6.13997           0          38.26      7.591           5.651      8.53
 fo2          14     9.54764           0          52.01      10.32           5.651    14.973
 bw1          14     42.8823           0            161      31.94           13.38        78
 bw2         -14     42.8823           0            161      31.94           13.38     78.02
>>> chopcal.lib.wavelength_limits(list(chopcal.bifrost(wavelength_max=3.0).values()))
(1, (1.2126790814161053, 2.984239670498114))

Mind the sense of it: the number you give is the slow end of the band. Asking for energy_min=7 does not cap the incident energy at 7 meV — it puts 7 meV at the bottom and lets everything faster through, up to about 31 meV.

The choppers come back by name in beam order. Each of BIFROST's is a single opening centred on the beam, so beam is 0 and edges is [-width/2, +width/2]; the table's open column is how long the beam spends inside the openings each turn, and delay is when the disk point at beam is on the beam path. A chopper therefore passes neutrons for open seconds around delay, and again every 1/speed after that — a little longer than that, in fact, because aperture is how wide the beam is on the disk and every window opens half of it early and closes half of it late.

Note the second window in the output above. It is not a band BIFROST delivers; see The band list, and a caveat.

Installation

pip install chopcal

Or, from the source repository to get the latest development version

pip install git+https://github.com/g5t/chopper-calculations.git

Describing a chopper

A Chopper is speed, delay, beam, edges, path and aperture — the disk described the way the NeXus NXdisk_chopper standard and McStas' CollectorDiskChopper component describe one, so the disk's own numbers go in unchanged.

edges are the slit edges in degrees from the disk's zero mark, two per opening, strictly increasing, opening edge first; beam is the angle from that mark to where the beam crosses the disk. An edge at angle a is on the beam at

t(a) = delay + (beam - a) / (360 * speed)

and every 1/|speed| seconds after that. Only beam - a appears, so there is no conversion to do. A slit straddling the mark is written past 360 — [350, 370] rather than [350, 10] — so the pairs stay ordered and each width is a difference. Negative angles are fine: one opening astride the mark reads better as [-85, 85].

Chopper.centred(speed, delay, width, path) writes a single opening centred on the beam, which is what the angle field of chopper-lib 2.x meant and all chopcal.bifrost produces.

aperture is the angular width of the beam where it crosses the disk, which widens every window by half of it at each end; 0, the default, is the point beam every other field describes. It is an angle rather than a width in metres because that is what the disk sees — see chopper-lib's header for the conversion, which is not a plain division.

delay is a time, so it means the same thing at any speed and in either direction of rotation — which is how a real chopper is set, and what McStas' DiskChopper acts on. It replaces the phase in degrees that chopper-lib took before its version 2.0.0, and which it converted to a delay at every point of use anyway.

A disk that is not turning is open or shut for good, and the two mean opposite things: the window and mask functions leave one parked open out of the calculation, having no period to constrain anything, while one parked shut is a beam stop and empties the answer. Chopper.parked_is_open tells them apart, which an empty result cannot.

Masking and sampling

inverse_velocity_time_mask returns which (inverse velocity, time) bins the train passes, as an array shaped (time, inverse velocity) — time is the slow axis, the opposite of the argument order. It takes bin edges, so the mask is one smaller in each direction.

For BIFROST that mask allows well under a percent of the plane, so a source that draws uniformly and discards what the mask excludes spends almost its whole ray budget. A MaskSampler draws from the allowed cells instead, never rejects, and carries the acceptance every drawn ray's weight must be multiplied by to keep the estimate unbiased. Over the same 200 × 300 grid that is a ~300× smaller error for the same number of rays.

acceptance is a ratio of areas, so it is exact — but only while both coordinates are drawn uniformly and independently of each other and of everything else the caller samples. A source that picks its emission time from a window centred on the neutron's own velocity, as McStas' ESS_butterfly does under time focusing, breaks that, and no single number corrects it. It is also not unmasked_probability, which is the allowed fraction of a weighted signal — the transmission an instrument sees, and the wrong normalisation for this.

These three need numpy, and only when called; nothing else in chopcal does.

The band list, and a caveat

wavelength_windows gives the bands the train passes; wavelength_limits gives the envelope around them. Reach for the first when there is more than one, because an envelope spanning two bands also spans the gap between them.

Both are an over-approximation. chopper-lib works out each disk's admissible inverse velocities letting the emission time range over the whole pulse independently, then intersects those ranges across disks — and an intersection of projections is a superset of the projection of the intersection. It can therefore report a band that no single emission time actually delivers.

It is tight wherever the disks leave wide, overlapping emission windows, which is the usual case and includes every band an instrument is set up to pass. BIFROST with its real beam is the exception: a ~0.03 Å sliver near 38 Å, where the frame overlap disks each leave only tens of microseconds and those windows do not overlap. Three independent checks agree nothing gets through — a brute-force scan over emission time, an exact interval intersection, and inverse_velocity_time_mask, which keeps the two coordinates coupled and cannot make the error.

So: where the mask and the windows disagree, the mask is right — and Region is righter still, being neither projected nor quantised. test_mask.py records this one.

The transmitted region, exactly

wavelength_windows over-reports and inverse_velocity_time_mask quantises. Since chopper-lib 4.2.0 the region a train actually transmits can be built instead, as a union of convex polygons in (inverse velocity, emission time):

>>> import chopcal
>>> from chopcal.lib import Region
>>> train = list(chopcal.bifrost(wavelength_max=3.0).values())
>>> source = Region.from_wavelengths(0.4, 45.0, 0.0, 3e-3)   # AA, AA, s, s
>>> region = source.transmit(train)
>>> len(region), [len(p) for p in region.polygons]
(1, [5])
>>> region.wavelength_ranges()
[(1.145836, 3.051543)]
>>> region.area / source.area          # the acceptance, exactly
0.003929118918304158

A neutron emitted at inverse velocity a and time t reaches path L at t + L*a, so a disc open on [lower, upper] accepts lower <= t + L*a <= upper — a slab between two parallel lines. A train's acceptance is an intersection of unions of such slabs, and intersection distributes over union, so the region is a union of convex pieces: one per choice of which opening and which turn of each disc a neutron goes through. Nothing is approximated, and for the six BIFROST discs it comes out as a single polygon of five vertices.

transmit returns a new Region and leaves the one it was asked of alone, so a source rectangle can be reused.

The same train through wavelength_windows reports that band and a 0.03 Å sliver near 38 Å that nothing passes — see the caveat above. The region does not, because it never separates the two coordinates.

A guide that is not a straight line

transmit takes one extra flight path per chopper, in metres, for a guide in which a neutron travels further than the straight line. The deviation is in path, so its effect on an arrival time is deviation * inverse_velocity — larger for a slow neutron, and nothing at all to the inverse velocity — so it does not grow the region evenly; it opens each slab into a wedge:

>>> spread = [1e-4 * chopper.path for chopper in train]
>>> source.transmit(train, spread).wavelength_ranges()
[(1.145722, 3.051543)]

Only the slow edge moves: a longer route rescues a neutron that would have arrived early and can do nothing for one arriving late. Like Chopper.aperture it gives the support under that uncertainty rather than a distribution over it. A spread wide enough to reach from one turn of a disc into the next would make the transmitted pieces overlap and their areas count twice; that raises rather than returning a number.

Sampling it

>>> sampler = region.sampler(source.area)
>>> sampler.count, sampler.acceptance
(3, 0.0039291189183041555)
>>> inverse_velocity, time = sampler.sample(100_000)    # needs numpy

The same job as MaskSampler and the same shape — three uniform deviates, one binary search, never rejects — over the exact region. acceptance here is exact: both areas are known in closed form rather than counted in cells, where the grid can only over-estimate, since a partly covered cell is weighted whole.

Writing it out

region.to_dict(sampled=source) gives plain data; region.write_json(path, sampled=source) writes the same thing to a file, byte for byte what the Polygon_ESS_butterfly McStas component writes beside its data. Every number carries enough digits to read back bit-exact, and the file is a few hundred bytes where a mask of any useful resolution is megabytes.

Describing a disk's beam

Chopper.aperture is the angular width of the beam where it crosses the disk, and beam_aperture computes it from the geometry:

>>> from chopcal.lib import beam_aperture
>>> beam_aperture(radius=0.35, slit_height=0.09846, window_width=0.060, window_height=0.090)
13.379...

slit_height is the radial extent of the opening — McStas DiskChopper's yheight, which puts the beam centre at radius - yheight/2 — and the window is the beam inside it. The largest angle belongs to the window's inner corners, nearest the spindle, where a given width subtends the most angle. Taking the width over the beam-crossing radius instead misses the height entirely: 11.4° where the real figure is 13.4°.

chopcal.bifrost fills this in from the instrument's own disk geometry. Pass apertures=False for the pencil beam every release before this one described, and the narrower band that went with it — BIFROST_BANDWIDTH_POINT_BEAM rather than BIFROST_BANDWIDTH.

Which of the three to reach for

Region is exact and the other two are not, so it is the right default for anything that feeds a number into a calculation. The window functions stay the cheapest way to ask whether a train passes roughly the band you meant, and the mask stays the one that gives you an array to plot.

Developing

The build fetches chopper-lib from GitHub. To build against a local checkout instead — which you need when changing both together — point FetchContent at it:

pip install --no-build-isolation -ve . \
  -Ccmake.define.FETCHCONTENT_SOURCE_DIR_CHOPPER_LIB=/path/to/mcstas-chopper-lib
pytest test

Nothing else changes: the same sources are compiled into the extension modules, so a modification to chopper-lib.c shows up in chopcal on the next build. chopper-lib carries a CHOPPER_LIB_VERSION macro and src/chopper.h asserts on it, so building against a version that means something else by a chopper's fields fails rather than returning different numbers. chopcal needs 4.2.1 or newer — 4.2.0 has the polygons but merges overlapping ranges in a way that depends on the platform's qsort, so wavelength_windows on a real train answered differently on Windows than on Linux.

Only chopper-lib.c is compiled with CHOPPER_LIB_DEFINITIONS, not the whole target: those are V2K, K2V and PI as macros, and chopcal::constants declares C++ variables by the same names.