McStas (ii): An overview of components, their use, and advice for user contributions

2.1.Component basics

A McStas component receives an incoming neutron state (a ray) (x,y,z,vx,vy,vz,t,sx,sy,sz,p)_in from the McStas system, then transports, records or manipulates the state and finally returns the resulting neutron state (x,y,z,vx,vy,vz,t,sx,sy,sz,p)_out back to the McStas system.

The component coordinate system can be chosen freely as to whatever is the easy system of the device/problem at hand. However, the usual starting condition in McStas sources is to have z along the neutron beam and y vertical, and since our coordinate system is right-handed x is horizontally transverse to the neutron beam, pointing left as seen from the source. This convention is then naturally carried over to other components in the instrument, e.g. a neutron guide would typically be translated along z seen from the source, as shown in Fig. 2.

Fig. 2.

(top) Illustration of a neutron source (left) and the entrance of the primary guide (right), displaced 1 m along the positive z direction. (bottom) the McStas instrument file code corresponding to the source and the primary guide.

2.2.Component structure

To be a valid McStas component, a file needs to be written in the syntax of the C-language and have the following structure and properties

3.1.Sources

All McStas simulations start with a Source-type component, which can initialize the parameters of the simulated neutron ray, i.e. its position (x,y,z), its velocity (vx,vy,vz), its time (t), and in particular its initial weight, (p0).

In the first McStas review article [29], we presented the rule for intensities of a simulation at component j in the instrument, Ij=∑i=1Npi,j−1, where pi,j−1 is the weight factor of the i’th neutron ray, leaving component number j−1. Hence, if a source has a total integrated brilliance, Φ0, given in units of neutrons per steradian per second, the total emitted intensity in the full solid angle is Iexp=4πΦ0. If we require that all N simulated rays start off with the same weight, we reach Isim=Np0. Demanding that the simulation should represent the experimental reality, Isim=Iexp, we reach

This focusing was also discussed in Ref. [29], where we presented the fundamental equation

This is the initial weight used in the most basic source components, including the component Source_simple(), which is the one used in the first instrument example.

Different continuous sources can be provided in the same way, where only the flux function Ψ(λ) is modified. This is e.g. the case with the component Source_Maxwell_3(), where the flux is given as a sum of three Maxwellians with different temperature and intensity,

For pulsed sources, there is a further complication that the intensity varies with time. The component Moderator() is a simple example of this. Here, the emission time is chosen from an exponential distribution,

A much more elaborate time-of-flight component is ESS_butterfly(), which covers the most recent description of the moderators for the European Spallation Source (ESS). An in-depth treatment of this component is out of the scope of this review, and we refer instead to the online documentation [21].

When a McStas component description or parameter-file does not yet exist for your source, an option is to use the MCPL_input component which makes use of the MCPL [9] framework. MCPL provides a very nice interchange-format for particle lists, with backends available for McStas and McXtrace, SIMRES [18]-[19], Vitess [30], Geant4 [1], PHITS [20] and several variants of MCNP [27]. Thus you can perform your McStas simulations with particle data directly out of your neutron source simulations performed in e.g. MCNP. You should however bear in mind that your derived McStas simulation is in principle statistically limited to what the chosen input-file contains and avioding bias can be dificult. The MCPL_input component provides facilities for smoothing when repeating the particle data an integer number of times, which e.g. happens implicitly when running MPI-parallel simulations. This is done by making Monte Carlo choices within a small sphere of confusion around each event in the input file. When using these user-parameters to resample the MCPL data, please examine your MCPL data carefully and make conservative choices, e.g. based on inspection of the particle distributions. Please also consider the energy-ranges available in your MCPL-input and set relevant, physically meaningful E_min and E_max limits in MCPL_input. – McStas is for instance not suited for transport of eV, keV, or MeV range neutrons!

3.2.Monitors

The term “monitor” in McStas is rather loosely defined and describes any way to gather information about the beam passing a component (typically a specified surface). The neutron ray is propagated until it intersects the specified surface, and the monitor component will determine whether the ray is within the boundary of the monitor. If so, the neutron ray will be registered in a way particular for each monitor (see details below). The monitor will keep track of both number of rays, ∑j1, the sum of the weight factors, ∑jpj, and the sum of weight factor squared, ∑jpj2. These are used to estimate the simulated neutron intensity and its error, as described in Ref. [29].

