This page is a jupyter notebook. To download it, please visit: tutorials/first-simulation.ipynb

First simulation

A (Monte Carlo) Simulation consists in (1) options, defining all the necessary parameters to setup the simulation, and (2) results, containing all the outputs of a simulation. One or more simulations form a Project. A pyMonteCarlo project stored on disk has the extension .mcsim. It consists of a HDF5 file and can be opened in the HDFViewer or using any HDF5 library.

Setting up simulation options

The options are defined by the class Options. It contains all the parameters necessary to run one simulation. The parameters are grouped into four categories:

  • program
  • beam
  • sample
  • analyses

The beam, sample and analyses are independent of Monte Carlo programs. In other words, the same sample definition can be used for different Monte Carlo programs. For a given Options instance, only the program needs to change to run the same simulation with different Monte Carlo programs. That being said not all beam, sample and analyses are supported by all Monte Carlo programs. Supported parameters for each Monte Carlo program are listed in the supported options page.


The program is specific to a particular Monte Carlo program. Each program follows the contract specified by the base program class Program. One implementation is for Casino 2 as part of the package pymontecarlo-casino2. The program can be imported as follow:

[ ]:
from pymontecarlo_casino2.program import Casino2Program

The parameters associated with the program will depend on each Monte Carlo program. For Casino 2, the number of trajectories and the models used for the simulation can be specified. Here is an example with the default models and 5000 trajectories:

[ ]:
program = Casino2Program(5000)

Throughout pyMonteCarlo, a parameter can also be set/modified using its attribute inside the class:

[ ]:
program.number_trajectories = 6000

All parameters are completely mutable and are only validated before a simulation starts.


The second category of parameters is the beam. At the moment, three types of beam are implemented/supported:

  • a pencil beam: beam with no diameter,
  • a Gaussian beam: beam where the electrons are randomly distributed following a two dimensional Gaussian distribution, where the diameter is defined as full width at half maximum (FWHM),
  • a cylindrical beam: beam where the electrons are randomly distributed within a cylinder.

All beam implementations must define the energy and type of the incident particles as defined by the base Beam class. The type of incident particle is defined for future expansions, since all currently supported Monte Carlo programs only accept ELECTRON. Unless otherwise stated, all beams assume that the incident particles travel downwards along the z-axis, i.e. following the vector (0, 0, -1).

The pencil beam is the most supported by the different Monte Carlo programs as no diameter is defined. Here is an example of a pencil beam with a beam energy of 15keV:

[ ]:
from pymontecarlo.options.beam import PencilBeam
beam = PencilBeam(15e3)

Other parameters of the beam are the beam center position. By default, the beam is centered at x = 0m and y = 0m. The position can be changed using either attribute:

[ ]:
beam.x_m = 100e-9
beam.y_m = 200e-9


The sample parameter defines the geometry and the materials of the sample being bombarded by the incident particles. There are currently 5 types of sample implemented:

  • substrate (SubstrateSample): An infinitely thick sample.
  • inclusion (InclusionSample): An half-sphere inclusion in a substrate.
  • horizontal layered (HorizontalLayerSample): Creates a multi-layers geometry. The layers are assumed to be in the x-y plane (normal parallel to z) at tilt of 0.0°.
  • vertical layered (VericalLayerSample): Creates a grain boundaries sample. It consists of 0 or many layers in the y-z plane (normal parallel to x) simulating interfaces between different materials. If no layer is defined, the geometry is a couple.
  • sphere (SphereSample): A sphere in vacuum.

For all types of sample, the sample is entirely located below the z = 0 plane. While some Monte Carlo programs support custom and complex sample definitions, it was chosen for simplicity and compatibility to constrain the available types of sample. If you would like to suggest/contribute another type of sample, please open an enhancement issue or submit a pull request.

Before creating a sample, material(s) must be defined. A material defines the composition and density in a part of the sample (e.g. layer or substrate). After importing the Material class,

[ ]:
from pymontecarlo.options.material import Material

There are three ways to create a material:

  1. Pure, single element material:
[ ]:
material = Material.pure(14) # pure silicon
  1. A chemical formula:
[ ]:
material = Material.from_formula('SiO2')
  1. Composition in mass fraction. The composition is expressed as a dict where keys are atomic numbers and values, mass fractions:
[ ]:
composition = {29: 0.4, 30: 0.6}
material = Material('Brass', composition)

In all three cases the mass density (in kg/m3) can be specified as an argument or set from its attribute:

[ ]:
material.density_kg_per_m3 = 8400
material.density_g_per_cm3 = 8.4

If the density is not specified, it is calculated using this following formula:

\[\frac{1}{\rho} = \sum{\frac{m_i}{\rho_i}}\]

