Prompt: Probability-Conserved Cross Section Biasing Monte Carlo Particle Transport System

An open source software package for simulating thermal neutron propagation in geometry is presented. In this system, neutron propagation can be treated by either the particle transport method or the ray-tracing method. Supported by an accurate backend scattering physics engine, this system is capable of reproducing neutron scattering experiments in complex geometries and is expected to be used in the areas of instrument characterisation, optimisation and data analysis. In this paper, the relevant theories are briefly introduced. The simulation flow and the user input syntax to control it are provided in detail. Five benchmarking simulations, focusing on different aspects of simulation and scattering techniques, are given to demonstrate the applications of this simulation system. They include an idealised total scattering instrument, a monochromatic powder diffractometer, a neutron guide, a chopper and an imaging setup for complex geometries. Simulated results are benchmarked against experimental data or well-established software packages when appropriate. Good agreements are observed.


Introduction
Along with the instrument design activities at the powerful neutron spallation sources over the last decades [Tay06, MAA + 06, KYN + 19, LBC + 11, WCC + 09], Monte Carlo ray-tracing codes emerged as the de facto standard for neutron instrument characterisation and optimisation. There are many packages currently available, such as McStas [LN99], VITESS [ZLM02], and RESTRAX/SIM-RES [ŠK97,SK04]. Most of them employ the linear chain method, in which an instrument is decomposed into a series of components and the neutron beam is only allowed to transmit from an upstream component to the adjoining downstream component. This method is computationally efficient and proven to be highly valuable in estimating neutron instrument flux on the sample [LNA + 06, CKL + 18, SDJE06, OMR + 04] and resolution functions [ASL + 12, KNI + 11, EPN + 11]. A recent review on the state-of-the-art of these packages can be found in [WL20].
It is worth noting that the linear chain method is not enforced in NISP [DSTH99,SD04], a package developed in the 1990s. This package inherits a part of the Monte Carlo engine from an early version of a particle transport code [WL20].
Due to the absence of the linear assumption, neutrons in this package are not travelled in a pre-defined chain, hence the simulation is more realistic, especially when the components can not be arranged in an upstream-downstream manner.
For instance, in an assemble of a sample environment vessel containing a sample, the unwanted scatterings from the vessel can occur before or after neutrons passing by the sample.
Despite fewer assumptions being made, NISP is not currently among the most popular packages. Apart from the computational efficiency, license and manpower shortage problems outlined in [WL20], a limited number of scattering models in bulk materials could also be a major disadvantage. The departure between the model and actual physics introduces errors in the cross section and leads to inaccurate results. Since 2015, NCrystal [CK20], a general purpose open-source engine for neutron scattering has been under active development.
In this package, microscopic cross sections are calculated based on material basic properties [CKKM19, KC21,KC23]. Thanks to the implemented robust numerical methods, neutron scattering can be modelled in a wide energy range, from ultracold to tens of electronvolts.
The presented software package, probability-conserved cross section biasing Monte Carlo particle transport system (Prompt), is our first attempt to tackle the cross section challenge. Available neutron scattering models and the probability-conserved cross section biasing technique are briefly introduced in Section 2. Section 3 explains the concepts of the system along with code examples and shows a minimal but complete simulation input script. Section 4 demonstrates overall system performance by studying a few problems, including an idealised total scattering instrument, a monochromatic powder diffractometer, a neutron guide, an imaging setup for complex geometries and a disk chopper. This work is discussed and concluded in section 5.

Theories
The objective of the system is outlined in section 2.1. Section 2.2 and 2.3 introduce the neutron scattering length and scattering types, respectively. As the theory of neutron scattering is given elsewhere in more detail, e.g. [Lov84,Squ12,Sch14], these sections are dedicated to introducing the physical meanings of the relevant parameters and the underneath cross section models. Section 2.4 introduces the implemented probability-conserved cross section biasing method that is able to effectively control the interaction number of neutrons in a finite volume. Section 2.5 introduces the mirror and chopper models.  where r is the position, p 0 is the momentum, t is the time, s is the spin, and G function is the propagator.

The scope of the system
This higher-dimensional integral is typically evaluated using Monte Carlo particle transport systems, in which the propagator is decomposed into geometry and physics modules in general. In such systems, the physics module predicts the particle states by sampling the physical distributions for given prior states, materials and any applied external fields; while the geometry model traces individual particles in volumes along their trajectories.
On the other hand, the ray-tracing method simplifies the integration by a scaling function , which represents the weight adjustment of the particles. In addition, when the linear chain approximation is applied, the function is set to zero if the leaving particles do not intersect with the detector.
f det (r, p, t, s) = (r, p, t, s; r 0 , p 0 , t 0 , s 0 )f src (r 0 , p 0 , t 0 , s 0 ) This method describes neutron transmission behaviour reasonably well in many optical components. One of the simplest examples is the rotating neutron disk chopper. The function is either unity or zero, representing respectively the event of neutron transmission or absorption.
The objective of Prompt is to solve Eq. 1 in bulk materials and Eq. 2 to model moving optical components, i.e. choppers. The linear chain approximation is relaxed in the system. Only thermal neutrons are supported in this version of code.

Neutron scattering length
The introduction starts with a limiting case, where neutrons scatter with a stationary nucleus that is fixed at a position throughout the scattering process.
If the target nucleus remains in the ground electronic state, and the neutron wavelengths are much longer than the range of neutron-nucleus interaction, the scatterings are contributed only by the s-wave components in terms of partial waves [Sch14]. The scatterings in that case are essentially elastic isotropic and can be characterized by a single parameter b, the so-called bound scattering length [Lov84]. The cross section for this hypothetical 1 case is simply 4π|b| 2 .
In the opposite limiting case, where the motion of a nucleus is not hindered at all, the process can be viewed as elastic scattering in the centre of mass frame of reference. The so-called free scattering length for this case can be derived [Squ12] to be b f = bM/(m + M ), where m and M are the mass of the neutron and target nucleus, respectively. When the incident neutrons are far more energetic than the chemical bonds in the system, say above a few electronvolts, the scattering cross section is approaching the limiting value 4π|b f | 2 , for instance, see ref. [DGM14].
The angular distributions in the laboratory and centre-of-mass frames can be correlated as Eq. 3 [Sch14]. From the viewpoint of this equation, the angular distribution of the bound situation is a special case for the free situation when γ is approaching zero, i.e. the mass of the target nucleus is infinity in comparison with a neutron.
1 Fixing a stationary nucleus is only sound classically, but violates the uncertainty principle.
However, neither the bound nor the free cases can describe the neutron scattering process adequately accurately in a finite and correlated system, especially when incident neutron wavelengths are comparable with the shortest inter-atomic distances.
Notice that the scattering length b is complex in general.

Neutron scattering cross section
Neutral particles exhibit step-like behaviours in materials. In the Lambert-Beer law, the probability of observing a free particle at length x in a homogeneous material can be described as where ρ is the scatter number density and σ is the cross section. The cross section can be calculated from the double differential cross section.
Under the Born approximation and Fermi's golden rule, the thermal neutron scattering cross section is the double integral of the scattering function that scaled by the ratio of the scattered and incident wavenumbers, , to obey the conservation law of momentum and energy. Here Q and ω are the reduced momentum and energy transfers, respectively.
The evaluation and implementation of Eq. 6 depend on the orientation distributions of the statistically equivalent subsystems with respect to the laboratory frame of the simulation system [CK20,KC21]. The orientation should be aligned with the laboratory frame of the simulation, if it is not isotropic, see more discussions in section 3.5.
The neutron scattering function can be decomposed into a series of S j,j ( Q, ω), the partial scattering function, contributed by the j th and j th nuclei pairs [VH54].
Here N is the number of nuclei, b j b j emphasizes the averaging is performed over the isotopes and spins of all the subsystems.
For a more realistic treatment of the scattering, in the case that the value b for the j th nucleus has no correlation with the isotope and spin, the averaging Therefore, In Eq. 10, the first term is originated from the pairs that are build by different nuclei and is known as the distinct scattering. The second term is contributed from a single nucleus solely, hence called the self scattering.
where S coh and S inc are the coherent and incoherent scatterings, respectively.
Coherent scattering can be physically interpreted as the scattering that would arise if all the scattering lengths equal to their mean value. While the incoherent scattering, is proportional to the self scattering, represents the deviations of the scattering lengths from their mean value. It is worth emphasising that the coherent scattering is contributed by distinct scattering as well as a part of the self scattering.
Coherent and incoherent scatterings can be further categorised by whether the neutron exchanges energy with the system, i.e. elastic and inelastic with respect to the laboratory frame of reference.
W is the Debye-Wallet factor, v 0 is the unit cell volume, τ is the lattice point in the reciprocal space, and R is the equilibrium position.
For idealised materials, the orientations of crystal unit cells are described by distribution functions. Numerical procedures are implemented to compute the effective scattering functions averaged over the given distribution. The averaging is rather sophisticated and outside the scope of this paper, refer to Ref. [KC21] for the detailed methods. The backend scattering engine supports the elastic scatterings in isotropic powders, mosaic single crystals and mosaic layer crystals.

Inelastic scatterings
Limited by the memory footprint and initialisation time of typical Monte Carlo simulations, inelastic scattering is only supported in isotropic materials, e.g. powder and liquid.
It has been shown that the incoherent part of the inelastic scattering for ideal crystal powders in the harmonic approximation [Sjö58] and liquids in the Gaussian approximation under the fluctuation-dissipation theorem [RSS62] can arrive independently to the identical form, Here the scattering function depends on the modulus of Q due to the isotropic nature of the material, and Γ(t) is width function.
where ρ(ω) is the density of states, k B is the Boltzmann constant and T is the temperature. The expressions indicate that it is adequate to calculate the incoherent cross section in isotropic materials using only the density of states.
The numerical evaluation of Eq. 17 for crystals is computationally efficient when applying the phonon expansion method [Sjö58], therefore it is suitable to evaluate on-the-fly. On the other hand, for liquids, Γ becomes divergent when ω approaches zero, because of the diffusive motions. The scattering function, in that case, is typically pre-computed based on model-assisted evaluation methods [DGM14] or more expansive direct convolutions [DC22].
The coherent inelastic scattering function may also be included in the precomputed input cross section. In the case of liquid, the Sköld approximation [Skö67], i.e. Eq. 19, is often applied .
, ω S jj (Q) is the partial static structure factor and equals to S jj (Q, ω)dω. As suggested by the zeroth order sum rule of the scattering functions, it is unity when j = j , i.e. the case of self scattering function.
For numerical simplification, the isotropic incoherent approximation, or otherwise known as the Gaussian approximation [RSS62], is often employed to skip the computation of the coherent inelastic scattering function. In that case, the distinct scattering functions are vanished, hence the overall inelastic scattering is approximated entirely by the inelastic part of the self scattering function (see Eq. 10). The absence of distinct scattering results in a smoothened inelastic scattering function with respect to Q. Consequently, the calculated angular distributions are not expected to be directly comparable with those from measurements in every detail. Nonetheless, in practice, the incoherent approximation is known to be highly useful when the angular distribution is less of concern, in the sense that the final results of interest are the averaged values in a large angular range, such as when measuring phonon density of states (see section 9.16 of [Sch14]), characterising neutron moderation processes [DGM14] and estimating the intensities of multi-phonon scattering background (see the section 3.10 of [Squ12]).
In additional to the incoherent approximation, a Debye model can be used for crystals, when only the structural information is mainly concerned, e.g. powder diffraction of a strong coherent sample. The density of states in that case is proportional to the power of phonon energy and the upper cut-off is defined by the Debye temperature [SSP95].
In the situations where the correlations of the systems are relatively insignificant, e.g. in strong absorbing or gaseous materials, the short-collision-time approximation [MK10,CKKM19], can be used to describe the scattering processes, see section 3.5 for more discussions.

The cross section biasing mechanism
According to the probability expressed in Eq. 4, the step length of a free particle can be sampled according to where ξ is a number generated in a uniform distribution [0, 1).
If an arbitrary positive scaling factor g is applied to the cross section to form an artificial cross section σ g = gσ(E), the sampled step length becomes In comparison with l, l g is increased or decreased, when the factor is smaller or greater than unity, respectively. That may effectively control the particle interaction number in a finite volume.
However, altering such an important parameter distorts the step length distribution in Eq. 4 and produces wrong statistical results.
It can be seen that by differentiating Eq. 4, and then dividing the result by When applying an artificial cross section, this equation becomes By the quotient rule of derivative, the impact of applying a biasing factor to the probability along the path can be quantified as It has been shown [MW12] that the distribution distortion can be eliminated when the particle weight is corrected by a compensation factor at the end of each interaction as Here, W o is the original weight of the particle; n is the number of interactions, the step length of the m th step is x m , N p is the number of reaction processes, the biasing factor for the i th process at the m th step is g m,i , the l th process occurs at the end of the m th step, and the corresponding biasing factor is b m,l .
The initial aim of the implementation is to facilitate the study of scatterings between an incident neutron beam and a small sample of a scattering instrument. In such a scenario, multiple scatterings are often considered to be a type of weak but important noise. Applying a biasing factor greater than unity, not only the statistics of the single scattering signal can be improved, but also the variances of the multiple scattering noises can be significantly reduced. Detailed validation and performance of this method are discussed in section 4.1 and 4.4.

Neutron mirror and disk chopper
The empirical reflectivity model in McStas [WFK + 22] is used for neutron mirrors. The reflection probability is a piecewise function, where R 0 is the low-angle reflectivity, Q c is the critical scattering wavenumber, α is the reflectivity slope, m is the m-value of the mirror, and W is the cut-off width. The default values retaken from McStas are given in Table 1. When a neutron is reflected, its energy remains unchanged, so Q only depends on the reflection angle. The reflected direction, n f , can be calculated using the neutron direction n i , and the normal of the contact point n r as The neutron disk chopper is a neutron optics component that modelled in ray-tracing method. The neutron weight is modified by the ray-tracing process instead of detailed scattering and propagation in the transport mode.
When a neutron arrives at the disk chopper, the ray-tracing method is used.
The corresponding function in Eq. 2 is a mask of position and time. It is unity if a neutron hits a opening window, and is zero otherwise. In the implementation, the chopper is considered to be infinitely thin.

The program
The executable, prompt, reads user input files in the Geometry Description Markup Language (GDML) format [CMPS06] and launches the simulation. The language schema is introduced in the manual of the standard in detail [gdm22].
And the available options are summarized in Table 2.
In a typical simulation, statistical distributions of particle histories are iteratively accumulated in a Monte Carlo loop of events. Each event begins with the generation of the particle, and ends with the absorption of the particle in

Option
Default Description 4096 Set the seed for the random generator.

-n [number ]
100 Set the number of primary events.

-v
The flag to activate the visualisation of geometries and trajectories. materials or escape of the particle from all the volume boundaries. The accumulators are referred to as scorers in this system and are realised by histograms that are attached to volumes. The accumulation processes are triggered by the concept of particle tracing states (PTS) on the volume.
The simulation flow is illustrated in section 3.1. Available particle generators, i.e. particle guns, are presented in section 3.2. The volume and scorers are introduced in section 3.3 and 3.4, respectively. The syntax for initialising physics models is described in section 3.5. Finally, a short example of input script is listed and explained in section 3.6.