Unlike a physical monitor or detector, the neutron ray is not affected by the monitor and may continue after the interaction. In fact, if the flag Restore_neutron is set to 1, the neutron ray is even propagated back to its position from before the interaction with the monitor, just to ensure that the monitor is transparent to the other parts of the simulation.

Monitors are as ubiquitous as source components, since this is the only way McStas simulations can produce results. Typically, each monitor will produce one data file for every simulation run. Monitor data are almost always histogrammed. The corresponding bins are either one-dimensional (for example in the L_monitor() and TOF_monitor() described below), or two-dimensional (such as the PSD_monitor() and Divergence_monitor() described below). The physical size, the binning range and the number of bins are taken as input parameters for all monitor components. Routines for formatting one- and two-dimensional data to file are present in the McStas library, to ensure uniformity in data formats.

A number of monitors exist that can perform rudimentary data treatment, for example by transforming time-of-flight and scattering angle into scattering vector, given proper information on the instrument geometry. The use of such components may speed up the design of full instruments, since simulated data are immediately in a human understandable format. However, a detailed description of these more complex monitors is outside the scope of this review.

Below, we describe in more detail some of the most used monitor components.

L_monitor(). This component is a rectangular shaped one-dimensional monitor that registers the true neutron wavelength and histograms the incoming intensity on a number of wavelength bins. A typical output from this monitor is shown in Fig. 3.

TOF_monitor(). This component is very reminescent of L_monitor(), only the neutron rays are distributed according to their arrival time at the detector surface. As such, it can be seen as a simple version of the physical time-of-flight detectors. A typical output from this monitor is shown in a later example, Fig. 9.

Fig. 3.

The simulation results obtained from shooting neutron rays from Source_Maxwell_3() into an L_Mon(), in effect showing the wavelength distribution of the source intensity. The Maxwellian source is run with two non-zero intensities, one corresponding to T=300 K, and one corresponding to T=30 K. The red points show the outcome of the source set with Lmin=1 and Lmax=9; the blue points show the outcome where the two wavelength values are 2 and 6, respectively.

PSD_monitor(). This monitor is another example that actually has a physical counterpart. It produces a two dimensional image of the neutron intensity in real space. It is often used to characterize the beam in a guide or just before the sample – or as the detector for neutron imaging instruments. An example of output from this component is shown in a later example, Fig. 7.

Divergence_monitor(). This two-dimensional monitor ditributes the neutrons according to their divergence in the two dimensions perpendicular to the detector normal. This unphysical monitor is usually used to inspect the phase space distribution of the beam inside a neutron guide system or at the sample position.

Monitor_nD(). This monitor can measure almost anything you can imagine, by means of a flexible options parameter string. Monitor_nD is very powerful, but using it correctly can take some effort and describing all features is out of the scope of this review. Instead, we refer to the relevant section of the McStas component manual, the Monitor_nD mcdoc page and the many McStas example instruments that make use of this component. The monitor can assume many geometrical shapes including flat panels, cylinders, spheres and a surface-description in OFF format, where each surface polygon becomes a pixel on the detector. The monitor generates 1-dimensional histograms (or multiples of these), 2-dimensional histograms and event lists with user-selectable columns of neutron data. This means that Monitor_nD can be used to measure anything the other monitors can do, and more. One powerful speciality is the use of the special user1,2,3 variables which can measure any c variable in your instrument file. An example of this can be found in the ILL_IN6 instrument by Emmanuel Farhi, found in the McStas examples: A variable measures which of the 7 IN6 monochromator blades has scattered the neutron, to be later correlated against the beam wavelength at the sample position. Further McStas-Mantid interface [17] to the Mantid [2] data reduction framework is implemented via a combination of mcdisplay.pl and Monitor_nD. Lastly, Monitor_nD can measure neutron states elsewhere upstream in your instrument file using the previous keyword or in combination with PreMonitor_nD (see below).

PreMonitor_nD(). The PreMonitor_nD component can store the neutron properties for use later in the simulation work flow, e.g. measure the originating state of the neutron at the source to be studied at the sample position. The component has a single input, namely the instance name of the Monitor_nD which will perform the measurement. (That Monitor_nD should correspondingly have premonitor in its options string.) The dedicated example instrument Test_PreMonitor_nD available in the McStas installation demonstrates its use.

3.3.Optics

We now turn to the optical components that in general are used to modify the beam characteristics. The geometry of the components in this section is shown in Fig. 5.

