McStas (i): Introduction,use,and basic principles for ray-tracing simulations

Abstract

We present an overview of, and an introduction to, the general-purpose neutron simulation package McStas. We present the basic principles behind Monte Carlo ray-tracing simulations of neutrons performed in the package and present a few simple examples. We present the implementation of McStas, the status of the package and its use in the neutron community. Finally, we briefly discuss the planned development of the package.

1. Introduction

Over the last decades, Monte Carlo ray-tracing simulation has become increasingly important for investigating the performance of neutron scattering instrumentation. Few neutron instruments are being built before undergoing a thorough investigation by simulation. In particular, the powerful spallation neutron sources in North America (SNS) [54] and Japan (J-PARC MLF) [49], the upgrade of the UK source (ISIS TS-2) [48], and the coming European Spallation Source (ESS) [47] have heavily utilized neutron simulation software, thereby also pushing the development of the simulation methods.

At the end of last century, neutron instrumentation ray-tracing simulations were performed using home-made monolithic codes. While fairly successful, this approach suffered from limited development ressources – with potential problems of quality assurance – while still creating a vast amount of duplicated effort between different projects and facilities. One exception from the monolithic approach was the general-purpose neutron simulation package NISP [7,52,57] released 1994/95, based on the earlier MCLIB subroutine library from 1978/1980 [17,18]. However, NISP never reached a critical mass of user community worldwide.

Almost simultanously, at the turn of the century, four general simulation packages were launched, following the NISP philosophy, but on more modern software platforms:

The HZB-based VITESS [68]

The Řež-based RESTRAX/SIMRES [38,39]

The SNS-based IDEAS [23]

The Risø-ILL based McStas [26]

Most of these packages are still maintained (to a higher or lower degree) and the status of their projects can be found on their respective homepages [50,53,55].

In the view of the authors of this paper, the state of the art of currently available neutron Monte Carlo ray-tracing packages are summarized in Table 1, including points of view on differences and main advantages and disadvantages.

Table 1

Overview of state of the art in Monte Carlo ray-tracing, with main advantages and disadvantages of the most widely used packages

Package	Plusses	Minusses
VITESS	∙ Easy for beginners – no code to write.∙ Use of pipes gives automatically parallel execution.∙ Is on track to become hosted at FZJ.	∙ Modules historically written without common structure requirement/standard.∙ Hard to modify at the physics/module level as a user.∙ Open source, but without open software repository.∙ Key staff has at periods had limited backing from hosting institution.
RESTRAX SIMRES	∙ Has mode for “upstream” simulations from sample to source. Gives a clear speed boost for sparse phase-space simulations.∙ Historically specialized for TAS.	∙ Historically specialized for TAS.∙ Open source, but without open software repository.∙ Only few developers
NISP	∙ No assumption that instruments should be linear. Any volume can scatter to any volume.∙ Monte Carlo engine derived from early MCNP.	∙ No assumption that instruments should be linear. (Intrinsically demanding/slow.)∙ Monte Carlo engine derived from early MCNP.∙ Legacy Fortran program with unspecified license terms.∙ Only few developers
McVine	∙ Automatically enherits much functionality from McStas.∙ Instrument description written in pure Python. (Industry standard.)∙ Interfaces with physics models written in Python/c++∙ Object-oriented with separated physics/geometry∙ Open source (BSD) with open repository.	∙ Instrument description written in pure Python. (Not a dedicated solution for neutron scatterers.)∙ Interfaces with physics models written in Python/c++∙ Due to object-oriented nature, components only interact via the neutron∙ Not widely used outside US neutron scattering.∙ Only few developers.
RAMP	∙ New kid on the block.∙ Intrinsically GPU-capable.∙ Written in pure Python.∙ Open Source with open repository.	∙ New kid on the block, i.e. short history so far.∙ Functionality so far limited.∙ Order of magnitude GPU speedup is not yet achieved.
McStas	∙ Flexible, users can interact at instrument, component and c-code level.∙ True Open Source from the beginning, has open repository, many user contributions.∙ New functionality can be added through short, modular c-code.∙ Several developers working in parallel for almost full duration of project history.∙ Parallelised using industry-standard MPI solution. (GPU developments are under way.)∙ Interfaces with many other codes e.g. MCNP, Mantid, iFit.	∙ Extremely flexible, users can interact at instrument, component and c-code level. (Which may be error prone if many advanced features are combined.)∙ True Open Source from the beginning, many user contributions. (There are however requirements for code inclusion.)

Based on the need for instrument developments, the early use of the software packages concentrated on simulating neutron optics, especially guide systems, and concepts for full instruments, including instrument resolutions.

A huge effort was performed in the beginning of this century with systematic comparison between simulation results of guides and optics in different packages and to some extent with actual measured data [58,69]. This work served strongly to track down tricky numerical errors, to increase confidence in the packages themselves, and to establish the ray-tracing technique as such.

With the aim of being able to simulate complete virtual experiments, the last decade has seen a growing demand for – and development of – a suite of scattering sample components.

