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

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 (1) δ E δ t ⩾ ℏ 2 ; δ x δ p x ⩾ ℏ 2 , 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 v ≈ 100 − 10000 m/s. For example we can imagine a neutron beam with the typical value λ = 4 Å ( v = 1000 m/s) and a low divergence of 0.06 ∘ (0.001 radians) passing a tight slit of 0.1 mm width. Here we have δ x = 10 − 4 m and δ p x = m n δ v x = 0.001 m n v x ≈ 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 δ E ≈ 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, δ t ≈ 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.

OPEN IN VIEWER

Fig. 1.

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

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. p = p ′ + p ″ , where 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: (2) I j = ∑ i = 1 N p i , j − 1 , where i is the ray index, N is the total number of rays, and j is the index of the given component. The index 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 p i , j − 1 = 0 . The weight of the neutrons after interacting with component j, is expressed by (3) p j = w j p j − 1 , where the ray index, i, is omitted for simplicity. The weight multiplier of the j’th component, w j , is calculated by the important probability rule (4) f MC , b w j = P b , where P b is the physical probability for event b, and f MC , 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 ( f MC = 0 ), while for the other case (transmission) f MC = 1 . The probability rule thus dictates that w j = P b .

When a Monte Carlo branch point is reached (selection between several events), we have f MC , b < 1 for each branch, b. However, since f MC is a probability function, we must have (5) ∑ b f MC , b = 1 . If a variance reduction scheme is used where some branches are simulated less frequently than the physical probability, we have w j = P b / f MC , 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: (6) f MC unfocus d Ω = d Ω 4 π . 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, Δ Ω – but with uniform probability within Δ Ω . Then the Monte Carlo sampling probability is: (7) f MC focus d Ω = d Ω Δ Ω . To avoid systematic errors in the simulation results, Δ Ω 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 (8) w focus = w unfocus Δ Ω 4 π , 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 (9) R = I 0 = ∑ i = 1 N p i , 0 . Assuming uniform distribution of ray weights, we find the initial weight to be p 0 = R / N . If focusing of the emitted rays is taken into account, we reach the expression for the initial weight (10) p 0 = R N Δ Ω 4 π . Most moderators have wavelength-dependent emission rates, where R = ∫ ϕ ( λ ) d λ . Imagine that we sample uniformly a wavelength band of width Δ λ . Then the sampling frequency would be f λ = 1 / Δ λ . With the same argument as the focusing, we reach the final expression for the moderator weight factors (11) p 0 ( λ ) = ϕ ( λ ) N Δ Ω 4 π Δ λ , 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, P d , of reaching the detector. The number of observed rays will be binominally distributed with a mean value N ¯ = M P d and variance σ 2 ( N ) = M P d ( 1 − P d ) [3]. Very often, we will have P d ≪ 1 , but still sufficiently large values of M, so that N ¯ = M P d ≫ 1 . In this case, the distribution will be approximately Gaussian with a variance M P d – or a standard deviation (12) σ ( N ) = M P d ≈ N , where the observed value, N, is used as the best estimate for the true mean value, N ¯ . Recalling the constant weight factors of the rays, the total simulated intensity becomes (13) I = N p , with a variance (14) σ 2 ( I ) = σ 2 ( N ) p 2 = N p 2 .

We now allow the rays to have different, but discrete, weight factors, p i . The simulated number of rays with one particular weight factor is denoted n i (standard deviation n i ). The total simulated result is now (15) I = ∑ i n i p i . Using the property of uncorrelated events, σ 2 ( a + b ) = σ 2 ( a ) + σ 2 ( b ) , we can rewrite the sum to reach the statistical variance of the simulated result: (16) σ 2 ( I ) = ∑ i n i p i 2 .

Now, let the simulations occur with rays of arbitrary weight factors, p j , different for each ray, i.e. n i ≡ 1 , in practice allowing a continuous distribution of p i . If the distribution of these weight factors is reasonably well behaved, we can generalize the equations above to reach (17) I = ∑ i p i , with a variance: (18) σ 2 ( I ) = ∑ i p i 2 , 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 p i and p i 2 are being accumulated during the simulation to avoid having to store the full list of p 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 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 (19) ⟨ N ⟩ = I T , with I being the simulated intensity. Now, care must be taken when evaluating the standard deviations. Counting statistics require that the variance is (20) σ 2 ( N ) = ⟨ N ⟩ , while the simulated value will have the simulated variance σ 2 ( I T ) = σ 2 ( I ) T 2 . We must require that the percieved counting statistics cannot be better than the simulated variance, or N ⩾ σ 2 ( I ) T 2 , leading to (21) T 2 ⩽ ⟨ N ⟩ σ 2 ( I ) . 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 (22) T max 2 = Min j { ⟨ N j ⟩ σ 2 ( I j ) } . 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 N , one will for all points, except the one with the smallest value of ⟨ N j ⟩ / ( σ 2 ( I j ) ) , need to add a random error, Σ j 2 : (23) Σ j 2 = ⟨ N j ⟩ − σ 2 ( I j ) T max 2 , and then round up or down to an integer.

In conclusion, the count-error-true simulation numbers are (24) N j = Round ( I j T max + E j ) , where E j is a normally distributed stochastic variable with mean zero and variance Σ 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.