Arm(). This is possibly the most used McStas component, yet it does not affect the neutron ray at all(!) The component is to be used as a reference for the internal coordinate system in McStas. Typically, an Arm() is placed at an axis of rotation, and following components – including further Arm()’s – are to be placed relative to the first Arm(). Used in this way, the Arm() component has the same function as an optical bench.

Slit() aka. Diaphragm(). This component is a simple slit, which in effect is an infinitely wide, 100% effective absorber with a single central opening, being either circular or rectangular. As rectangular, the slit opening width is given by the parameters xmin and xmax. Alternatively, a symmetric opening is given by instead setting the parameter xwidth. The corresponding parameters setting the opening height are denoted ymin and ymax, alternatively yheight. These variable names are generally used for heights and widths of components within McStas – as described in the (rarely studied) NOMENCLATURE document on the McStas home page [21].

If Slit() should rather have a circular opening (in this case centered at the origin of the local coordinate system), this is done by setting the input parameter radius. This way of controlling the component functionality via setting of input variables is a frequent feature within McStas library components.

The Slit() functions in the following way: First the neutron ray is propagated to the plane of the component. Next, a check is performed to see if the neutron has hit the slit opening. If not, the neutron ray is discarded, using the McStas command ABSORB.

Collimator_linear(). This component models a linear Soller collimator, physically made from a stack of absorbing, parallel blades. It is box-shaped with the physical dimensions xwidth, yheight, and length, as shown in Fig. 5. The degree of collimation is parametrised as the maximal divergence allowed through, using the parameters for the two directions: eta_x and eta_y – just as for a Soller collimators on an experimental set-up. If any of these parameters is zero, this is interpreted as if there is no collimation restriction in the corresponding direction (the default value).

The neutron ray is first propagated to the entry plane, then to the exit plane of the component. If the ray does not pass the collimator opening at both planes, it is ABSORB’ed. If the neutron ray has a larger divergence than the maximally allowed, in either the x or the y direction, it is also ABSORB’ed.

In a physical Soller collimator, the transmission of a certain neutron is either 0 or 1, depending whether the neutron trajectory intersects one of the absorbing plates. The test for collision for each of the plates is somewhat tedious, so in this McStas component we use a stochastic approach. If we assume that we do not know exactly where the absorber blades are positioned (but only the distance between them), we can calculate the probability of being transmitted as a function of divergence (for collimation in one dimension):

Guide(). McStas contains a number of neutron guide components. We here describe the simplest of them, Guide(). This is a neutron guide with rectangular openings. The entry has the dimensions h1 × w1, while the exit has dimensions h2 × w2, i.e. it is possible to describe a linearly tapering guide by making the entry and exit rectangles having different sizes. The length of the guide is given by l.

The inner surface of Guide() is imagined to be covered by a neutron-reflecting layer or layers, where the outermost layer is Ni, which has a critical scattering vector of qc,Ni=0.0219 Å⁻¹. This means that neutron rays with a sufficiently low incident angle, θ to the surface can be specularly reflected with probability R=1.

To quantify this, we need to define the actual scattering vector, q, in the conventional way, q=2ksin(θ). To model a supermirror guide, we define the parameter m, representing that the extended critical scattering vector is qc=mqc,Ni. However, the probability of reflection decreases for q>qc,Ni. This is parametrised by the value α, which modifies the reflectivity (for small q values) as

The component functions by propagating the neutron ray to the entry plane. If the rectangular opening is not hit, the neutron is ABSORB’ed. Otherwise it continues into the guide. Each time the ray intersects one of the four guide sides, a specular reflectivity event is assumed. The scattering vector, q, is calculated, and the neutron weight multiplier becomes

For a consistent description of the mirror reflectivity, it can be recommended in stead of R0 and alpha to use the library ref-lib, which contains updated descriptions of reflectivity for mirrors of different values of m. This library is used for all McStas components that decribe optics with reflecting surfaces, such as Guide_gravity(), Guide_wavy(), Mirror() and several others. The current version of the reflectivity profiles is illustrated in Fig. 4.

Fig. 4.

The reflectivity profiles, R(q), used by mirror and guide components alike. The figure shows plots for m=1, m=2, and two SwissNeutronics m=4 models derived by Klaus Lieutenant and Henrik Jacobsen, respectively. The Klaus Lieutenant models (2020) are distributed as text files in the McStas data folder, the Henrik Jacobsen (2013) models are easily selected by setting α=0 and W=0 for any m>3 mirror.