This is the first article in a series of short reviews of the McStas package – with the planned future articles described in Section 6. For this reason, we will here review the general expressions used for Monte Carlo ray-tracing simulations of neutron scattering, In addition, we will discuss the implementation of the McStas system, its utilization and user community, and the planned future development of the package.

We have aimed this article for neutron scatterers and software developers alike, and we thus begin with a few introductory paragraphs on the nature of neutrons and on the ray-tracing technique. Most material presented here, and other related information, can be found also in the user and component manual for McStas [51].

2. The mathematics of Monte Carlo ray-tracing simulations for slow neutrons

Monte Carlo simulations are in general a way to perform approximate solutions to complex problems by use of random sampling, thereby performing numerical experiments. The first known example of such a method is that of the Buffon needle problem [46] first posed in the 18th century by Georges-Louis Leclerc, Comte de Buffon:

“Suppose we have a floor made of parallel strips of wood, each the same width, and we drop a needle onto the floor. What is the probability that the needle will lie across a line between two strips?”

As explained in a wikipedia article [56], the thinking behind that numerical experiment can in fact be used to construct a Monte Carlo algorithm for the estimation of π. In terms of an actual mathematical algorithm, and later code for computers, the first formulation was in fact to facilitate calculations of neutron physics for the Manhattan project [29].

As with all stochastic methods, Monte Carlo methods are prone to statistical variances. To reduce the statistical errors, a number of variance reduction methods have been designed, of which we here mention two. In Stratified sampling, one divides the parameter space up into mutually exclusive strata, and sample a given amount of each stratum. In Importance sampling, one will choose to sample more often in regions more important for the final result. As an analogy, to calculate the average height over sea level of a landscape, one would map the ragged mountain slopes more carefully than the still water in the mountain lake.

The method relevant for neutron simulations is known as Monte Carlo ray-tracing. This is used to study objects which travel along a path (a ray), and can be (partially) absorbed or scattered into another direction, but not converted into other types of radiation. The most known example of this is visible light, and this type of ray-tracing is frequently used to generate realistic illumination in 3D computer graphics, see e.g. Pov-Ray [33].

2.1. The semiclassical representation of neutron coordinates

In neutron ray-tracing simulations a neutron is represented semiclassically by simultaneously well-defined position, r, velocity, v, time, t, and all three components of the neutron spin vector, s.

Formally, this approach violates the laws of quantum mechanics, in particular the Heisenberg uncertainty relations [28], given for pairs of conjugate position/momentum variables by $$\begin{array}{}\text{(1)}& \mathit{\delta }\mathit{E}\mathit{\delta }\mathit{t}⩾{\displaystyle \frac{\hslash }{2}};\mathit{\delta }\mathit{x}\mathit{\delta }{\mathit{p}}_{\mathit{x}}⩾{\displaystyle \frac{\hslash }{2}}\mathrm{,}\end{array}$$ and similarly for the conjugate variables in the y- and z-directions. However, the semiclassical approximation works really well for describing instruments that use “typical” neutrons with velocities of the order $\mathit{v}\approx 100-10000$ m/s. For example we can imagine a neutron beam with the typical value $\mathit{\lambda }=4$ Å ( $\mathit{v}=1000$ m/s) and a low divergence of ${0.06}^{\circ }$ (0.001 radians) passing a tight slit of 0.1 mm width. Here we have $\mathit{\delta }\mathit{x}={10}^{-4}$ m and $\mathit{\delta }{\mathit{p}}_{\mathit{x}}={\mathit{m}}_{\mathrm{n}}\mathit{\delta }{\mathit{v}}_{\mathit{x}}=0.001{\mathit{m}}_{\mathrm{n}}{\mathit{v}}_{\mathit{x}}\approx 1.7·{10}^{-26}$ kg m/s. Hence, the product of the two uncertainties is around $1.7·{10}^{-30}$ kg m²/s, or 16000 ℏ, so the uncertainty relation holds with a large margin. As another example, the best energy resolution of a backscattering instrument is around $\mathit{\delta }\mathit{E}\approx 0.5$  μeV for neutrons of 6.6 Å, and the best time resolution is the time is takes one of these neutrons to transverse a small (2 mm) sample, $\mathit{\delta }\mathit{t}\approx 3$  μs. The product of this is $1.5·{10}^{-12}$ eVs, or 400 ℏ. Hence, the uncertainty relation also here holds well.

The validity of the semiclassical approach for neutron scattering is also discussed in detail in Ref. [30]. The similar issues with the semiclassical approximation of the neutron spin will be discussed in the future review on polarized scattering.

Fig. 1.

Neutron ray interacting with a reflecting surface: 1. The neutron begins with parameters $\mathit{r}$ , v, p (black) at a distance from the surface. 2. The neutron is propagated to the surface and now has parameters ${\mathit{r}}^{\prime }$ , v, p (blue). 3. The neutron is reflected and achieves updated parameters ${\mathit{r}}^{\prime }$ , ${\mathit{v}}^{\prime }$ , ${\mathit{p}}^{\prime }$ (red).

2.2. The weight factor

The physical unit of a neutron ray is neutron rate, or neutron intensity (s⁻¹). On a typical neutron source, the emission rate can amount to 10¹⁵/s; an order of magnitude larger than the number of rays that can typically be simulated on a laptop in one year. For this reason, to simulate realistic values for neutron count numbers, a neutron ray in general represents more than a single physical neutron (per second). To accommodate this, each ray contains an additional parameter, the weight factor, p, with a unit of neutrons per second. When the ray begins at the source, p is typically 10³/s to 10⁹/s. The intensity of the given neutron ray is updated on interaction with components of the neutron instrument, see Fig. 1 for an illustration.

As one example of importance sampling, the weight factor is manipulated through the simulation. For example, when some physical neutrons are “lost” due to e.g. finite reflectivity or absorption, the simulated ray will in general continue in the simulations, while p is adjusted to reflect the correct average physical behaviour. When (if) the neutron ray reaches the detector, p may in practice be less than 1/s. Turning again to Fig. 1, using importance sampling by describing the reflected beam only will not neglect the physical influence of absorption, since we ensure conservation of intensity, i.e. $\mathit{p}={\mathit{p}}^{\prime }+{\mathit{p}}^{″}$ , where ${\mathit{p}}^{″}$ is the intensity of the non-reflected beam (which we in this example do not transport).

We proceed along this line of reasoning to calculate the neutron intensity in a beam. Consider a neutron component at a surface perpendicular to the beam. The neutron intensity is given by the sum of all simulated rays reaching this surface: $$\begin{array}{}\text{(2)}& {\mathit{I}}_{\mathit{j}}={\displaystyle \sum }_{\mathit{i}=1}^{\mathit{N}}{\mathit{p}}_{\mathit{i}\mathrm{,}\mathit{j}-1}\mathrm{,}\end{array}$$ where i is the ray index, N is the total number of rays, and j is the index of the given component. The index $\mathit{j}-1$ on p indicates that we consider the intensity being emitted by the previous component. If the ray does not reach this point, we have ${\mathit{p}}_{\mathit{i}\mathrm{,}\mathit{j}-1}=0$ . The weight of the neutrons after interacting with component j, is expressed by $$\begin{array}{}\text{(3)}& {\mathit{p}}_{\mathit{j}}={\mathit{w}}_{\mathit{j}}{\mathit{p}}_{\mathit{j}-1}\mathrm{,}\end{array}$$ where the ray index, i, is omitted for simplicity. The weight multiplier of the j’th component, ${\mathit{w}}_{\mathit{j}}$ , is calculated by the important probability rule $$\begin{array}{}\text{(4)}& {\mathit{f}}_{\mathrm{MC}\mathrm{,}\mathrm{b}}{\mathit{w}}_{\mathit{j}}={\mathit{P}}_{\mathrm{b}}\mathrm{,}\end{array}$$ where ${\mathit{P}}_{\mathrm{b}}$ is the physical probability for event b, and ${\mathit{f}}_{\mathrm{MC}\mathrm{,}\mathrm{b}}$ is the Monte Carlo sampling probability.

Often, there is only one non-trivial event to consider. This may, e.g., be the case for a neutron beam being attenuated by absorption. Here, each neutron ray either survives or is lost. Absorbed neutrons are typically not simulated ( ${\mathit{f}}_{\mathrm{MC}}=0$ ), while for the other case (transmission) ${\mathit{f}}_{\mathrm{MC}}=1$ . The probability rule thus dictates that ${\mathit{w}}_{\mathit{j}}={\mathit{P}}_{\mathrm{b}}$ .

When a Monte Carlo branch point is reached (selection between several events), we have ${\mathit{f}}_{\mathrm{MC}\mathrm{,}\mathrm{b}}\mathrm{<}1$ for each branch, b. However, since ${\mathit{f}}_{\mathrm{MC}}$ is a probability function, we must have $$\begin{array}{}\text{(5)}& {\displaystyle \sum }_{\mathit{b}}{\mathit{f}}_{\mathrm{MC}\mathrm{,}\mathrm{b}}=1.\end{array}$$ If a variance reduction scheme is used where some branches are simulated less frequently than the physical probability, we have ${\mathit{w}}_{\mathit{j}}={\mathit{P}}_{\mathrm{b}}\mathrm{/}{\mathit{f}}_{\mathrm{MC}\mathrm{,}\mathrm{b}}$ , which may exceed unity. Hence, the weight factor may also increase during the simulation.

2.3. Neutron moderators and focusing

Consider a component that may send neutrons out in all directions, e.g. the moderator or the sample. Often in McStas, we choose uniform sampling when simulating all neutron directions. The sampling probability per solid angle is then: $$\begin{array}{}\text{(6)}& {\mathit{f}}_{\mathrm{MC}\text{ unfocus}}\mathit{d}\mathrm{\Omega }={\displaystyle \frac{\mathit{d}\mathrm{\Omega }}{4\mathit{\pi }}}.\end{array}$$ However, this simulation scheme may cause huge inefficiencies. In many simulations, only rays emitted in certain “interesting directions” have any chance of being detected. In such cases, one will employ the technique of focusing, where the simulated ray will be emitted only within a certain solid angle, $\mathrm{\Delta }\mathrm{\Omega }$ – but with uniform probability within $\mathrm{\Delta }\mathrm{\Omega }$ . Then the Monte Carlo sampling probability is: $$\begin{array}{}\text{(7)}& {\mathit{f}}_{\mathrm{MC}\text{ focus}}\mathit{d}\mathrm{\Omega }={\displaystyle \frac{\mathit{d}\mathrm{\Omega }}{\mathrm{\Delta }\mathrm{\Omega }}}.\end{array}$$ To avoid systematic errors in the simulation results, $\mathrm{\Delta }\mathrm{\Omega }$ should contain all the directions contributing to the neutron intensity to be simulated. When using (4) to obtain the weight factor transformations, we see that the focusing alone will bring the amount $$\begin{array}{}\text{(8)}& {\mathit{w}}_{\text{focus}}={\mathit{w}}_{\text{unfocus}}{\displaystyle \frac{\mathrm{\Delta }\mathrm{\Omega }}{4\mathit{\pi }}}\mathrm{,}\end{array}$$ which is the general focusing contribution to the weight factor. Comparing to the uniform case, the focusing method gives smaller weight factors per ray, but a larger number of rays traveling towards the detector. This gives, on average, the same final result, but with a smaller statistical error. The focusing technique is a typical example of importance sampling.

Imagine a McStas moderator component, simulating a moderator that emits R neutrons per second, divided into N rays. The total intensity in the simulation is found from $$\begin{array}{}\text{(9)}& \mathit{R}={\mathit{I}}_0={\displaystyle \sum }_{\mathit{i}=1}^{\mathit{N}}{\mathit{p}}_{\mathit{i}\mathrm{,}0}.\end{array}$$ Assuming uniform distribution of ray weights, we find the initial weight to be ${\mathit{p}}_0=\mathit{R}\mathrm{/}\mathit{N}$ . If focusing of the emitted rays is taken into account, we reach the expression for the initial weight $$\begin{array}{}\text{(10)}& {\mathit{p}}_0={\displaystyle \frac{\mathit{R}}{\mathit{N}}}{\displaystyle \frac{\mathrm{\Delta }\mathrm{\Omega }}{4\mathit{\pi }}}.\end{array}$$ Most moderators have wavelength-dependent emission rates, where $\mathit{R}=\int \mathit{\varphi }\mathrm{(}\mathit{\lambda }\mathrm{)}\mathit{d}\mathit{\lambda }$ . Imagine that we sample uniformly a wavelength band of width $\mathrm{\Delta }\mathit{\lambda }$ . Then the sampling frequency would be ${\mathit{f}}_{\mathit{\lambda }}=1\mathrm{/}\mathrm{\Delta }\mathit{\lambda }$ . With the same argument as the focusing, we reach the final expression for the moderator weight factors $$\begin{array}{}\text{(11)}& {\mathit{p}}_0\mathrm{(}\mathit{\lambda }\mathrm{)}={\displaystyle \frac{\mathit{\varphi }\mathrm{(}\mathit{\lambda }\mathrm{)}}{\mathit{N}}}{\displaystyle \frac{\mathrm{\Delta }\mathrm{\Omega }}{4\mathit{\pi }}}\mathrm{\Delta }\mathit{\lambda }\mathrm{,}\end{array}$$ which is the typical equation being used in McStas moderator components.

2.4. Estimates of simulation uncertainty

In a simulation, like for real experiments, it is important to be able to estimate the statistical uncertainty. We here present a simple derivation of the uncertainty in simulations with weight factors.

Imagine M rays, all with weight factor p. Each ray is imagined to have an overall probability, ${\mathit{P}}_{\mathrm{d}}$ , of reaching the detector. The number of observed rays will be binominally distributed with a mean value $\overline{\mathit{N}}=\mathit{M}{\mathit{P}}_{\mathrm{d}}$ and variance ${\mathit{\sigma }}^2\mathrm{(}\mathit{N}\mathrm{)}=\mathit{M}{\mathit{P}}_{\mathrm{d}}\mathrm{(}1-{\mathit{P}}_{\mathrm{d}}\mathrm{)}$ [3]. Very often, we will have ${\mathit{P}}_{\mathrm{d}}\ll 1$ , but still sufficiently large values of M, so that $\overline{\mathit{N}}=\mathit{M}{\mathit{P}}_{\mathrm{d}}\gg 1$ . In this case, the distribution will be approximately Gaussian with a variance $\mathit{M}{\mathit{P}}_{\mathrm{d}}$ – or a standard deviation $$\begin{array}{}\text{(12)}& \mathit{\sigma }\mathrm{(}\mathit{N}\mathrm{)}=\sqrt{\mathit{M}{\mathit{P}}_{\mathrm{d}}}\approx \sqrt{\mathit{N}}\mathrm{,}\end{array}$$ where the observed value, N, is used as the best estimate for the true mean value, $\overline{\mathit{N}}$ . Recalling the constant weight factors of the rays, the total simulated intensity becomes $$\begin{array}{}\text{(13)}& \mathit{I}=\mathit{N}\mathit{p}\mathrm{,}\end{array}$$ with a variance $$\begin{array}{}\text{(14)}& {\mathit{\sigma }}^2\mathrm{(}\mathit{I}\mathrm{)}={\mathit{\sigma }}^2\mathrm{(}\mathit{N}\mathrm{)}{\mathit{p}}^2=\mathit{N}{\mathit{p}}^2.\end{array}$$