where \(\rho_i\) and \(m_i\) are respectively the elemental mass density and mass fraction of element \(i\).

Each sample has different methods and variables to setup the materials. Here is an example for the substrate sample:

[ ]:
from pymontecarlo.options.sample import SubstrateSample
from pymontecarlo.options.material import Material

copper = Material.pure(29)
substrate = SubstrateSample(copper)

and here is an example for the horizontal layered sample. The substrate is set to copper and two layers are added on top, forming from top to bottom: 100nm of SiO2, 50nm of brass and then copper:

[ ]:
from pymontecarlo.options.sample import HorizontalLayerSample
from pymontecarlo.options.material import Material

copper = Material.pure(29)
sio2 = Material.from_formula('SiO2')
brass = Material('Brass', {29: 0.4, 30: 0.6})

sample = HorizontalLayerSample(copper)
sample.add_layer(sio2, 100e-9)
sample.add_layer(brass, 50e-9)

One trick to make sure the sample is properly setup is to draw it. pyMonteCarlo uses matplotlib to draw the sample in 2D along the XZ, YZ or XY perspective. Here is an example:

[ ]:
import matplotlib.pyplot as plt
from pymontecarlo.figures.sample import SampleFigure, Perspective

fig, axes = plt.subplots(1, 3, figsize=(10, 3))

samplefig = SampleFigure(sample, beams=[beam])

for ax, perspective in zip(axes, Perspective):
    samplefig.perspective = perspective


The analyses define which results from the Monte Carlo simulation will be processed and stored by pyMonteCarlo. To see a list of the supported analyses, please refer to supported options page.

Here is how to store the X-ray intensityies emitted from the sample. First we need to define a photon detector. Each detector requires a name, and the photon detector, an additional argument specifying its elevation, i.e. the angle between the detector and the XY plane.

[ ]:
from pymontecarlo.options.detector import PhotonDetector
import math
detector = PhotonDetector(name='detector1', elevation_rad=math.radians(40))

The photon detector is then used to create a new analysis.

[ ]:
from pymontecarlo.options.analysis import PhotonIntensityAnalysis
analysis = PhotonIntensityAnalysis(detector)


The final step is to put together the program, beam, sample and analysis and create an Options. Note that the options can take several analyses, but in this example we only specified one.

[ ]:
from pymontecarlo.options import Options
options = Options(program, beam, sample, [analysis])

Running simulation(s)

We are now ready to run the simulation. pyMonteCarlo provides an helper function to run several simulation options. These options and their results are automatically stored in a Project object, which can be stored on disk and viewed in either the HDFViewer or the graphical interface of pyMonteCarlo. The results can also be processed programatically, as it will be demonstrated in this tutorial.

[ ]:
from pymontecarlo.runner.helper import run_async
project = await run_async([options])

import os
import tempfile
project.write(os.path.join(tempfile.gettempdir(), 'project1.h5'))

Interpreting simulation results

Let’s now explore the results. Each simulation gets stored in the simulations attribute of the project object:

[ ]:

Each simulation consists in the options used to setup the simulation and the results, which is a list of result objects.

[ ]:
simulation = project.simulations[0]
print('Simulation at {}keV contains {} result(s)'.format(simulation.options.beam.energy_keV, len(simulation.results)))

The PhotonIntensityAnalysis returns an EmittedPhotonIntensityResult, which essentially consists of a dictionary, where the keys are the emitted X-ray lines and the values, their intensities. Here is a quick way to list all X-ray lines. The attribute siegbahn can be replaced with iupac if this notation is preferred. As shown below, the total X-ray intensity of a family of lines (e.g. K, L) is automatically calculated.

[ ]:
result = simulation.results[0]
print('Available X-ray intensities:')
for xrayline in result:

And here is how to retrieve the intensity of one line.

[ ]:
import pyxray
xrayline = pyxray.xray_line('Si', 'Ka')
intensity = result[xrayline]
print('X-ray intensity of {}: {} +/- {}'.format(xrayline.siegbahn, intensity.n, intensity.s))

Another way to analyze the result is to convert them to a pandas dataframe. The project has two methods to create a data frame for the options using create_options_dataframe(...) and for the results create_results_dataframe(...). Each row in these data frames corresponds to one simulation. Both have one required settings argument to specify the X-ray notation and units used.

[ ]:
from pymontecarlo.settings import Settings, XrayNotation
settings = Settings()
settings.preferred_xray_notation = XrayNotation.SIEGBAHN
[ ]:
[ ]:
project.create_results_dataframe(settings, abbreviate_name=True)

This concludes the first tutorial on how to run a single simulation.

[ ]: