McStas (web), (Manual) is a general tool for simulating neutron scattering instruments. It is based on a compiler that reads a high-level specification language defining the instrument to be simulated and produces C code that performs the Monte Carlo Simulation. This makes it very fast, typical figures are 500000 neutron histories per second on a fast PC.
It is common to test calculations in more than one software to check their validity; an alternative to McStas is Vitess.
McStas shares codebase with its X-Ray analogous, McXtrace, with a similar syntax.
It is recommended to install McStas with miniforge. Once miniforge is installed, create a conda environment and install McStas:
conda create --name mcstas
conda activate mcstas
mamba install mcstas
To use McStas in HPC clusters like Hyperion, it should be installed in the /scratch folder.
Replacing gila by your own username,
module load Python
conda create --prefix /scratch/gila/.conda/envs/mcstas
conda activate /scratch/gila/.conda/envs/mcstas
conda install mamba --channel conda-forge
mamba install mcstas
Once installed, it can be accessed anytime as
module load Python
conda activate /scratch/gila/.conda/envs/mcstas
McStas comes with a comprehensive library of well-tested components that include most standard elements of neutron scattering instruments. New components are constantly created by the community. Neutron sources can also be imported from other programs.
McStasScript is the McStas Python API, really useful to automatize calculations and to make things easier. The ESS workflow using McStas is presented in the ESS DMSC Summer School.
McStasScript provides really useful help commands that can guide you when designing instruments in Jupyter Lab:
import mcstasscript as ms
instrument = ms.McStas_instr("MyInstrument")
instrument.available_components()
# Shows all available component categories
instrument.available_components("sources")
# All available components from a specific category
instrument.component_help("Source_custom")
# Help for a specific component
You can also see the parameters of a component that you already placed,
source = instrument.add_component("MySource", "Source_custom")
monitor = instrument.add_component("monitor", "L_monitor")
source.show_parameters()
# Shows current parameters used in the component
source.set_parameters(xwidth=0.1, ...)
You should configure the instrument settings before running the simulation:
instrument.settings(ncount=1e7, mpi=4, ...)
instrument.show_settings()
# Shows current instrument settings
Finally, run the simulation and plot the data from the monitors that you placed on the instrument:
data = instrument.backengine()
ms.make_plot(data)
The handy functions ms.make_plot() and ms.make_sub_plot() have some customization options, but if you want more control you can just plot the data with Matplotlib:
data_L = ms.name_search("L_monitor", data)
wavelength = data_L.xaxis
intensity = data_L.Intensity
plt.plot(wavelength, intensity)
plt.show()
Since we mostly measure histograms of neutron counts per wavelength or energy range, the binning will be very important to properly normalise the data. This implies not only the number of bins, but also the range of wavelengths or energies that we are measuring. A proper normalisation of the intensities is performed by dividing the intensity by the size of the individual bin. This is the same as multiplying by the number of bins and dividing by the whole measuring range. Additionally, we should also normalise by the detector area, in cm units.
measurement = ms.name_search("L_monitor", data)
I = measurement.Intensity * binning / (Lmax-Lmin) / (xw*100 * yh*100)
If we normalise the intensity from a neutron energy measurement (E_monitor) by the size of the energy bin, the intensity will be in units of energy. On the other side, If we normalise the intensity of a wavelength measurement (L_monitor) by the size of the wavelength bin, and then we convert the x axis to energy, it will be in units of lethargy. Remember the De Broglie conversion from neutron wavelength to energy, \(E = \frac{h}{2m\lambda^2}\) \(E[meV] = \frac{81.82}{(\lambda[AA])^2}\) \(\lambda[AA] = \sqrt{\frac{81.82}{E[meV]}}\)