We now allow the rays to have different, but discrete, weight factors, ${\mathit{p}}_{\mathit{i}}$ . The simulated number of rays with one particular weight factor is denoted ${\mathit{n}}_{\mathit{i}}$ (standard deviation $\sqrt{{\mathit{n}}_{\mathit{i}}}$ ). The total simulated result is now $$\begin{array}{}\text{(15)}& \mathit{I}={\displaystyle \sum }_{\mathit{i}}{\mathit{n}}_{\mathit{i}}{\mathit{p}}_{\mathit{i}}.\end{array}$$ Using the property of uncorrelated events, ${\mathit{\sigma }}^2\mathrm{(}\mathit{a}+\mathit{b}\mathrm{)}={\mathit{\sigma }}^2\mathrm{(}\mathit{a}\mathrm{)}+{\mathit{\sigma }}^2\mathrm{(}\mathit{b}\mathrm{)}$ , we can rewrite the sum to reach the statistical variance of the simulated result: $$\begin{array}{}\text{(16)}& {\mathit{\sigma }}^2\mathrm{(}\mathit{I}\mathrm{)}={\displaystyle \sum }_{\mathit{i}}{\mathit{n}}_{\mathit{i}}{\mathit{p}}_{\mathit{i}}^2.\end{array}$$

Now, let the simulations occur with rays of arbitrary weight factors, ${\mathit{p}}_{\mathit{j}}$ , different for each ray, i.e. ${\mathit{n}}_{\mathit{i}}\equiv 1$ , in practice allowing a continuous distribution of ${\mathit{p}}_{\mathit{i}}$ . If the distribution of these weight factors is reasonably well behaved, we can generalize the equations above to reach $$\begin{array}{}\text{(17)}& \mathit{I}={\displaystyle \sum }_{\mathit{i}}{\mathit{p}}_{\mathit{i}}\mathrm{,}\end{array}$$ with a variance: $$\begin{array}{}\text{(18)}& {\mathit{\sigma }}^2\mathrm{(}\mathit{I}\mathrm{)}={\displaystyle \sum }_{\mathit{i}}{\mathit{p}}_{\mathit{i}}^2\mathrm{,}\end{array}$$ which is, in the proper limits, consistent with (14) and (16). This is the principle for calculating uncertainties of simulation results in the McStas package. In practice, the sum of both ${\mathit{p}}_{\mathit{i}}$ and ${\mathit{p}}_{\mathit{i}}^2$ are being accumulated during the simulation to avoid having to store the full list of ${\mathit{p}}_{\mathit{i}}$ ’s for rays reaching the detectors.

2.5. Scaling to real-world statistics

In order to compare with real experiments, in particular when using simulated data as input to data analysis programs, simulated intensities must be converted to average integer counts using an imagined counting time, T. However, this time cannot be chosen to be arbitrarily high, since the data analysis program will assume $\sqrt{\mathit{N}}$ statistics, and for high values of T, the counting statistics will be better than that of the simulation, causing unreal errors in the data treatment.

Let us now we quantify the effect of scaling with a counting time. The average count number is easily obtained by $$\begin{array}{}\text{(19)}& ⟨\mathit{N}⟩=\mathit{I}\mathit{T}\mathrm{,}\end{array}$$ with I being the simulated intensity. Now, care must be taken when evaluating the standard deviations. Counting statistics require that the variance is $$\begin{array}{}\text{(20)}& {\mathit{\sigma }}^2\mathrm{(}\mathit{N}\mathrm{)}=⟨\mathit{N}⟩\mathrm{,}\end{array}$$ while the simulated value will have the simulated variance ${\mathit{\sigma }}^2\mathrm{(}\mathit{I}\mathit{T}\mathrm{)}={\mathit{\sigma }}^2\mathrm{(}\mathit{I}\mathrm{)}{\mathit{T}}^2$ . We must require that the percieved counting statistics cannot be better than the simulated variance, or $\mathit{N}⩾{\mathit{\sigma }}^2\mathrm{(}\mathit{I}\mathrm{)}{\mathit{T}}^2$ , leading to $$\begin{array}{}\text{(21)}& {\mathit{T}}^2⩽{\displaystyle \frac{⟨\mathit{N}⟩}{{\mathit{\sigma }}^2\mathrm{(}\mathit{I}\mathrm{)}}}.\end{array}$$ For many virtual experiments, the value of T can be surprisingly low, of the order seconds, since in particular the highest counting rates tend to be underrepresented in the simulations.

If a series of simulated values need to be scaled by the same time, T, the maximum time is given by the smallest value of T that fulfills (22), or $$\begin{array}{}\text{(22)}& {\mathit{T}}_{max}^2={Min}_{\mathit{j}}\left\{{\displaystyle \frac{⟨{\mathit{N}}_{\mathit{j}}⟩}{{\mathit{\sigma }}^2\mathrm{(}{\mathit{I}}_{\mathit{j}}\mathrm{)}}}\right\}.\end{array}$$ Now, when performing a scaling of this type, some data points will likely obtain actual errors smaller than counting statistics. To make sure that every point in the series has an errorbar that complies with $\sqrt{\mathit{N}}$ , one will for all points, except the one with the smallest value of $⟨{\mathit{N}}_{\mathit{j}}⟩\mathrm{/}\mathrm{(}{\mathit{\sigma }}^2\mathrm{(}{\mathit{I}}_{\mathit{j}}\mathrm{)}\mathrm{)}$ , need to add a random error, ${\mathrm{\Sigma }}_{\mathit{j}}^2$ : $$\begin{array}{}\text{(23)}& {\mathrm{\Sigma }}_{\mathit{j}}^2=⟨{\mathit{N}}_{\mathit{j}}⟩-{\mathit{\sigma }}^2\mathrm{(}{\mathit{I}}_{\mathit{j}}\mathrm{)}{\mathit{T}}_{max}^2\mathrm{,}\end{array}$$ and then round up or down to an integer.

In conclusion, the count-error-true simulation numbers are $$\begin{array}{}\text{(24)}& {\mathit{N}}_{\mathit{j}}=Round\mathrm{(}{\mathit{I}}_{\mathit{j}}{\mathit{T}}_{max}+{\mathit{E}}_{\mathit{j}}\mathrm{)}\mathrm{,}\end{array}$$ where ${\mathit{E}}_{\mathit{j}}$ is a normally distributed stochastic variable with mean zero and variance ${\mathrm{\Sigma }}_{\mathit{j}}^2$ .

This transformation to integer count numbers is not used in McStas, but must be applied as a post processing of the data after the simulations. Further work is, however, needed on the sampling strategies for future more efficient transformation of errorbars.

3. The concept and implementation of McStas

As also outlined in the introduction, the philosophy and mindset behind starting the McStas project, as laid out by K.N. Clausen, was to a large extent to minimize duplication of code and efforts by allowing code to be shared [26]. The main ideas that were implemented since the first version of McStas were:

• Licensing should follow Open Source standards. Since version 1.8 McStas is released under the GNU General Public License v. 2.0. Earlier versions were licensed under a McStas- and RISØ specific license that was open, but did not allow redestribution.

• The code should be modular to allowed sharing models of instrument parts between users. This evolved into the so-called components, explained below.

• The syntax of the user-supplied instrument definition should be simple and powerful. For the development of the syntax details, K. Nielsen got input to the syntax, using pen and paper before writing any code, from the RISØ staff, in particular K. Lefmann and H.M. Rønnow wrote the first instrument and handful of components in this way.

• The software concepts should first of all make sense to instrument scientists, allowing them to work in a way that felt almost like runnning an actual instrument

• Documentation should as far as possible be embedded within the code

• Simplicity and readability of the code should be preferred over performance.

McStas is implemented in a three-layered structure. The core system and run-time libraries take care of the compiling and execution of the program, user interface, and visualization of the results. The core and run-time c-code libraries typically come in a set of header and c-code files, such as

• mccode-r.h and mccode-r.c, common to both McStas and McXtrace and containing functions for generating random numbers, routines for intersecting particle ray and geometrical objects etc.

• mcstas-r.h and mcstas-r.c, containing physical constants for neutron physics, propagation routines with/without gravitation etc. (McXtrace has a similar set of files implementing X-ray physics)

• interoff-lib.h and interoff-lib.c for handling of surface-polygon based geometries in GeomView OFF file format [41]

• interpolation-lib.h and interpolation-lib.c, routines for interpolating in gridded and sparse datasets

• ref-lib.h and ref-lib.c, routines for describing the physics of reflective super-mirror surfaces

• read_table-lib.h and read_table-lib.c, routines for reading tabulated data from files

• pol-lib.h and pol-lib.c, routines for handling polarized neutron transport

• adapt_tree-lib.h and adapt_tree-lib.c, routines to define adaptive importance sampling schemes

These levels of code are only modified by a handful of system specialists. The core system also manages the main Monte-Carlo loop that in turn produces a specified number of neutron rays for the simulation.

Supplementing the core system run-time codes, a number of c-files are present to support individual McStas components, such as

• monitor_nd-lib.h and monitor_nd-lib.c, supporting the very general Monitor_nD monitor component from Emmanuel Farhi, ILL

• ESS_butterfly-lib.h, ESS_butterfly-geometry.c and ESS_butterfly-lib.c, supporting the ESS_butterfly source component

• Geometry_functions.c, Union_functions.c and Union_initialization.c, supporting the Union component set by Mads Bertelsen, ESS DMSC

• nxs.h and sginfo.h, supporting the Sample_nxs component from Mirko Boin, HZB

• ess_source-lib.h and ess_source-lib.c, supporting the legacy ESS_moderator source component

• chopper_fermi.h, chopper_fermi.c, general.h, general.c, intersection.h, intersection.c, vitess-lib.h and vitess-lib.h, supporting the Vitess_Chopperfermi component, ported from VitESS

The components are modular pieces of code. The McStas components are written in a unique domain-specific language, which is augmented using ANSI-C. A McStas component will typically be equivalent to a beam-optical component, such as a chopper, a collimator, a particular sample, or a detector (often called “monitor” in McStas). Sources and moderator are in McStas implemented as one combined component. Components are typically parametrised, e.g. with respect to size, in order to be of general use. The components will use the matematical tools for e.g. weight transformation and focusing, presented in Section 2.

The instrument is a description of a particular neutron scattering set-up. McStas operates at one instrument at a time. The instrument file contains a selection of components and their spatial position and orientation. The instrument is written fully in the McStas domain specific language, although ANSI-C may be used to augment the functionality. Also the instruments can be parametrised, for example to allow the scanning of a triple-axis instrument by changing a number of angles. Figure 2 illustrates the roles of system, components, and instrument file for a very simple time-of-flight instrument, consisting of a pulsed source, a guide, a chopper, and a time-of-flight detector.

Fig. 2.

The figure illustrates how the McStas (and McXtrace) codes are layered, including how user, developers and code generation mechanism interact with the various layers of code. The geometrical lay-out of the instrument to simulate is shown at the top panel, while the bottom panes represent the different interacting layers of McStas. The instrument file is shown fully in the large panel to the left; the figure is not to scale. Reproduced from [65].

When McStas is executed on an instrument file, the system compiles the information from the instrument file and the component library and produces an ANSI-C file. This is in turn compiled by a standard C-compiler to an executable code, which is called by the system. After the process has terminated, McStas collects the resulting data files, stores them in specified directories, and visualises the data upon request. Everything can be controlled from a user-friendly GUI interface, or by a scripting language.

McStas has a users manual and a component manual [51], both of which are very comprehensive documents. A more readily accessible documentation may be found in the project home page. Documentation of component functionality is furthermore embedded in the header of the component code. This documentation is e.g. used for the generation of online component help at the McStas home page [50].

4. Use and users of McStas

In its first two decades, McStas has seen a steady increase of users. We have chosen to illustrate this with a literature search. We have analyzed the 350 articles citing one (or more) of the basic McStas release articles, and by far most of these articles in fact represent simulation work done by the use of McStas.

In Fig. 3, we show the number of McStas-using articles as a function of publication year. We see that after a slow start (except for a spike in 2002), the publication rate reached 10/y in 2006 and 20/y in 2011. Since 2014, the number has been approx. 30/y, which is partly based on the large instrumentation work related to ESS.

Fig. 3.

The annual number of articles citing McStas shown as a function of publication year.

Fig. 4.

The number of articles citing McStas, divided into use categories, shown as a function of 5-year intervals.

It is highly interesting to study what McStas has been used for, judging from this literature base. Figure 4 shows the six main use categories of the packages, distributed on 5-year periods. Although the overall use of McStas has increased by a factor 4 over the years, the distribution between the categories has remained remarkably constant. Let us present the themes one by one (all will be described in much more detail in later issues of this article series):

Components. These papers describe the simulations of individual beam-optical components, such as choppers, monochromators [36], lenses [14], or mirrors [59]. Around 5% of the articles fall in this category.

Guides. A particular application of McStas has been to study different guide systems, both generic guide geometries [20,40], and guide constructions for actual instruments [15,22]. Recently, tools have become available for semi-automatic optimization of guide geometries [4] and guide coatings [32], strongly enhancing the capabilities and speed of simulation-guided optimizations for guide designs. Around 20% of the McStas articles are in this category.

Instrument papers. This denotes the important papers describing the performance of an actual instrument. Such papers will typically obtain many citations, as they will be the standard reference for users of the instrument. Often, the use of McStas in these papers is related to prediction of flux on sample, and the calculation of resolution functions, both to be compared to the actual measured values. Approx. 10% of the McStas papers are instrument papers and represent instruments at the major facilities around the world, for example:

• SNS (e.g. CNCS and ARCS) [1,11]

• NIST (MACS) [35],