Monochromator_flat(). This component models a simple crystal monochromator of dimensions yheight × zwidth. The coordinate z is here used to denote the width, since the unrotated monochromator lies in the (y,z)-plane. This ensures that the ROTATE angle of the monochromator becomes equal to the angle θ in the Bragg law

The Monochromator_flat() is assumed to be an ideally imperfect mosaic crystal, i.e. it is assumed to have a Gaussian mosaic spread of tiny crystallites, each too small to have significant primary extinction. This, in turn, gives rise to a Gaussian distribution of the scattering from the monochromator. The divergence of the scattered beam is characterized by the parameters mosaich and mosaicv, for the horizontal and vertical directions, respectively. In other words, mosaich is the effective mosaic in the horisontal direction, such that an in-plane oriented monochromator will display a rocking curve width (FWHM) with the value mosaich. Likewise, mosaicv is the effective mosaic in the vertical direction, such that a pencil-beam reflected from the same monochromator will have a vertical spread (FWHM) of 2sin(θm) mosaicv.

Two other monochromators exist in the McStas package: Monochromator_curved() is a multi-slab composite of crystals, focusing the scattering to a particular point, while Monochromator_pol is a flat, neutron polarizing monochromator. Polarized neutrons are described in a later article in this review series.

Monochromators of perfect (or bent perfect) single crystals are at present not present in McStas, as primary extinction is not implemented. An approximation does however exist in the Single_crystal component, where the crystal lattice can be bent locally, keeping the external geometry of the crystal unchanged. In this approximation, curvature is spherical along the outer vertical and horizontal axes of the crystal.

DiskChopper(). This component models a disk chopper; a neutron-absorbing flywheel spinning with its axis along the beam direction, and with wedges cut into the wheel to let through neutrons at particular times. A sketch of the component is shown in Fig. 5.

The component takes the geometry parameters radius, nslit, yheight, theta0 – the latter representing the angular width of each slit (in degres). In addition, we have the chopper frequency nu (in Hz) and the phase of the choppers (at time zero), phase.

A special feature of this component is its effect on a continuous beam. Here, the intention is that the beam will be pulsed after hitting the chopper. However, naively used, this could result in a huge loss in simulated neutron rays. This we have remedied by letting all neutron rays pass the first chopper, while its time is being adjusted to a (random) value where the chopper was indeed open for this ray. This functionality is activated by setting the flag isfirst. It should, however, be used on the fist chopper only, since otherwise very misleading results could arise.

McStas also features other rotating optical components, such as the descriptions of spinning velocity selectors, Selector() and V_selector(), and two different descriptions of Fermi choppers, FermiChopper() and Vitess_ChopperFermi(). For descriptions of these components, we refer to the McStas documentation.

Fig. 5.

Geometries of the five described beam-optical elements with the standard naming of geometrical parameters shown. Top: Guide(), Monochromator(), Collimator() and Slit(). Bottom: DiskChopper().

3.4.Samples

At this place, we only briefly touch upon scattering samples in McStas, since a later article in this series is dedicated to exactly samples. For this reason, we here initially present two very simplified examples and then only list all presently available samples in McStas.

Incoherent(). This sample scatters all incoming neutron rays incoherently, i.e., uniformly in the 4π solid angle (or by focusing in a smaller solid angle, ΔΩ). In addition, Incoherent() can scatter inelastically with a fixed energy transfer (for calibration purposes). The sample can take many different shapes.

Powder1(). This sample introduces powder diffraction from a single value of q. The neutron rays are Bragg scattered out into different directions, forming a Debye-Scherrer cone. The scattering and absoption cross sections for the sample are parametrized by F2 and sigma_abs, respectively. The sample has a larger sibling, PowderN(), which can tackle an arbitrary number of powder lines, read from an specified crystallographic input file.

The McStas sample suite. For completeness, in Table 1 we list the presently available scatering samples in McStas. All samples take as parameters their physical size, their absorption (and possibly incoherent) cross section, and then a parametrized or file-based description of their scattering cross section.

3.5.Misc

Our miscellaneous component category contains the components that do not naturally fit into the other categories where physical properties are a common denominator, e.g. sources or samples. The most notable misc components are these:

3.6.Contrib

The contrib component category is where a user-contributed component is placed for distribution with McStas. In principle we offer the same technical support for contributed components as the official components of the other categories, but the McStas authors did not write these component themselves and can therefore naturally not take full responsibility for their behaviour or performance. Yet, the contributed components are generally speaking both relevant and working as they should. There are currently 117 components in the contrib section, of different types. For author information, please consult the individual component definition.

