docs/tutorials/first-simulation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# First simulation\n",
    "\n",
    "A (Monte Carlo) Simulation consists in (1) [options](#Setting-up-simulation-options), defining all the necessary parameters to setup the simulation, and (2) [results](#Interpreting-simulation-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](https://www.hdfgroup.org/) file and can be opened in the [HDFViewer](https://www.hdfgroup.org/) or using any HDF5 library."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setting up simulation options\n",
    "\n",
    "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: \n",
    "\n",
    "* program\n",
    "* beam\n",
    "* sample\n",
    "* analyses\n",
    "\n",
    "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](../supported-options.rst) page."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Program\n",
    "\n",
    "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](http://www.gel.usherbrooke.ca/casino) as part of the package **pymontecarlo-casino2**. The program can be imported as follow:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pymontecarlo_casino2.program import Casino2Program"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "program = Casino2Program(5000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Throughout **pyMonteCarlo**, a parameter can also be set/modified using its attribute inside the class:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "program.number_trajectories = 6000"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "All parameters are completely mutable and are only validated before a simulation starts."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Beam\n",
    "\n",
    "The second category of parameters is the beam. At the moment, three types of beam are implemented/supported: \n",
    "\n",
    "* a pencil beam: beam with no diameter, \n",
    "* 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),\n",
    "* a cylindrical beam: beam where the electrons are randomly distributed within a cylinder.\n",
    "\n",
    "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)`.\n",
    "\n",
    "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:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pymontecarlo.options.beam import PencilBeam\n",
    "beam = PencilBeam(15e3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "beam.x_m = 100e-9\n",
    "beam.y_m = 200e-9"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sample\n",
    "\n",
    "The sample parameter defines the geometry and the materials of the sample \n",
    "being bombarded by the incident particles.\n",
    "There are currently 5 types of sample implemented:\n",
    "\n",
    "* substrate (`SubstrateSample`): An infinitely thick sample. \n",
    "* inclusion (`InclusionSample`): An half-sphere inclusion in a substrate.\n",
    "* horizontal layered (`HorizontalLayerSample`): Creates a multi-layers geometry. \n",
    "  The layers are assumed to be in the x-y plane (normal parallel to z) at tilt of 0.0°.\n",
    "* vertical layered (`VericalLayerSample`): Creates a grain boundaries sample.\n",
    "  It consists of 0 or many layers in the y-z plane (normal parallel to x) simulating interfaces between different materials.\n",
    "  If no layer is defined, the geometry is a couple.\n",
    "* sphere (`SphereSample`): A sphere in vacuum.\n",
    "    \n",
    "For all types of sample, the sample is entirely located below the ``z = 0`` plane.\n",
    "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.\n",
    "If you would like to suggest/contribute another type of sample, please open an enhancement [issue](https://github.com/pymontecarlo/pymontecarlo/issues) or submit a [pull request](https://github.com/pymontecarlo/pymontecarlo/pulls).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before creating a sample, material(s) must be defined.\n",
    "A material defines the composition and density in a part of the sample (e.g. layer or substrate).\n",
    "After importing the `Material` class, "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pymontecarlo.options.material import Material"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are three ways to create a material:\n",
    "\n",
    "1. Pure, single element material:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "material = Material.pure(14) # pure silicon"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. A chemical formula:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "material = Material.from_formula('SiO2')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3. Composition in mass fraction. \n",
    "   The composition is expressed as a *dict* where keys are atomic numbers and values, mass fractions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "composition = {29: 0.4, 30: 0.6}\n",
    "material = Material('Brass', composition)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In all three cases the mass density (in kg/m3) can be specified as an argument or set from its attribute:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "material.density_kg_per_m3 = 8400\n",
    "material.density_g_per_cm3 = 8.4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If the density is not specified, it is calculated using this following formula:\n",
    "\n",
    "$$\\frac{1}{\\rho} = \\sum{\\frac{m_i}{\\rho_i}}$$\n",
    "\n",
    "where $\\rho_i$ and $m_i$ are respectively the elemental mass density and mass fraction of element $i$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Each sample has different methods and variables to setup the materials. \n",
    "Here is an example for the substrate sample:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pymontecarlo.options.sample import SubstrateSample\n",
    "from pymontecarlo.options.material import Material\n",
    "\n",
    "copper = Material.pure(29)\n",
    "substrate = SubstrateSample(copper)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and here is an example for the horizontal layered sample. \n",
    "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:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pymontecarlo.options.sample import HorizontalLayerSample\n",
    "from pymontecarlo.options.material import Material\n",
    "\n",
    "copper = Material.pure(29)\n",
    "sio2 = Material.from_formula('SiO2')\n",
    "brass = Material('Brass', {29: 0.4, 30: 0.6})\n",
    "\n",
    "sample = HorizontalLayerSample(copper)\n",
    "sample.add_layer(sio2, 100e-9)\n",
    "sample.add_layer(brass, 50e-9)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One trick to make sure the sample is properly setup is to draw it.\n",
    "**pyMonteCarlo** uses [matplotlib](http://matplotlib.org) to draw the sample in 2D along the XZ, YZ or XY perspective.\n",
    "Here is an example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from pymontecarlo.figures.sample import SampleFigure, Perspective\n",
    "\n",
    "fig, axes = plt.subplots(1, 3, figsize=(10, 3))\n",
    "\n",
    "samplefig = SampleFigure(sample, beams=[beam])\n",
    "\n",
    "for ax, perspective in zip(axes, Perspective):\n",
    "    samplefig.perspective = perspective\n",
    "    samplefig.draw(ax)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Analyses\n",
    "\n",
    "The analyses define which results from the Monte Carlo simulation will be processed and stored by **pyMonteCarlo**. \n",
    "To see a list of the supported analyses, please refer to [supported options](../supported-options.rst) page."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pymontecarlo.options.detector import PhotonDetector\n",
    "import math\n",
    "detector = PhotonDetector(name='detector1', elevation_rad=math.radians(40))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The photon detector is then used to create a new analysis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pymontecarlo.options.analysis import PhotonIntensityAnalysis\n",
    "analysis = PhotonIntensityAnalysis(detector)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Options\n",
    "\n",
    "The final step is to put together the program, beam, sample and analysis and create an `Options`. \n",
    "Note that the options can take several analyses, but in this example we only specified one."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pymontecarlo.options import Options\n",
    "options = Options(program, beam, sample, [analysis])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Running simulation(s)\n",
    "\n",
    "We are now ready to run the simulation. \n",
    "**pyMonteCarlo** provides an helper function to run several simulation options.\n",
    "These options and their results are automatically stored in a `Project` object, which can be stored on disk and viewed in either the [HDFViewer](https://www.hdfgroup.org/) or the graphical interface of **pyMonteCarlo**. \n",
    "The results can also be processed programatically, as it will be demonstrated in this tutorial."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pymontecarlo.runner.helper import run_async\n",
    "project = await run_async([options])\n",
    "\n",
    "import os\n",
    "import tempfile\n",
    "project.write(os.path.join(tempfile.gettempdir(), 'project1.h5'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Interpreting simulation results\n",
    "\n",
    "Let's now explore the results. \n",
    "Each simulation gets stored in the `simulations` attribute of the `project` object:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "project.simulations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Each simulation consists in the `options` used to setup the simulation and the `results`, which is a list of result objects."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "simulation = project.simulations[0]\n",
    "print('Simulation at {}keV contains {} result(s)'.format(simulation.options.beam.energy_keV, len(simulation.results)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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.\n",
    "Here is a quick way to list all X-ray lines.\n",
    "The attribute `siegbahn` can be replaced with `iupac` if this notation is preferred.\n",
    "As shown below, the total X-ray intensity of a family of lines (e.g. K, L) is automatically calculated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "result = simulation.results[0]\n",
    "print('Available X-ray intensities:')\n",
    "for xrayline in result:\n",
    "    print(xrayline.siegbahn)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And here is how to retrieve the intensity of one line."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pyxray\n",
    "xrayline = pyxray.xray_line('Si', 'Ka')\n",
    "intensity = result[xrayline]\n",
    "print('X-ray intensity of {}: {} +/- {}'.format(xrayline.siegbahn, intensity.n, intensity.s))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another way to analyze the result is to convert them to a [pandas](https://pandas.pydata.org/) dataframe.\n",
    "The project has two methods to create a data frame for the options using `create_options_dataframe(...)` and for the results `create_results_dataframe(...)`.\n",
    "Each row in these data frames corresponds to one simulation.\n",
    "Both have one required `settings` argument to specify the X-ray notation and units used."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pymontecarlo.settings import Settings, XrayNotation\n",
    "settings = Settings()\n",
    "settings.set_preferred_unit('eV')\n",
    "settings.set_preferred_unit('nm')\n",
    "settings.set_preferred_unit('deg')\n",
    "settings.preferred_xray_notation = XrayNotation.SIEGBAHN"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "project.create_options_dataframe(settings)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "project.create_results_dataframe(settings, abbreviate_name=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This concludes the first tutorial on how to run a single simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}