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

Run several simulations

One advantage of pyMonteCarlo is the possibility to run several simulations. In this example, we will try to study the influence of the beam energy and the thickness of the carbon coating on the k-ratios of O K\(\alpha\) and Al K\(\alpha\) in orthoclase.

We start by defining the beam energies and carbon thicknesses.

beam_energies_keV = [5.0, 10.0, 20.0]
carbon_thicknesses_nm = [10, 20, 50]

We import the required classes of pyMonteCarlo and Python’s built-in module, itertools. For this example, we will run the simulations with Casino 2, but this could be replaced by any support programs.

import math
import itertools

from pymontecarlo_casino2.program import Casino2Program
from pymontecarlo.options.beam import PencilBeam
from pymontecarlo.options.material import Material
from pymontecarlo.options.sample import HorizontalLayerSample
from pymontecarlo.options.detector import PhotonDetector
from pymontecarlo.options.analysis import KRatioAnalysis
from pymontecarlo.options import Options
from pymontecarlo.runner.helper import run_async

We create the materials.

material_orthoclase = Material.from_formula('KAlSi3O8')
material_carbon = Material.pure(6)

Using itertools.product(...), we create options for every combination of beam energy and carbon thickness and store them in a list.

list_options = []
for beam_energy_keV, carbon_thickness_nm in itertools.product(beam_energies_keV, carbon_thicknesses_nm):
    program = Casino2Program(number_trajectories=5000)

    beam = PencilBeam(beam_energy_keV * 1e3) # Convert to eV

    sample = HorizontalLayerSample(material_orthoclase)
    sample.add_layer(material_carbon, carbon_thickness_nm * 1e-9) # Convert thickness to meters

    detector = PhotonDetector(name='detector1', elevation_rad=math.radians(40))

    analysis = KRatioAnalysis(detector)

    options = Options(program, beam, sample, [analysis])
    list_options.append(options)

With 3 beam energies and 3 carbon thicknesses, we get 9 options.

len(list_options)

We can now run these simulations in parallel. By default, the number of simulations that will be run concurrently depends on the number of CPUs. This can be also modified by the argument, max_workers of the run_async(...) function.

project = await run_async(list_options)

We first checks the number of simulations. We should have 9 simulations from the 9 options provided, but also one simulation for each standard used to calculate the k-ratios. There are 4 elements in orthoclase and one element in the coating. The same standard can be reused for the different carbon coating thicknesses, but not for the different beam energies. So the expected number of simulations should be: 9 + 3 * (4 + 1) = 24. This example shows that pyMonteCarlo is aware which additional simulations it needs to calculate the k-ratios and only simulate these once.

len(project.simulations)

Now to the analysis of the results. There are several ways to extract results from the simulations, but perhaps the easiest one is to convert all the results into a pandas DataFrame. To reduce the number of columns of the DataFrame, it is useful to pass the result class, KRatioReslt in this case, and request only the columns with different information.

from pymontecarlo.results import KRatioResult
from pymontecarlo.settings import Settings, XrayNotation

settings = Settings()
settings.set_preferred_unit('keV')
settings.set_preferred_unit('nm')
settings.set_preferred_unit('deg')
settings.set_preferred_unit('g/cm^3')
settings.preferred_xray_notation = XrayNotation.SIEGBAHN

df = project.create_dataframe(settings, result_classes=[KRatioResult], only_different_columns=True)
df

The pandas DataFrame above contains one row for each simulation, including the standards. Since we are only interested in the k-ratios of the unknowns, we can filter the DataFrame and drop some columns.

# Keep only standard rows
df = df[df['standard'] != True].dropna(axis=1)

# Remove columns which all contain the same values
df = df[[col for col in df if df[col].nunique() != 1]]
df

We can then use this DataFrame to plot the results using matplotlib.

import matplotlib.pyplot as plt

xraylines = ['O Kα', 'Al Kα']

fig, axes = plt.subplots(1, len(xraylines), figsize=(5 * len(xraylines), 4))
fig.subplots_adjust(wspace=0.3)

for ax, xrayline in zip(axes, xraylines):
    for beam_energy_keV, df_beam_energy in df.groupby('beam energy [keV]'):
        ax.plot(df_beam_energy['layer #0 thickness [nm]'], df_beam_energy[xrayline], 'o-', label=f'E0={beam_energy_keV:.0f} keV')

    ax.set_xlabel('Carbon coating thickness (nm)')
    ax.set_ylabel(f'{xrayline} k-ratio')

axes[0].legend(loc='best')

plt.show()