• J-PARC (e.g. AMATERAS and 4SEASONS [19,31]

• CARR (HIPD) [67]

• ANSTO (TAIPAN) [8]

• ILL (D7) [60]

• ISIS (ENGIN-X) [37]

• FRM-2 (KWS-2) [34]

• HZB (FLEX) [22]

• PSI (RITA-II) [24]

Instrument simulations. These articles represent the design work behind instruments, some of which are not yet built, and a few of which may never be built at all. This is a main category for McStas work, containing 35% of all articles. The simulations cover all facilities listed above, as well as:

the Algier reactor

the Argentine reactor

CSNS

Dhruva

ESS

HANARO

Dubna

LLB

PIK

Most of the instrument simulation publications are related to the newest sources: ESS (24 papers), J-PARC (20 papers), and FRM-2 (13 papers), but also much activity is related to ILL (12 papers) and SNS (7 papers). Examples are too plenty to be given at this stage, but we will elaborate on the topic later in this review series.

As a separate task, McStas has been used in two demanding studies of the complete instrument suite for ESS, first to optimize the pulse structure [25], then to optimize the moderator geometry [2].

Virtual experiments. These articles represent the most advanced use of simulations, namely the performance of a full virtual experiment on a simulated instrument [13,27]. This has been used for the determination of accurate resolution functions [9,63], tools for teaching practical neutron scattering without the investment of beam time [62], for correcting for multiple scattering events in samples [5,12,13], and even for analysis of complete data sets [6,16,61,66]. Around 15% of the articles fall in this category.

Packages. These articles describe simulation packages, typically McStas, but also VITESS and RESTRAX have cited McStas, partially for the reason that these three packages have inspired and cross-fertilized each other over the last two decades. In addition, McStas has spawned an X-ray cousin (McXtrace) [21] and a He-scattering cousin [10]. Around 15% of the articles are in this category.

In addition to this bibliometric study, we have performed a survey of the McStas user community during the autumn of 2018. Ten questions were posed via the SurveyMonkey platform [42] in the period from November 6th 2018 to December 3rd 2018. All questions were voluntary/optional and contained the possibility of free-form input via an ‘other’ text field. In total 52 responses were received, out of which 47 filled demographic data. The full McStas user base is estimated to be approximately 400, based on the current 252 members of the McStas user mailinglist and knowing from experience that not all users register.

The main findings of the survey are that:

The McStas users are to a high degree running the most recent version of McStas, at the time of survey v. 2.4.1.

Linux is the dominant operating system, followed by Windows 10 and macOS. Some users are now running the Linux binaries under Windows 10 via the Windows Subsystem for Linux [44].

McStas users are distributed over all of the continents that have access to neutron facilities.

McStas users are mostly based at neutron facilities, but the package also finds use at universities and in industry.

Most users are now routinely using our modernised Python based gui and command-line utitlites.

The vast majority of McStas users are very happy about the McStas installation process, daily use and the available user support.

For those interested in the full details, a full report of the survey is published at the McStas website [43].

5. The future of McStas

As our user survey demonstrates, McStas is as a project alive and healthy with a 20 year long history of development and a stable user community. What separates McStas from the alternative simulation packages is mainly the applied code-generation technique, the quality of documentation and the strongly collaborative nature of the package, allowing users to contribute. To stay alive and healthy however, a sustained effort is needed in terms of manpower, resources, support and development.

The main focal points for the development in the coming years are:

Lower the typical release cycle of 12–18 months to 3–6 months. Updating more often lowers ‘time to market’ for new ideas, solutions and bug-fixes. McStas 2.5 was released in December 2018 and 2.5.1 will be released in the spring of 2019.

Make our test/validation effort more stringent and use more of a ‘continuous integration’ approach to testing. All components should have output test data that we can monitor against on a continuous basis. Today only a scalar numerical output defines a ‘test’ – in the future we would like to compare at a more detailed level of data in individual bins of individual monitors.

Enrich the sample library even further. A natural layer of growth is the Union [5] concept by Mads Bertelsen which allows to separate geometry and physical processes by defining an inner Monte Carlo loop (e.g. for describing complicated arrangements of matter in/around the sample area). The Union components should become complete in terms of describing all possible geometries and all possible physics, which is today not fully the case.

A modernization of the code generation layer is under way, which will allow us to move to a C-function based grammar rather than today’s #define-based approach.1

¹

The classical McStas code uses c-preprocessor [45] define instructions to locally change the meaning of a variable or parameter name. This means that e.g. a geometrical parameter symbol xwidth can be used by multiple components within the same instrument, but with a local meaning and without inconsistency.

This will break certain backward-compatibility and move us to McStas v. 3.0

The modernized code-generation will ease the process of targeting modern GPU-based architectures.

The reader is further encouraged to consult our standard references [65] and [64] for a relatively recent and complete description of McStas.

6. The McStas review series

This article should be seen as an introduction to a series of McStas review papers. Planned themes in this series cover the main uses of McStas described above:

Components

Guide systems

Instrument simulations and virtual experiments

In addition to a review of these important cases, we also will describe a few more technical matters that deserve a more thorough presentation than what has been given in literature so far:

• Modeling of scattering from samples

• Simulation of polarized neutrons, including spin precession

Is is our aim that this article series on McStas will serve to share knowledge of the package and its utilization, as well as enhance the capabilities for designing and modeling neutron instrumentation worldwide.