4.1.A two-axis diffractometer on a continuous source

Fig. 6.

The McStas instrument file describing a simple two-axis diffractometer on a continuous source.

As the first example, we have assembled a very simple two-axis diffractometer on a continuous source, see Fig. 6. The instrument consists of a Source_simple() that feeds a neutron beam into a 20 m long, 12×3 cm² cross section, m=2 Guide(). Just after the guide is placed a 15×20 cm² Monochromator_flat() with parameters relevant for pyrolytic graphite. The monochromator is rotated with respect to the main beam by the angle A1. Next is placed an Arm() rotated by A2 with respect to the main axis, and on the arm is placed an L_monitor() and then a cylinder-shaped (3 cm tall, 1 cm radius) Incoherent() as the scattering sample. Incoherent() is here used in its default configuration where all neutron are scattered, but where the scattering weight is given by the cross section of metallic vanadium. (A detailed description of sample scattering will appear as a later article in this review series.) After the sample, another Arm() is placed, turned at an angle A4 with respect to the first Arm(), and as the detector, a 20×20 cm² PSD_monitor is placed on this second Arm().

Typical results of simulations with this test instrument are shown in Fig. 7. The simulation is run with its default parameters: A2 value of 74.17° (corresponding to scattering of 5 meV neutrons on a PG monochromator), A1 being the half of this. Finally, we use as scattering angle, A4, a (somewhat arbitrary) value of 60°. The L_monitor shows a clear peak at λ=4.04 Å, corresponding to 5 meV as presumed. The PSD_monitor shows a uniformly illuminated area of 10×10 cm² with zero counts outside this. This limited illumination is an effect of the focusing applied in Incoherent(), where the focusing area is 10×10 cm² at target_index 2, that is 2 components downstream from Incoherent(), which is indeed the PSD_monitor. Hence, this example serves as well to illustrate the effect of focusing in McStas.

Fig. 7.

The simulation results obtained from the simple two-axis diffractometer in Fig. 6. (left) the wavelength distribution on the sample position as measured by L_monitor(). (right) The beam profile on the detector, measured by PSD_monitor().

4.2.A powder diffractometer on a pulsed source

The second example illustrates a very simple powder diffractometer at a pulsed source, as shown in Fig. 8. This instrument starts with a Progress_bar() that does nothing but generate information on the progress of the simulation. The neutron source is Moderator(), which is set to emit neutrons with energies between 2.5 meV and 6 meV. This is pulsed source, set to a pulse length of 100 μs. the neutron rays emitted from Moderator() are focused on a 3×3 cm² square located at the position of the 3rd component downstream (target_index = 3), which is the Slit() described below. The flux of the moderator is given in units of neutrons per second per square centimeter per steradian per meV energy.

After the source is placed a 50 Hz DiskChopper(), whose phase and opening angle can be controlled by external parameters, in this case PHASE and THETA. Next, the divergence of the beam is controlled by a Collimator_linear(), set to a maximum divergence of 60′ (60 arc minutes =1°). The beam is further restricted in size by a Slit() of size 3×3 cm². The sample is the simple powder scatterer Powder1(), a 3 cm tall cylinder with radius of 8 mm, which is set to scatter all neutron rays at the first Bragg peak of Al₂O₃, at q=1.8049 Å⁻¹. the weight of the neutron is adjusted according to the default cross section, which is exactly the one for this Al₂O₃ peak. The function of samples will be detailed further in a coming article in this review series. Finally, an Arm() is set to the fixed scattering angle of 90°, and on the arm a 25 mm wide TOF_monitor() is used as a detector, placed 2 m from the sample.

Figure 9 shows the outcome of this simple virtual experiment. A peak in intensity is seen close to a time of 28.7 ms. With a total flight path of L=23 m, this corresponds to a wavelength of λ=t/(αL)=4.93 Å, with α=252.7μs/m/Å. Using the Bragg law, Eq. (17), for a scattering angle θ=45°, we reach dexp=3.486 Å- quite in agreement with the assumed value dinput=2π/1.8049=3.481 Å. A more detailed simulation and data analysis, taking into account the time structure of the moderator, would be able to give a higher precision on the estimated lattice spacing. We encourage the readers to test this out for themselves.

Fig. 8.

The McStas code for a simple powder diffractometer on a pulsed source.

Fig. 9.

The simulation results obtained from the simple pulsed-neutron powder diffractometer in Fig. 8: the time dependence of the Braggdiffracted pulse on the detector position, as seen by a TOF_monitor().