Simulation flow
The flow of the simulation is divided into two phases, the tracing phase and the interaction phase, as illustrated in Fig. 2.
The system first looks for the volume that contains the last particle in the particle stack. The successful determination of a valid volume indicates that the particle is inside the world, i.e. the largest volume that contains all the child volumes in the simulation. Otherwise the simulation for this particle is terminated. The program then asserts the particle alive status before picking the simulation mode for the interaction phase. If the particle is not just generated at a random position by a particle gun, it must arrive at the boundary of a volume at this stage. In that case, the system checks whether a ray-tracing model is defined for the volume, before branching out into the ray-tracing mode or the transport mode.
In the transport mode, the particle is first treated by the mirror, if it is specified. Next, the macroscopic total cross section is evaluated, and the step length is sampled using Eq. 21. If the sampled step is longer than the maximum distance that allowed by the volume, the particle is moved to the geometry intersection point and enters the relocation subprocess. Otherwise, the particle is moved forward to the location where an interaction takes place. In this case, the particle enters the propagation subprocess or the absorption subprocess, depending on the sampling on the cross sections. In any cases, whenever the particle is moved, the weight is adjusted according to Eq. 25, if the cross section is biased.
In the ray-tracing mode, the particle is simulated in Eq. 2. The main difference between this mode and the transport mode is that the particle states in the volume is completely manipulated by the implemented model. The system views the ray-tracing volume as a black box and has no knowledge about the particle states, until the model decides to eliminate the particle or sends the particle to the boundary of the volume. Notice that, a volume containing child volumes is not supported by the ray-tracing mode, as the internal structure of the volume can not be modelled systematically.
The sub-routines of the interaction phase are presented in Fig. 3. They can modify particle parameters such as location, energy and direction of motion, as well as changing the particle tracing states that can trigger available scorers.
More details of the scorers are introduced in section 3.4.
In the relocation sub-routine shown in Fig. 3a, the particle is moved to the boundary of the next volume, with the momentum remains unchanged. In the propagation sub-routine shown in Fig. 3b, the scattering angle and energy of the particle are sampled. In the absorption sub-routine shown in Fig. 3c, the alive status of the particle is set to false.
Particle(s) being generated with particle gun and stored in stack  Two basic particle guns, which can create neutron point sources, are provided, namely, SimpleThermalGun and IsotropicGun. There are three ways to specify neutron energies in these two guns, using a single value, sampling from uniform distribution or a Maxwell-Boltzmann distribution of coefficient kT .

Geometry
Prompt uses VecGeom [A + 05] v1.2.0 as the geometry engine and GDML [CMPS06] as the language to define the geometry of a simulation. A geometry is a hierarchy of volumes, each of which a shape and a composition should be assigned to. There are three essential parts to describe a volume, namely solids for the shape description, structure for the hierarchy relationship and material for the composition. In this section, the methods to define the solids and structure are given, and leaving the materials part in section 3.5.
The solids currently supported in Prompt are retaken from the GDML user's guide [gdm22] and presented in Appendix A. A solid defines solid type and the corresponding parameters. As an example, the lines presented in Listing 6 define a box of the identical length each side. The length unit, lunit, is millimeter by default. It can be otherwise chosen to be nanometer (nm), micrometer (um), centimeter (cm), meter (m) or kilometer (km). The default unit for angle, aunit, is degree (deg), and radian (rad) can also be used, alternatively. <solids> <box lunit="mm" name="Cube" x="10.0" y="10.0" z="10.0" /> </solids> Listing 6: A box named Cube of 10 mm for each side length The structure is a hierarchy of volumes. A volume is called logical when it refers to a solid and a material. A logical volume becomes physical when it is positioned into a mother logical volume. Notice that any boundary of a daughter volume should contain within its mother volume and boundary overlap test is not implemented in this version of the system. All volumes in a simulation are physical because they have to be positioned, except for the world, the largest volume on the top of the hierarchy. An example is presented in Listing 7, where a hierarchy of two levels is defined. The first level corresponds to the world of this simulation, while the second is a volume named phy 2ndLevel placed in the center of the world.

Scorer
Statistical Information

Scorer
Eight scorers that are designed to accumulate particle statistical information in histograms are listed in Table 3. The variables of interest are generally fixed to one of the attributes of the particles, such as energy, position, and time. For a DeltaMomentum scorer, the method used to calculate the variable should be chosen from elastic or inelastic scatterings. In the elastic case, the momentum transfer Q is calculated as While in the inelastic case, the momentum transfer is To set up, a scorer is attached to a volume from which the statistical information is expected to be extracted. For instance, in Listing 8, a PSD scorer is attached to the volume named Detector. The parameters of scorers are listed in Appendix C.
<volume name="Detector"> <materialref ref="Vacuum"/> <solidref ref="DetectorSolid"/> <auxiliary auxtype="Scorer" auxvalue= "Scorer=PSD; name=NeutronHistMap; xmin=-250; xmax=250; numBins_x=100; ymin=-250; ymax=250; numBins_y=100; ptstate=SURFACE; type=XY"/> </volume> Listing 8: A PSD scorer attached to the volume named Detector The ptstate in Listing 8 denotes the particle tracing state. The particle tracing states reflect the relationships between the particles and volumes, and also act as triggers allowing scorers to accumulate particle weights. Each scorer should be assigned to a particle tracing state to listen to. In Prompt five particle tracing states are defined, as follows: • SURFACE: A particle reaches the geometry boundary of a volume, • ENTRY: A particle entries a volume, • ABSORB: A particle is absorbed in a volume, • PROPAGATE: A particle propagates in a volume, • EXIT: A particle exits a volume.
All of them are available for all scorers, except that MultiScat is immutably set to PROPAGATE.
To illustrate the triggering mechanism of scorers, the transport process of a particle in the geometry of Listing 7 is shown in Fig. 4. The particle transports four steps in this simple geometry. In step 1, when the particle reaches point A at the boundary of volume Sphere, it simultaneously exits volume World and Figure 4: Transport process of a particle in a simple geometry entries volume Sphere. Therefore, it can trigger the scorer that is listening to the EXIT PTS of the volume World and the SURFACE and ENTRY PTS of the volume Sphere. In step 2, when the particle arrives at point B, it is possible to trigger the PROPAGATE or ABSORB PTS of the volume Sphere. In this example, it is chosen to be PROPAGATE so that the simulation continues. Similar to step 1, when the particle reaches point C at the end of step 3, it can trigger the scorer on volume Sphere set to EXIT and the scorer on volume World set to ENTRY.
When the particle reaches point D at the end of step 4, it is at the boundary of volume World and will be out of volume World. Hence, it is killed and can trigger the scorer on volume World set to EXIT.
It is worth noting that the SURFACE and ENTRY can be used interchangeably for scorers. These two PTS are used in the core of the system to indicate whether the interaction is treated using the surface or bulk material models.
See more discussions on surface type physics at the end of section 3.4.

Total scattering
A heavy water sample that is surrounded in a perfect spherical detector is simulated. Fig. 7 shows the visualisation of the geometry and 100 neutron trajectories. In the simulation, a SimpleThermalGun is used to emit neutrons in the positive z-axis direction at 3.635 eV. At 10 m downstream of the gun, a 12 mm radius spherical sample is located in the center of the detector, of which the inner radius is 2000 mm and the thickness is 1 mm. A MultiScat scorer is Hence the simulated structure factor is related to S j,j (Q), the partial static structure factor, used in the cross section generation.
In Fig. 8, the single scattering momentum transfer distribution of the nonbiased sample is compared with the interference differential scattering cross section (DCS) for heavy water measured by Soper [Sop13]. The structures of the simulated result are broadly similar to those in measurement. However, the widths of the peaks are slightly different. Also, the raising of the curve towards the lowest Q is not reproduced by the simulation. The observed discrepancies can be explained by the discrepancies between the static structure factor used in the cross section calculated and that from measurement.
As indicated by Eq. 19, the partial static structure factors should be used for polyatomic systems. To obtain that, the measured data need to be fitted to a microscopic model to determine the atomic structure candidate. Therefore, the peak width discrepancies are likely introduced by the fitting procedure.
Secondly, the raising of the experimental data at low Q is caused by the liquid density fluctuation, which may not be included in the fitting model, hence not included in the cross section.   The results of two runs for the same scattering number are overlapped, hence the probability is conserved when the cross section biasing technique is applied. Although the number of incident neutrons of the biased run is only 10% of that of the non-biased run, their results are almost statistically equivalent. Therefore, the cross section biasing technique is able to improve the computational efficiency significantly.   in between 20°and 160°take-off angles. To eliminate any anisotropic effects, the sample is also made spherical with a 2.5 mm radius. Fig. 11 shows the simulated powder diffraction pattern for the Al 2 O 3 sample. The simulated and measured data are normalised by aligning the 143 peak at 124.1°. The intensity and resolution behaviours of the pattern are well captured by the simulation, and the observed agreements are very good. However, the experimental data is rising towards zero degrees. Similar behaviour is not observed in the simulation. In the experiments, the scattering environment between the sample and the detector is filled with helium; while in the simplified simulation, the environment is effectively vacuum. Likely, this discrepancy is introduced by the fact that the scattering in the gas is not modelled.  Fig. 13. The square shape of the guide is visible and the four corners are highly distinguishable.

Neutron guide
The pattern from the McStas simulation of the same detector appears to be identical to that from Prompt, hence is not shown.
The integral intensity of the rear detector, normalised by the front detector count, is compared with that from McStas. Not surprisingly, as the reflectivity models are identical, the curves from these two packages are overlapped.

Stanford bunny
A Stanford bunny [TL94] is simulated using Prompt in an imaging setup.
The geometry visualisation is shown in Fig. 15. The bunny is in a surface mesh format originally, and the TetGen [Han15] code is used to generate a tetrahedral mesh of 2441 solids in total. In the simulation, the bunny is made of vanadium and approximately 15 cm in height. The bunny is irradiated in a parallel monochromatic 2Å neutrons beam. The beam size is 0.2 m×0.2 m.
Behind the bunny, an ENTRY type position detector is set up to measure the transmission pattern.
Two simulations are performed with two different absorption biases, 1 and 0.1. The images recorded by the detector are shown in Fig. 16. In both Fig. 16a and Fig. 16c, the shape of the object can be clearly observed. The biased simulation produces visibly identical results as the unbiased case. To compare the statistical uncertainties in these two simulations, the areas included in the rectangle of Fig. 16a and Fig. 16c are then zoomed-in to Fig. 16b and Fig. 16d, respectively. Evidently, there are much less fluctuation for the simulation where the absorption bias is set to 0.1 (Fig. 16c and Fig. 16d) than that with bias 1 (Fig. 16a and Fig. 16b), indicating that the cross section biasing mechanism shows its capability in variance reduction in this example.

Outlook
A thermal neutron simulation system is presented. As a standalone system, it is expected to be applicable in the areas of neutron instrument characterisation, optimisation as well as data analysis.
This system is the Monte Carlo module of the China Spallation Neutron Source Simulation System (Cinema). Two associated cross section calculation modules, Trajectory Analysis Toolkit(Tak) [DC22] and PiXiu, are currently under active development. They are essentially the post-analysis modules for molecular dynamics and density functional theory calculations, respectively. It is expected that the overall system should provide a full tool-chain to predict data that are directly comparable with instrument detector raw data from sample atomic structures and potentials. In parallel with Cinema, a machine learning framework is currently in the conceptual design phase. As the simulation tool-chain can describe the instruments and samples as parameters, the initial aim of this framework is for the optimisation of these parameters. To cope with this framework, Prompt is expected to become purely pythonic at the very top level shortly. It could be implemented as a python GDML manipulator or a direct interface to the VecGeom code. In any cases, the system will be backward compatibility with the GDML user input introduced in this paper.